Measuring and crust-correcting finite-frequency travel time residuals – application to southwestern Scandinavia

We present a data-processing routine to compute relative finite-frequency travel time residuals using a combination of the Iterative Cross-Correlation and Stack (ICCS) algorithm and the Multi-Channel Cross-Correlation method (MCCC). The routine has been tailored for robust measurement of Pand S-wave travel times in several frequency bands and for avoiding cycle-skipping problems at the shortest periods. We also investigate the adequacy of ray theory to calculate crustal corrections for finite-frequency regional tomography in normal continental settings with nonthinned crust. We find that ray theory is valid for both P and S waves at all relevant frequencies as long as the crust does not contain low-velocity layers associated with sediments at the surface. Reverberations in the sediments perturb the arrival times of the S waves and the long-period P waves significantly, and need to be accounted for in crustal corrections. The data-processing routine and crustal corrections are illustrated using data from a network in southwestern Scandinavia.


Introduction
Since the pioneering work of Aki et al. (1977), regional tomographic models have had a considerable impact on our understanding of the upper mantle.A tremendous increase in data quantity and improvements in inversion techniques have greatly enhanced tomographic resolution power in recent years.The quality of the models is however fully dependent on the quality of the travel time residuals the models are based upon: better measurements will ultimately translate to better models, and good routines for travel time measurements are therefore crucial.
Relative travel times, in most cases defined as relative to the mean arrival time of each event, or in some cases relative to a reference station, are usually measured using a crosscorrelation method.These methods exploit the similarity of seismic waveforms measured at different stations to estimate precise travel times in a seismic network (e.g.Bungum and Husebye, 1971;VanDecar and Crosson, 1990;Sigloch and Nolet, 2006;Pavlis and Vernon, 2010).Inherent in all cross-correlation methods is a risk of cycle skipping during analysis (e.g.VanDecar and Crosson, 1990;Chevrot, 2002), meaning that an earlier or later wiggle in the seismogram is interpreted as the arrival of the phase of interest.Initial alignment of traces is therefore necessary such that a narrow window can be used in the search for the maximum of the cross-correlation function.Predictions from theoretical calculations of arrival times (e.g.Crotwell et al., 1999) are usually not sufficient to solve this issue (e.g.Lou et al., 2013), as the travel time residuals are commonly larger than or similar to wave periods in the highest frequency range.Picking of approximate phase arrival times is therefore usually performed manually before cross-correlation analysis.Lou et al. (2013) recently proposed an alignment procedure to avoid manual picking.This Iterative Cross-Correlation and Stack algorithm (ICCS) and other recent tools (e.g.Chevrot, 2002;Rawlinson and Kennett, 2004;Pavlis and Vernon, 2010), are however tailored for measuring high-frequency signals.With the development of finitefrequency theory (Dahlen et al., 2000;Hung et al., 2000), it has become essential to measure travel times in several narrow frequency bands.We find from experience that the measurement of travel times in more narrow high-frequency bands increases the risk of cycle skipping during analysis in both the Lou et al. (2013) and VanDecar and Crosson (1990) algorithms, especially for lower quality events from great epicentral distances ( > 75 • ).An attempt to remedy this problem in finite-frequency measurements of P-wave arrivals has been presented by Bonnin et al. (2014).By introducing a weak constraint that arrival times should be close to the theoretical ones in a reference model, they eliminate the data with, among other problems, cycle skipping.This might however discard extreme but good data.We therefore rather aim at correcting and using the data suffering from cycleskipping problems.
We propose here a very simple approach that manages to remedy this problem to a large extent.Inspired by waveform inversion of surface waves (Cara and Lévêque, 1987), we measure arrival times in the lowest frequency band first and then in successively higher bands, after adjusting for the arrival times found in the previous frequency band.We apply the ICCS algorithm followed by the classical Multi-Channel Cross-Correlation method (MCCC) of VanDecar and Crosson (1990) along with several rounds of automated data rejection that increase the robustness of the initial stack in the ICCS algorithm.The application of the MCCC ensures robust results and provides quantitative uncertainty estimates important in the further tomographic inversion.
In ACH-type body-wave tomography (ACH refers to the authors in Aki et al. (1977)), the structure of the crust is usually not inverted for, but its influence on the travel times is taken into account by a crustal correction applied to each station individually (e.g.Allen et al., 2002;Martin et al., 2005).Crustal corrections are usually computed using ray theory (e.g.Tian et al., 2007), providing frequency-independent results.It is well-established that the reverberations in the crust affect the travel times at different frequencies in different ways, and that crustal delay times of high-frequency and lowfrequency data are significantly different (Obayashi et al., 2004;Yang and Shen, 2006;Ritsema et al., 2009).As finitefrequency tomography plays on the different sensitivity of high-and low-frequency signals to the structure, it is essential that the differences between the high-frequency and lowfrequency residuals are not biased by inadequate crustal corrections.Obayashi et al. (2004) and Yang and Shen (2006) have shown the importance of this frequency-dependence for P waves with applications mainly in thin oceanic crustal environments.We analyse it here for the continental crust and in the framework of relative travel time tomography, where only relative residuals within the network are important.We analyse the discrepancy between the crustal corrections calculated with ray theory and those calculated with a reflectivity algorithm that takes into account the frequency-dependent effect of the crust on the travel times.We find that the discrepancy is particularly significant in the presence of sedimentary layers and that, depending on the geological context and used frequency range, one should always consider the need to apply frequency-dependent corrections.
Our residual measurement procedure and the frequencydependence of the crustal corrections are illustrated using data from a network in southwestern Scandinavia (Weidle et al., 2010).

Measurement of travel time residuals
We apply our analysis to seismological data recorded by the temporary MAGNUS network between September 2006 and June 2008 (Weidle et al., 2010), by the temporary DANSEIS network from April 2008 to June 2008, by the temporary CALAS stations (Medhus et al., 2009) and by permanent stations in the study area (Fig. 1).Direct P and S phases were studied from earthquakes with epicentral distances of 30-95 • (Fig. 2).To have an even azimuthal coverage, we used events with smaller magnitude in the southern and western quadrant, decreasing the data quality and increasing the need for a robust algorithm for travel time measurement.In short, the data-processing sequence to measure residuals is performed in four steps.A major step to make the procedure robust to cycle skipping is to perform the processing first on the lowest frequency band, and then on successively higher frequency bands.Corrections for ellipticity and crustal structure are calculated after step 4 and will be presented in the next section.We use data from a 6.4 M event at the Andreanoff Islands ( = 68 • , back azimuth 5 • ), 15 April 2008, to illustrate the data processing.

Step 1: preprocessing
This step, which is detailed in Kolstrup (2015), consists of filtering, calculation of theoretical phase arrival time, and windowing according to the frequency content of the band and the phase in question.Frequency bands for P waves are 0.03-0.125Hz (33-8 s) and 0.5-2 Hz (2-0.5 s) to isolate the secondary noise peak in southern Norway around 0.2 Hz (5 s).S waves are bandpass-filtered in the ranges 0.03-0.077Hz (33-13 s) and 0.077-0.125Hz (13-8 s).These frequency ranges are similar to those used generally in finite-frequency tomography (Hung et al., 2004;Schmandt and Humphreys, 2010;Hung et al., 2011), with the exception that S waves are usually measured up to slightly higher frequencies (up to 0.5 Hz in Hung et al., 2004).Schmandt and Humphreys (2010) note that their highest frequency for S waves, 0.4 Hz, produces few data.

Step 2: data rejection based on the envelope of data
The first round of data rejection excludes traces based on an iterative comparison of their envelopes with the envelope of the stack (see Kolstrup, 2015 for details).This ensures that strong outliers do not ruin the computation of the array stack before applying the ICCS algorithm in the next step.

Step 3: ICCS algorithm tailored for several frequency bands
To avoid cycle skips in the Multi-Channel Cross-Correlation analysis (MCCC, step 4), it is essential that traces are initially aligned, either by manual picking or by an automatic procedure (VanDecar and Crosson, 1990;Lou et al., 2013), and that noisy or incoherent traces are excluded from the analysis.We use the Iterative Cross-Correlation and Stack algorithm (ICCS, Lou et al., 2013) to initially align traces, but we do not pick the absolute arrival time of the stack as in Lou et al. (2013), since our focus is relative travel times.The ICCS algorithm (Lou et al., 2013) was recently developed to replace the unavoidable initial phase-marking part of the cross-correlation procedures, and is part of the AIMBAT tool (Automated and Interactive Measurement of Body-wave Arrival Time) (Lou et al., 2013).In the ICCS algorithm, each individual trace is correlated with the stack and the time lag associated with the maximum of the cross-correlation function is found.Then each individual trace is aligned according to the time lag and relative to the stack in an iterative  procedure where the stack is updated for each iteration.This iterative alignment procedure is very similar to the adaptive stacking in Rawlinson and Kennett (2004), except that Rawlinson and Kennett (2004) use an L3 misfit criterion to determine the time lags of the traces with respect to the stack instead of the maximum of the cross-correlation function.
Figure 3 shows the high-frequency P-wave traces before and after alignment by ICCS.
After alignment by ICCS, the cross-correlation coefficient and mean spectral coherence between each trace and the stack are calculated.A weighted average of the two is computed and traces with a value lower than a user-defined cutoff (usually about 0.5) are excluded.This procedure rejects data with a significantly different shape than the array stack.Figure 4 shows the low-and high-frequency P-wave stacks after alignment by the ICCS algorithm with a few greencoloured incoherent traces that will be rejected from further analysis.
Cycle skipping is a recurrent problem in cross-correlation analysis when the period of the data and the travel time residuals are of the same order of magnitude.For the lowfrequency bands of both the P and S waves, the travel time residuals are much smaller than the periods of the waves, and the risk of cycle skipping in the pick of cross-correlation maxima is therefore very low.In contrast, there is a significant risk of cycle skipping in the ICCS algorithm in case of high-frequency low-quality data.This mostly happens in the first iterations where the array stack is less robust, because it is calculated from traces aligned according to the theoretical arrival time.With a poor stack there will be little difference between local and global maxima in the crosscorrelation functions, and the algorithm might therefore misplace a few traces, possibly permanently.
This problem of cycle skipping in the ICCS algorithm is probably enhanced in our data set because most events have epicentral distances > 70 • , giving relatively strong attenuation of the high-frequency content of the data.For events with shorter epicentral distances, even with magnitudes M < 6.0, the problem of cycle skipping in the ICCS algorithm is almost non-existent.
To remedy the risk of cycle skipping in the first iterations of the ICCS algorithm for the higher frequency bands, we use the low-frequency travel time residuals from the MCCC analysis in step 4 as input time lag to initially align the traces in the stack.Only reliable travel time residuals with small standard deviations are used as input time lags.This ensures a better initial stack in the high-frequency bands and reduces the number of iterations in the ICCS algorithm.
The effect of the improved initial stack is especially important for the more band-limited S-wave data.Figure 5 shows the improvement of the initial stack, when the low-frequency S-wave residuals are used to align the traces of the highfrequency band, compared to using the theoretical arrival time from AK135.

Step 4: Multi-Channel Cross-Correlation (MCCC) analysis
Finally, we use the classical Multi-Channel Cross-Correlation method (MCCC) of VanDecar and Crosson (1990) to measure the relative travel times t i and evaluate their uncertainty.
The MCCC algorithm is a robust method for relative travel time measurement that also provides quantitative uncertainty estimates.In the MCCC method, relative delay times t ij between all pairs of stations i, j are measured using the cross-correlation of seismograms.These relative delays are then inverted in a least-squares sense for the optimized relative arrival times, t i , at each station, under the constraint of zero mean, i.e. t i = 0.The equation residuals between ob- served relative delays and predicted delays are calculated as res ij = t ij −(t i −t j ).The standard deviation of the equation residuals associated with each trace, σ i , provides an estimate of the timing uncertainty associated with the relative arrival time t i (VanDecar and Crosson, 1990).An evaluation of the uncertainty is particularly important in order to be able to invert P-and S-waves arrival times in a consistent manner and interpret jointly the obtained P-and S-wave models.
If the mean value is subtracted from the final time lags determined by the ICCS algorithm, they are approximately equal to the relative travel time residuals calculated by the MCCC algorithm (Lou et al., 2013).The main reason for applying MCCC to the traces aligned by ICCS is therefore the calculation of uncertainty estimates, σ i , from the equation residuals res ij .
Working with high-frequency P waves, VanDecar and Crosson (1990) found that equation residuals res ij > 0.5 s are usually the result of cycle skipping in the crosscorrelation analysis.For the low-frequency P waves we find that a threshold of 0.8 s is appropriate to target the (very few) cycle skips, and for the S waves we use larger thresholds of 2 s for the low-frequency band and 1.5 s for the highfrequency band.The cross-correlation functions associated with large equation residuals are recalculated in a smaller window near the time lag predicted by the least-squares solution (t i − t j ), as in VanDecar and Crosson (1990).The data set of relative delay times ( t ij ) is updated with the new measurements before a weighted reinversion using the crosscorrelation coefficients between traces as weights.The dif-ferences between the unweighted relative travel times and the weighted relative travel times are very small.

Crustal corrections
Before performing a tomographic inversion, the residuals need some corrections.We compute relative travel time corrections for the ellipticity of the Earth (Kennett and Gudmundsson, 1996) using the MATLAB implementation of Euler (2014).These corrections are in general not very important as they are very small compared to the residuals (−0.05 to 0.09 s for P waves and −0.10 to 0.14 s for S waves in our case).
Crustal corrections are much more important (e.g.Allen et al., 2002;Martin et al., 2005).They are usually computed using ray theory (e.g.Tian et al., 2007), providing frequency-independent results.It is well-established that the reverberations in the crust affect the travel times at different frequencies in different ways, and that crustal delay times of high-frequency and low-frequency data are significantly different (Obayashi et al., 2004;Yang and Shen, 2006;Ritsema et al., 2009).Obayashi et al. (2004) shows in particular that the analysis of PP waves reflecting at oceanic locations with a thin crust requires a frequency-dependent correction and Yang and Shen (2006) shows the importance of the frequency-dependence for P waves, in particular for a rather thin oceanic crust.In the case of regional body-wave tomographies, the absolute difference between the travel times of the short-and long-period waves is not an issue since the average travel time, calculated separately in each frequency band, is subtracted.The important element here is if the frequency-dependence differs from one station to the other.

Frequency-dependent crustal corrections
In order to verify the adequacy of ray-theory corrections, we have computed the delays of long-period and short-period waves across the crustal structures in our study area using both ray theory and the reflectivity method.The reflectivity method assumes a 1-D model below each station and does not take into account lateral variations, but does take into account all reverberations and conversions within the stack of layers.It generates the three-component impulse response of the layered structure, assuming a plane wave incident from below with a prescribed incidence angle.We computed the response of the crustal layer for P and SH waves using the reflectivity software of Levin and Park (1997) and defined the frequency-dependent crustal travel times as the time of the maximum of the impulse response after filtering with the same zero-phase bandpass filter as applied to the data.In order to ensure a common reference time and include the effect of the topography, the crustal corrections are computed at each station by calculating the travel times from 50 km below sea level and up to the free surface, for waves with the prescribed incidence.As the relative travel time residuals are zero mean, final corrections also have to be zero mean and the actual corrections applied for each event are the travel times from 50 km depth minus the average correction for the selected data for the particular event and frequency range.In order to keep a common reference, we will however show here crustal corrections relative to the average correction calculated with ray theory.
Our study area comprises a deep sedimentary basin above a rather thin crust to the south, typical of younger continental regions.Most of the area to the north is an older cratonic region with a purely crystalline crust of typically 45 km to the east, thinning to smaller values to the west, closer to the Atlantic continental margin.Its lateral variations are thus representative of many settings for regional body-wave tomography.At each station, we use a detailed local crustal model compiled in Kolstrup (2015) from a range of sources, especially the work of Stratford and Thybo (2011); Lassen and Thybo (2012) and Kolstrup and Maupin (2013).The model is specified as a 1-D model below each station.The total thickness of the sedimentary layers and the thickness of the crust, which are the two main controlling factors for crustal corrections, are shown in Fig. 6, and the S-wave velocity down to 50 km depth is shown at a selection of stations in Fig. 7.
Figure 8 shows the impulse responses of the crustal structures for P and SH waves incident at two different stations together with their bandpassed versions.At station NWG12 (left plots), located on crystalline crust, the secondary arrivals associated with the crustal reverberations are late enough not to perturb the location in time of the maximum of the impulse response after filtering.At station DK08 (right plots), located on 5 km of sedimentary layers, the first arrival is followed closely by large arrivals in both the P and SH cases, shifting the arrival time of the maximum in the bandpassed signals significantly.The long-period waves arrive usually earlier than predicted by ray theory.The upper panels of Fig. 9 shows P-wave crustal corrections calculated by ray theory and reflectivity at short and long periods for all stations.The examples shown here are for slownesses corresponding to 75 • epicentral distance, but the same results are found at all other relevant slownesses.The corrections calculated with the reflectivity method are not significantly different from those calculated with ray theory, except at low frequency in Denmark, in the southern part of the study region, where the discrepancy reaches about 0.3 s.This discrepancy is clearly associated with the presence of sediments in the area (Fig. 6).We can notice that the crustal corrections at long periods follow rather closely the thickness of the crust, even in the south.This simply shows, not unexpectedly, that the relatively thin sedimentary layers are not seen by the 100 km wavelength long-period P waves.A good approximation would therefore be to calculate the corrections at long periods with ray theory but replace the sedimentary layers by crystalline crust in the upper part of the model.
Frequency-dependence in the sedimentary region is also apparent for the SH waves (lower panels of Fig. 9).The main difference with the P-wave case is that we notice a discrepancy with ray theory already in the high-frequency range, which has a central frequency of about 10 s here.The socalled high frequencies are therefore also affected by the reverberations in the sedimentary layer.The discrepancy is however highly non-linearly correlated with the sediment thickness: it is very small at the stations to the west in Denmark, where we have the thickest sediments (more than 5 km), but similar to the values of 0.4 to 0.7 s that are relevant for the lower frequencies at the outskirts of the basin, with sediments of 1 to 5 km thickness.In our case, the effect of the reverberations on the SH waves is more complex than for the P waves and cannot easily by modelled using ray theory only.In classical continental environments, and for the frequencies that we have used in our study, the frequencydependence of the crustal variations is not associated with variations in crustal thickness, but with variations in the lowvelocity layers in the upper crust.This dependence is of the same order as the corrections themselves, and cannot usually be neglected.

Residuals before and after crustal corrections
Figure 10 shows high-frequency and low-frequency P-wave residuals before and after crustal correction.The residuals have been measured according to the procedure described in Sect. 2 before crustal correction.In order to have a good station coverage, we do not show residuals from a particular event but average residuals for a cluster of events in back azimuths from 0 to 65 • .Averages in other back azimuths are shown in Kolstrup (2015) and show the same features as those presented here concerning the frequency-dependency of the residuals.The variation of travel time residuals with back azimuths provides preliminary information on the spatial distribution of velocity anomalies prior to inversion, as discussed in Kolstrup (2015).We focus here only on the effect of the crustal corrections on the residuals.
The dominant pattern in the residual variation is not strongly affected by the crustal corrections, showing that the signal in this example study is dominated by mantle heterogeneity.The residuals corrected using ray theory (middle panels) have somewhat reduced small-scale variations  compared to uncorrected residuals, especially at long periods.Comparison of the spatial patterns in the low-and high-frequency residuals show that the raw and ray-corrected residuals are significantly different in the southern part of the study area, where the low-frequency waves arrive early, whereas the high-frequency waves arrive late.After correction with the reflectivity method (right panels), this difference between high-and low-frequency residuals has disappeared.This clearly show that the low frequencies arrive earlier than the high frequencies in Denmark because they are not sensitive to the presence of thick sediments, and this discrepancy should not be interpreted in terms of frequencydependent sensitivity of the travel times to the mantle structure.
For the S waves (Fig. 11), as the two frequency ranges are closer to each other and the crustal corrections are smaller relative to the residuals, the trend in the variation of the residuals from high to low frequencies and from ray-corrected to reflectivity-corrected residuals is not as clear.The residuals at high and low frequency do not get closer to each other in Denmark after reflectivity-based corrections.As noticed already, the reflectivity-based correction in the highest fre-quency range depends on the thickness of the sedimentary layer in a very non-linear way: thin layers of 1 to 5 km lead to discrepancies of 0.4 to 0.7 s between the two types of corrections, whereas thicker layers lead to small discrepancies.Uncertainties in the crustal model or neglecting its lateral variations may therefore lead to errors in the crustal corrections that are of the same order of magnitude as the error related to using ray theory.
We have done the crustal correction after measurement of the residuals.In the present case, it has the advantage of ensuring that the similarity of the residuals in the two frequency bands after reflectivity-based correction is not introduced by our measurement procedure based on measuring low-frequency bands first.In another setting, it would make sense to introduce crustal correction before measuring the residuals, at the same time as the correction for theoretical arrival times, as this would improve the initial alignment for the ICCS algorithm.One should however then also take into account the difference between the low-frequency and the highfrequency corrections when aligning the short-period traces according to the low-frequency residuals.As the frequencydependent discrepancies are smaller than the periods used, −0.5 0 0.5 Figure 10.Average P-wave travel time residuals in the back-azimuthal range 0 to 65 • .Left panels: uncorrected residuals; middle panels: residuals corrected for topography and crustal structure with ray theory; right panels: residuals corrected for topography and crustal structure with reflectivity.The upper panels are for the high-frequency range and the lower panels are for the low-frequency range.
we do not expect a large benefit from this procedure, but just a better coherency.

Discussion and conclusion
A prerequisite for high-quality regional body-wave tomography is high-quality measurement of the residuals.As finitefrequency body-wave tomography requires the measurement of residuals in different frequency bands on bandpassed traces, a number of new issues arise in the measurement and quality assurance of residuals for this kind of tomography.
We have addressed two of these issues here: how to measure the residuals in a robust way, and whether the crustal correction should be made frequency-dependent.

Measuring residuals
We propose an automated processing procedure combining the Iterative Cross-Correlation and Stack (ICCS) and Multi-Channel Cross-Correlation method (MCCC) algorithms and tailored for estimating relative travel time residuals in several frequency bands.Using the low-frequency travel time residuals to initially align the high-frequency traces increases the reliability of the automated alignment in ICCS and reduces the risk of cycle skipping in the MCCC analysis, especially for distant low-quality data.As the only human interference is the choice of various parameters and a quality check at the end, there is little risk of drift in the arrival time picking and a high degree of objectiveness in the data rejection procedures.Combining the ICCS and MCCC algorithms, we benefit from the advantages of the two methods, the ICCS providing a robust way of aligning traces for measuring their travel times and the MCCC providing an estimate of the uncertainties in the measurements.
Using the long-period residuals as a priori information to adjust the high-frequency traces carries the risk of propagating errors due to long-period noise into the high-frequency residuals.To reduce this risk, only low-frequency residuals with small standard deviations (in the MCCC sense) are used to align the high-frequency traces.The alignment at these high-quality stations benefits measurements at all stations as it improves the quality of the stack in the ICCS algorithm.As long as spurious cycle skips are not introduced by bad longperiod data, our procedure of shifting traces in time does not modify the waveform or their cross-correlations, and there-  fore does not reduce the independency of the low-and highfrequency measurements.Medhus et al. (2012) and Wawerzinek et al. (2013) used basically the same database as we do, but employed different methods for measuring travel times.The P-wave residuals of Medhus et al. (2012) were estimated by a combination of cross-correlation measurements and final manual picking of high-frequency P waves (0.125-4 Hz), whereas the S-wave travel time residuals of Wawerzinek et al. (2013) were found by manual picking of S waves on magnified waveforms in a lower frequency band (0.03-0.125 Hz).An advantage of using the same cross-correlation procedure on both the Pand S-wave data set, is that uncertainty estimates are calculated in a consistent manner for both data sets, facilitating a quantitative comparison of the data sets and their joint inversion, as we have done in Kolstrup et al. (2015) with the data set used here.In a later inversion for seismic velocities, these uncertainties can be used to weight the data objectively, suppressing the influence of individual noisy data and lowquality events.The two data sets can be inverted either jointly for e.g.V and V S , or for V P and V S separately.If the same inversion technique and model parameterization is used to invert for V P and V S , the resulting models can be compared quantitatively as in e.g.Chou et al. (2009) and Hung et al. (2011).

Crust-correcting residuals
Another factor contributing to the quality of the residuals is how the influence of the non-resolved upper part of the Earth is dealt with.Station correction terms can be used to absorb the effect of the crust.In order to avoid absorbing a significant part of the mantle heterogeneity as well, these terms need to be damped.If the structure of the crust is known well enough, it is therefore an advantage to correct for the crust even if a station term is used in the inversion.Although it is well-established that the crust affects travel times in a frequency-dependent manner (Obayashi et al., 2004;Yang and Shen, 2006;Ritsema et al., 2009), crustal corrections in regional body-wave tomography are usually computed independently of frequency, using ray theory (Hung et al., 2004;Schmandt and Humphreys, 2010;Kolstrup et al., 2015).As finite-frequency tomography plays on the different sensitivities of different frequencies to the Earth's structure, and as the differences in the residuals at different frequencies are usually not very large (35 % of the short-period residuals for P waves in Hung et al. (2004) and similarly in Hung et al. (2011), less for S waves in Hung et al. (2011) and in the present study), biases in the crustal corrections may harm the finite-frequency approach or at least limit its usefulness.
We have analysed, in a typical continental setting, the discrepancy between crustal corrections calculated with ray theory and those calculated with the reflectivity method, which takes into account all crustal reverberations.In our continental model of southern Scandinavia, where crustal thickness varies from 27 to 48 km, we find that Moho depth variations do not produce any significant discrepancy between the two types of corrections.This is true for both P and S waves, in any of the rather conventional frequency ranges that we have used.This is in agreement with Obayashi et al. (2004) and Yang and Shen (2006) who show that significant effects occur for thinner crust, e.g.oceanic crust.
The presence of 1 to 8 km thick sedimentary layers in the upper part of the crust has, on the other hand, large effects.We find discrepancies with ray theory of about 0.3 s for long-period P waves and 0.4 to 0.7 s for S waves.For P waves, a valid approach in our case is to consider that the short-period waves are sensitive to the sedimentary layers but the long periods "do not see them" and calculate the crustal corrections with ray theory, but using two different models: one including the sediments, the other one with sediments replaced by upper crustal rocks.This approach may not be usable if an intermediate frequency band is introduced, as in Hung et al. (2004).It is not usable for S waves either, which are usually measured in frequency ranges such that they interfere in a more complex way with the sedimentary layers.For sedimentary layers of 1 to 5 km, the S waves arrive early by 0.4 to 0.7 s compared to ray-theory predictions in both frequency ranges 0.03-0.077Hz (33-13 s) and 0.077-0.125Hz (13-8 s).The shortest periods are consistent with ray theory for larger sedimentary thicknesses, but the longest periods are not.This complex and non-linear behaviour makes it a challenge to compute accurate crustal corrections for S waves: uncertainties in the thickness of the sedimentary layers or 3-D effects not accounted for here may lead to large inadequacies in these corrections.Frequency-dependent station terms may in some cases help absorb some of the potential errors introduced by ray theory or by uncertainty in the model.
Although our main purpose is to analyse the crustal corrections required in finite-frequency tomography, our results are also relevant for traditional regional body waves, especially in the S-wave case.Due to their lower frequency content and the higher noise level on the horizontal components, S waves are usually filtered to rather low frequency in traditional body-wave tomography.Wawerzinek et al. (2013), for example, use the frequency range 0.03-0.125Hz in their study of the same data set as ours.Their data are therefore also affected by crustal reverberations.For data filtered in a large frequency range, the effect of the crustal reverberations will depend on the dominant frequency in the data and may change from event to event.Obayashi et al. (2004) used the waveform of the direct P wave to correct the PP waves from the crustal corrections at the reflection point, ensuring that the correct spectrum is used for each event.Ideally, this should also be done in regional body-wave tomography, but, as opposed to the PP case, we do not have a reference wave that has not propagated through the crust.An average waveform may be used as reference, and the delays may be measured by convolving the reference waveform with the impulse responses at the different stations and cross-correlating with a reference station.For the more narrow-filtered data used in finite-frequency tomography, Yang and Shen (2006) have shown that variations in the crustal corrections due to varying spectral content within each band can be neglected, and finite-frequency crustal corrections might therefore be done with a simpler procedure than in the broadband case.
We conclude that in classical continental environments, the frequency-dependence of the crustal corrections is not associated with variations in crustal thickness, but with variations in low-velocity layers in the upper crust.This dependence may have the same order of magnitude as the crustal corrections themselves and cannot normally be neglected.Ray theory produces valid corrections in continental regions free from low-velocity sedimentary layers and non-thinned crust.

Figure 1 .
Figure 1.Topographic map of the study area with location of seismological stations and national borders.Numbers xx are short for stations NWGxx in the MAGNUS network.

Figure 2 .
Figure 2. Maps of earthquake sources for P-and S-wave travel times.The events displayed have given travel time residuals in at least one frequency band.Left: events for P waves.Right: events for S waves.

Figure 3 .
Figure 3. High-frequency P-wave traces before and after alignment by ICCS.(a) Seismograms after alignment by theoretical arrival times from AK135 and data rejection in step 2 (Sect.2.2).(b) Seismograms after alignment by ICCS and additional rejection of incoherent traces in step 3 (Sect.2.3).The colours are used to distinguish the different traces better and no significance is attached to them.

Figure 4 .
Figure 4. P-wave stack and individual traces after alignment by ICCS.(a) Low-frequency band.(b) High-frequency band.The array stack is coloured blue, the accepted traces are red and green traces are the ones that will be rejected from analysis before step 4 (Sect.2.4).

Figure 5 .
Figure 5. High-frequency S-wave stack and traces in the first iteration of the ICCS algorithm.The array stack is coloured blue and the individual traces red.(a) High-frequency traces aligned according to the theoretical arrival time predicted with AK135.(b) High-frequency traces aligned according to the low-frequency travel time residuals from the MCCC analysis in step 4 (Sect.2.4).

Figure 6 .
Figure 6.Thickness of the sedimentary layer and of the crust in the study region.

Figure 7 .
Figure 7. S-wave velocity models down to 50 km depth at a selection of stations located along two east-west profiles.

Figure 8 .
Figure 8. Impulse response of the crust for P and SH waves at stations NWG12 and DK08, without (blue) and with bandpass-filtering (red for high-frequency band and green for low-frequency band).The incidence angle corresponds to 75 • epicentral distance.The time axis shows propagation time from 50 km below sea level.

Figure 9 .
Figure9.P (upper panels) and SH (lower panels) crustal corrections for all stations in the network calculated with ray theory (left window), reflectivity method in the high-frequency band (middle window) and reflectivity method in the low-frequency band (right window).All corrections are with respect to the average correction with ray theory.

Figure 11 .
Figure 11.Same as Fig. 10 but for the SH waves.