2.1 Seismic network

We use continuous raw seismic waveform data from the AASN (Hetényi et al., 2018a) that operated from January 2016 until April 2019 and together with permanent networks covered the broader European Alpine region (inverted triangles in Fig. 1). The main goal for the deployment of AlpArray was to provide an updated image of the crust and mantle structure beneath the European Alps orogen (e.g., Hetényi et al., 2018a).

We complement the AASN and permanent networks with data from temporary seismic deployments such as the Eastern Alpine Seismic Investigation (EASI; Hetényi et al., 2018b), the China–Italy–France Alps seismic transect (CIFALPS; Zhao et al., 2015, 2016) and the Pannonian–Carpathian–Alpine Seismic Experiment (PACASE; Hetényi et al., 2019) to expand the coverage and density of the seismometers. The EASI seismic network deployed 55 broadband seismometers along a north–south transect from the Bohemian Massif to the Adriatic coast (triangles in Fig. 1) from late 2014 until late 2015. The CIFALPS also deployed 55 broadband seismometers along a WSW-ENE transect across the southwestern Alps (squares in Fig. 1) for 13 months (2012–2013). The PACASE seismic network continued operating over 100 temporary broadband seismometers from AASN until 2022 and deployed over 50 additional temporary sites in Slovakia, Hungary, and southern Poland, to the east of the AASN (diamonds in Fig. 1).

Figure 1Distribution of three-component broadband seismometers of the AlpArray Seismic Network (AASN), the China–Italy–France Alps seismic transect (CIFALPS), the Eastern Alpine Seismic Investigation (EASI), and the Pannonian–Carpathian–Alpine Seismic Experiment (PACASE) used in this study (Hetényi et al., 2018a, b; Zhao et al., 2016; Hetényi et al., 2019). Seismic stations from AASN including the permanent stations, CIFALPS, EASI, and PACASE are respectively shown as inverted triangles, squares, triangles, and diamonds and are colored according to the number of waveforms used (Z-N-E component triplets).

2.2 Teleseismic events

Figure 2 shows the epicenters of the 1687 teleseismic earthquakes that we use for the receiver function calculations. We selected teleseismic events with M>5.5 and epicentral distance, Δ, between 30 and 90^∘. These epicentral distances and earthquake magnitude ranges were used to avoid the upper mantle triplications and regional and core phases. We exclude earthquakes with large magnitudes M>8.5 because their waveforms have long and complex source-time functions that can also be contaminated with signals from large aftershocks. In our case, this threshold only discarded one earthquake (M=8.6 earthquake that took place off the west coast of northern Sumatra in April 2012 and was recorded by the CIFALPS seismic network).

Figure 2Distribution of teleseismic events used in this study shown as circles colored according to hypocentral depths (USGS catalog; https://earthquake.usgs.gov/earthquakes, last access: October 2022). The study area is shown by the blue box. Dashed black lines show the 30 and 90^∘ epicentral distance from the center of the study area.

3.1 Receiver function calculations

We first cut continuous seismic waveform data around the theoretical P-wave arrival times (120 s before and after) of the teleseismic events (1687 events in total) from all available seismic stations. We downsample the cut waveforms to 20 samples per second. We discard any incomplete trace or any incomplete Z-N-E triplet. This provides us with 521 762 Z-N-E waveform triplets in total.

We then rotate the horizontal components into true N and E components based on the alignment stated in the station metadata (i.e., azimuth and dip values of each channel stored in the inventory file). This step is particularly crucial for the randomly oriented horizontal components of free-fall ocean bottom seismometers of the AlpArray Seismic Network in the Ligurian Sea.

Following Hetényi et al. (2018b), we apply a first quality-control (QC) step that will help us remove anomalous traces. Particularly, this step includes the calculation of two signal-to-noise parameters: (1) the ratio between peak amplitude to background noise amplitude, and (2) the ratio between peak amplitude and background noise root mean square (rms). We also calculate the rms value of all the available seismic traces for each teleseismic event and use the median rms value for that event, median(rms). Then we discard any trace with an rms that does not meet the following criteria for that event:

where c₁ and c₂ are sensitivity parameters and equal to 10 and 0.1, respectively. The c₁ and c₂ values are based on empirically defined values (for more details, see Hetényi, 2007). This quality control ensures that we discard noisy waveforms or waveforms that contain biases introduced by local noise or problems with the seismometers. We apply this quality-control step to all three (Z-N-E) component seismograms separately.

The horizontal seismogram components are then rotated to the radial (R) and transverse (T) components, using the theoretical back-azimuth value, ϕ. This step helps to isolate the energy from the direct P wave and the converted S waves.

At this stage, we apply a second quality-control step to ensure that we keep traces with strong signals (i.e., large amplitudes). To do this, we use an algorithm based on observing changes in the amplitudes of the seismic traces within two moving time windows (i.e., ratio of short-term average, STA, to long-term average, LTA, algorithm; Allen, 1978) that is most commonly used in earthquake detection applications (e.g., Michailos et al., 2019, 2021). We apply a low-pass filter of 1 Hz to the radial component waveforms and use STA and LTA values of 3 and 50 s, respectively. We keep the waveform traces that have STA $/$ LTA ratio values larger than 2.5 and discard traces with weak signals.

Prior to the deconvolution, we (1) trim the 240 s long radial and vertical component traces 40 s before and 60 s after the P arrival, (2) remove the mean, and (3) taper the traces using a 15 s Hann window at both ends to avoid filter artifacts. We then apply a second-order 0.05–1 Hz Butterworth bandpass filter that eliminates the noisiest frequencies.

We use the iterative time-domain deconvolution approach of Ligorría and Ammon (1999) with 200 iterations to remove the effect of source-time functions. In this approach, spikes are added iteratively to the trial RF that is convolved with a Gaussian with a half width of 1 Hz. The convolution of the series of Gaussian spikes with the vertical component approximates the radial component. The predicted final radial components have high cross-correlation values (RF fit values; see Figs. A1 and A2) when compared to the observed radial traces. We chose the approach of Ligorría and Ammon (1999) among other available approaches (e.g., Gurrola et al., 1994, 1995; Park and Levin, 2000; Helffrich, 2006) based on its stability (see discussion in Lombardi et al., 2008). It should also be noted, however, that for high-quality datasets, most deconvolution methods work well, and the advantages of one method over another are insignificant (e.g., Ligorría and Ammon, 1999).

We then apply the last quality control to the calculated RFs, similar to Hetényi et al. (2015), Singer (2017), Hetényi et al. (2018b), Scarponi et al. (2021), that is based on the signal-to-noise ratio and the amplitude of the direct wave signal. For the first part of this quality control, we calculate (1) the rms of the noise from 30 to 10 s before the P-wave arrival (rmsnoi), and (2) the rms of the signal between 2 and 30 s after the P arrival (rmssig). We keep the RFs that meet the following criteria:

In addition, we require that (1) the timing of the maximum amplitude peak should be between 0 and 2 s and that (2) its amplitude should be positive and have a value between 0.05 to 0.8. The first criterion ensures that the dominant arrival is either the direct wave, expected to appear at t=0 s in Z-R RF, or the conversion at the sediment–basement interface. The second requirement on the amplitude value also helps to discard seismic traces with incorrect gain values or deconvolution artifacts (e.g., Hetényi, 2007). We discard RF traces with overly large rms values (threshold equal to 0.07; Fig. A3) which represent the top ∼2 % the traces. The results are not sensitive to the choice of this threshold; however, the chosen value enables us to discard the noisiest RF traces that are expected to mostly affect and contaminate our final migrated images.

As the last check, we visually inspected the RF traces (radial components) plotted versus their back azimuths from each individual seismic station. This final step ensures that we identify and discard any traces with low-quality results.

Figure 3Workflow summarizing the processing steps followed for (a) the receiver function calculations and (b) for the time-to-depth migration. These steps are described in detail in Sect. 3.

Download

3.2 Common conversion point stacking

To examine the spatial variability of crustal structures, we use the common conversion point (CCP) stacking techniques that take advantage of dense arrays of seismic stations to provide images of discontinuities in the crust and upper mantle (e.g., Yuan et al., 1997; Kosarev et al., 1999; Ryberg and Weber, 2000; Zhu, 2000).

For each seismic trace, we backtrace the ray paths for S waves, assuming the theoretical P-arrival slowness. This is equivalent to assuming that all energy in the seismic trace has resulted from direct P–S conversions at horizontal boundaries. Because the study region is characterized by very significant topography and includes stations below sea level and on mountain tops, we take into account the elevation of each seismic station. Given that our study area spans approximately 2000 km, we also take into account Earth's sphericity by implementing a new CCP stacking code in a spherical coordinate system. This way the obtained profiles and maps should represent structures closer to reality.

For the time-to-depth migration, we use EPcrust (Molinari and Morelli, 2011), a 3D P- and S-wave reference velocity model based on a selection of previous geological and geophysical studies (e.g., surface-waves velocity models, active seismic experiments, receiver functions, noise correlations, and geological assumptions), that covers the entire study region. We opt to use EPcrust instead of a global 1D model (i.e., iasp91 velocity model; Kennett and Engdahl, 1991) in order to include the local structure variations within the study area. In regions with thick sedimentary layers, such as the Pannonian Basin or Po Plain, the absence of a sedimentary layer in the iasp91 model provides phase arrival times with significant time shifts. At the same time, both higher-resolution Vp and Vs models are only locally or sub-regionally available in the Alpine region. but none of these models cover the entire study area. Therefore, we have chosen to use EPcrust instead of compiling an ad hoc composite velocity model from various locally available models to ensure internally consistent results.

The grid spacing of EPcrust is 0.5^∘ in latitude and longitude and consists of three layers in depth (i.e., sedimentary, upper crust, and lower crust). We use a linear interpolation method that allows us to estimate P- and S-wave velocities at any given point within the velocity model. For the depths below the EPcrust Moho, we use the velocity values of the lower crust in order to avoid mismigration with mantle velocities when the observed Moho conversion phase is later than expected from EPcrust. Conversions from below the Moho will be mismigrated, which in the context of this study focusing on Moho depth, is not of concern. We then define a 3D spherical migration grid spaced 0.05^∘ in latitude and longitude and at 0.5 km in depth. Within this grid, we calculate the theoretical ray paths with a smaller, 0.25 km increment, taking into account the station elevation information. The P- and S-wave velocities along the ray paths are sampled from the interpolated EPcrust velocity model. We back-project the converted wave amplitudes of the receiver functions along these theoretical ray paths and stack them within the defined 3D spherical migration grid. This 3D grid of stacked RF amplitudes, therefore, highlights any increase or decrease in seismic velocity with depth, which primarily represents structures associated with sharp velocity changes such as the Moho interface.

4.1 Receiver functions

We obtain 112 205 quality-controlled RFs in total. Figure 4 depicts the number of RFs at each seismic station. As a general rule, we expect that there will be more RFs at stations that have longer operational times or have high-quality waveforms (e.g., stations located on bedrock). Consequently, seismic stations in southern France, the Po Basin, the Pannonian Basin, the OBS deployment in the Ligurian Sea, and the CIFALPS and EASI deployments have relatively low numbers of RFs compared to other stations because of shorter deployment times, the presence of sedimentary strata, which cause higher noise levels on the horizontal components, and because of high noise levels on the horizontal components of OBS. Furthermore, reverberations within the sedimentary layer or the water layer for OBS can obscure the Moho phase and make the RF difficult to interpret.

Figure 4Map showing the number of receiver functions calculated for each seismic station. Seismic stations from AASN including the permanent stations, CIFALPS, EASI, and PACASE are respectively shown as inverted triangles, squares, triangles, and diamonds and are colored according to the number of receiver functions available (see color scale). The seismic stations with blue edges marked with a, b, c, and d indicate the four stations from which we plot the receiver function stacks in Fig. 5.

We plot the calculated RFs as a function of their back-azimuth values for each station (see examples in Fig. 5) to observe any differential times, amplitudes, and polarities of the phases that may contain additional information on the inclination of the discontinuities (e.g., Cassidy, 1992). This kind of graph also helps in providing an indication of the quality of the dataset (i.e., how easily identifiable is the Ps-Moho phase, what is the total back-azimuthal coverage, how does the variation of the sharpness of the Ps-Moho phase for different back azimuths look, and is there a dipping Moho interface). Systematic variations of the timing and amplitude can indicate the presence of dipping boundaries, 360^∘ periodicity, or the effect of anisotropy, and 180^∘ periodicity for horizontal anisotropy (e.g., Cassidy, 1992; Levin and Park, 1997; Savage et al., 2007).

In Fig. 5, we show RFs at four different seismic stations that are located in different geologic environments in the broader European Alpine region (i.e., foreland/Molasse basin, the high Alpine region, and the Pannonian Basin) and that have a varying number of RFs due to their respective operational times during the AASN deployment. The back-azimuthal coverage is generally continuous with some gaps for the southerly and westerly directions. The amplitude of the direct P-wave is variable, presumably due to changes in the angle of incidence (i.e., epicentral distances) correlated with back azimuth. We observe sharp P and Ps-Moho phases for all seismic stations at around 0–2 and 2.5–4 s, respectively. For seismic stations located in regions with thick sedimentary layers (e.g., Pannonian Basin; Fig. 5d), these phase arrivals appear with a ca. 1–1.5 s of delay in time, similar to what Kalmár et al. (2019) observed. The delay of the initial P phase is caused by interference with the conversion at the sediment–basement interface, while the delay of the Moho phase should simply be a travel time effect, which can be corrected by a velocity model including sedimentary layer information. Seismic stations with relatively fewer RF traces (Fig. 5a) nevertheless have clear high-quality signals of comparable quality to those with longer operational times (Fig. 5b and c).

Figure 5Receiver function stacks from selected seismic stations plotted versus back azimuth. We apply a normal moveout correction and stack the traces in 20^∘ wide back-azimuth bins. The number of traces contributing to each stack is shown on the right side of each panel. Dashed vertical lines indicate time delays at 0, 5, and 10 s from the P wave arrival time. Panels (a) and (b) show the RFs from two seismic stations Z3.A196A and Z3.A115A in the European foreland/Molasse Basin region. Panel (c) depicts the RFs from station CH.BERNI in the high Alpine region, and panel (d) shows station Z3.A267A in the Pannonian Basin.

Download

4.2 Receiver function migration

To assess the coverage, we compare the station spacing with the horizontal offset of the piercing points. The lateral offset between the station location and the piercing points at the Moho is roughly half the Moho depth. Hence for a Moho ranging from ca. 22–25 km depth (beneath the Pannonian Basin) to ca. 55–60 km (beneath the Alps), the corresponding range of horizontal offsets of the piercing points is about 11–30 km. Since the nominal station spacing is 50 km, over most of the network the piercing points of neighboring stations are not significantly overlapping. We calculate the piercing points, using the EPcrust velocity model, at 35 km depth as a compromise between areas of expected shallower and deeper Moho depths (Fig. 6). The piercing points are well distributed and do not leave major gaps in coverage – with the exception of some seismic stations in the Po Basin, in the Pannonian Basin, and the OBS stations in the Ligurian Sea where the coverage is less dense due to either short operational times or the strong-quality criteria applied.

Figure 6Map depicting the piercing points (blue crosses) at 35 km depth computed for each seismic station (inverted triangles) using the EPcrust velocity model (Molinari and Morelli, 2011). The size of the crosses is based upon a Fresnel volume estimate of ∼9 km for a typical frequency of 1 Hz.

Having calculated a 3D grid of stacked migrated RF amplitudes, we are able to extract 2D migrated receiver-function cross-sections at any location within our study region. Figure 7 shows the crustal structure along three cross-sections. Namely, cross-section (A–A^′) is located beneath the western Alps collocated with the CIFALPS deployment to benefit from high station density, cross-section (B–B^′) beneath the central Alps, and cross-section (C–C^′) beneath the Pannonian Basin and the eastern Alps. In cross-section A–A^′, we observe strong migrated amplitude signals from the European Moho, which dips gently to the east–northeast. The European Moho (EU in Fig. 7) starts from depths of approximately 30 km in the west reaching gradually depths of ca. 50–55 km. The signal then weakens at around 3–3.5^∘ distance along cross-section A–A^′, and for the rest of the cross-section, the Adriatic Moho (AD in Fig. 7) is generally harder to interpret. The crustal geometry in cross-section A–A^′ is in good agreement with the result of Zhao et al. (2015) who used the CIFALPS seismic data and Monna et al. (2022) who jointly inverted P and S receiver functions using AASN seismic data. In cross-section B–B^′, we also observe strong amplitudes for the European Moho (∼30 km in the northwest that gently dips down to ∼55 km beneath the Alps). The Adriatic Moho is a bit easier to identify, with respect to cross-section A–A^′, with depths of around 30 km (around 4^∘ in B–B^′). Farther southeast, the B–B^′ profile continues through the Po Basin (PB in Fig. 7) where the signal is harder to interpret, as the Moho no longer appears as a continuous interface; possibly this is caused by interference from multiples of the sediment–basement interface, but also there are fewer high-quality RFs available for this area. Finally cross-section C–C^′ starts from the Po Basin where again no clear continuous interface is visible; the whole section is harder to interpret than previous ones, especially beneath the eastern Alps. The Adriatic Moho at offsets 1–3^∘ is between 35 and 40 km deep. Farther northeast, the European Moho lies between 25 and 30 km depth beneath the Vienna Basin and past the western Carpathians (VB and WC in Fig. 7).

Figure 7Migrated receiver-function cross-sections in spherical coordinates along lines A–A^′ (top panel), B–B^′ (middle panel), and C–C^′ (lower panel) as marked in the inset map. Note that cross-section D–D^′ is shown in Fig. 9. White and gray diamonds represent the manually determined certain and uncertain Moho picks, respectively. For comparison, the solid black and dashed black lines depict the Moho depths calculated by Grad and Tiira (2009) and Spada et al. (2013), respectively. EU: European Moho, AD: Adriatic Moho, PB: Po Basin, VB: Vienna Basin, and WC: West Carpathians. All cross-sections are smoothed and filtered and have no vertical exaggeration, profile curvature represents real Earth sphericity. The inset map shows the location of seismic stations used with inverted blue triangles. Solid lines indicate the cross-section locations with tracks every 1^∘ distance.

4.3 Moho depth map

To compile the new Moho depth map, we manually pick Moho depths on 21 equally spaced cross-sections along the meridians (3–23^∘) and another 10 equally intercepting cross-sections along the parallels (43–50^∘). Refer to the Supplement for all cross-sections and manually defined picks. We determine the Moho depth picks based on the strength and continuity of the amplitude signal. We assign the picks on the maximum of the signal. In cases where the signal appears to be unclear (weaker signal) or not completely continuous, we assign uncertain Moho depth picks. For examples where the signal is double (overlapping Moho discontinuities), we choose to pick the shallower signal. Double signals might be associated with underthrusted lower crustal slivers or subduction (e.g., Spada et al., 2013; Hetényi et al., 2018b; Mroczek et al., 2023); therefore picking the shallower Moho ensures the bottom of a continuous crustal stack is picked.

This analysis provides us with a dense grid of Moho depths in the broader European Alpine region (Fig. 8a). We interpolate the Moho depth points using a linear interpolation method to obtain a continuous image of the Moho depth distribution (Fig. 8b). We also test interpolating (1) the Moho depth picks separately for each plate (i.e., European, Adriatic, and Ligurian) and (2) using only the certain Moho depth picks (excluding uncertain) and obtained very similar patterns (Figs. A4 and A5). The plate boundaries are taken from Spada et al. (2013)'s Moho map, although in some locations our migrated RF profiles point to slightly different plate limit positions. Moho depths vary from ∼30 to ∼60 km in the broader European Alpine region, with the deepest values (>50 km) observed near the highest topography of the Alps close to the plate boundary between the European and Adriatic plates. The majority of the European plate has Moho depths of around 30 km, which is typical for the continental crust in western Europe. We are able to estimate very few Moho depth picks in the Ligurian Sea and no picks in the Adriatic Sea due to the lack of deployed stations.

Figure 8(a) Map showing the spatial variations of Moho depths for the broader European Alpine region. Circles and diamonds with black edges depict certain and uncertain manually picked Moho depths, respectively. Numbers indicate the names of the cross-sections that are used for making the Moho picks (see manually determined Moho picks for each cross-section in the Supplement). Thick black lines depict the plate boundaries between the European, Adriatic, and Ligurian plates adapted from Spada et al. (2013). (b) Contoured distribution of all Moho depths. Empty boxes depict regions where the interpolation method did not provide any results.

5.1 Comparison with previous Moho depth models

Our Moho depth estimates and migrated images are comparable with and generally supported by the results from previous studies in the region (e.g., Grad and Tiira, 2009; Spada et al., 2013), with some differences which we attribute to a more complete spatial data coverage and the different approaches used. In general, our Moho estimates are relatively shallower in the European plate and relatively deeper in the Adria plate when compared to the two previous studies (Fig. A6; Grad and Tiira, 2009; Spada et al., 2013). The main differences are located along the tectonic plate boundaries, especially along the Ivrea-Verbano zone where estimates around the protruding Ivrea geophysical body cause large differences. This is mostly due to differences in the coverage of seismic stations and the respective method of interpolation. Denser arrays can provide a higher local resolution for the Moho (e.g., 5 km spacing of IvreaArray; Scarponi et al., 2021). When all quantified Moho depth differences are considered (Fig. A7), our Moho depth estimates are on average 0.6 km deeper than those of Spada et al. (2013) and 0.4 km shallower than those of Grad and Tiira (2009). The standard deviation of the difference distribution – which includes the high values mentioned above – is respectively 7.6 and 6.4 km.

In profile A–A^′ of Fig. 7, our deepening Moho signal and picks of the European plate can be followed to 3.25^∘ distance at least, while that of Grad and Tiira (2009) starts transitioning to a shallower depth at 2.75^∘ distance (∼50 km to the west), and that of Spada et al. (2013) features a sharper jump at around 3.1^∘ distance. On the Adria side, the Moho depth estimates of the three models differ, which is expected given the complexity related to the shallow intruding Ivrea geophysical body that requires more detailed analyses (e.g., Scarponi et al., 2021). In profile B–B^′, the situation at the Europe–Adria limit is somewhat similar to that in profile A–A^′. Our European Moho signal and picks can be tracked until 3.7^∘ distance, while Grad and Tiira (2009)'s estimates become shallower already at 3^∘ distance (70 km to the NW), and Spada et al. (2013)'s jump is at 3.5^∘ distance. The continuity of the European Moho signal beneath Adria in both profiles in our images could be explained by the teleseismic waves imaging this zone from below, while a good portion of Spada et al. (2013)'s dataset is from active seismic imaging from above, hence determining the shallower discontinuity. In that sense, defining the boundary between two plates in a collision zone goes beyond the task of tracing a line on a map, and it requires 3D analysis and the representation of overlapping Mohos, something that was done by Spada et al. (2013) in some regions. From a geodynamic perspective, this is not surprising, as during collisional processes, different layers of the lithosphere decouple from each other, may deform independently from each other, and also imbricate between plates. In profile C–C^′, the three Moho models show fewer differences, estimated at ∼5 km and locally at 10 km in depth. The largest difference at around 4^∘ distance reveals a shallower Moho in our model compared to Grad and Tiira (2009), which could reflect a thinner crust in the transition zone of the easternmost Alps to the Pannonian Basin.

5.2 Comparison of RF migrated profiles using different velocity models

To examine the influence of the laterally heterogeneous velocity model used for migration, we compare migrated profiles using two different velocity models (i.e., iasp91 and EPcrust in our case). For this comparison, we choose a cross-section along the EASI profile. By doing this, we are also able to compare our results to the previous study of Hetényi et al. (2018b). Figure 9 depicts the migrated cross-section images from north to south along the EASI seismic network using the EPcrust velocity model (top panel) and the global 1D iasp91 velocity model (lower panel). Although the overall pattern is the same, there are significant differences of more than 5 km in the absolute inferred Moho depth that are mostly due to the absence of a sedimentary layer in the iasp91 model (depths of <5 km). In addition, these images are also in good agreement with the results of Grad and Tiira (2009), Spada et al. (2013), and Hetényi et al. (2018b) to the first order. A notable zone of difference is the root of the Alps, between 3.5 and 4^∘ distance (Fig. 9), where Hetényi et al. (2018b) interpreted the Moho as the bottom of a broad velocity gradient, while the two other models relied mostly on active seismic data and interpreted the top of that zone. Further discussion would require more detailed processing, such as RF multiples, initiated in Hetényi et al. (2018b) and continued in Mroczek and Tilmann (2021).

Figure 9Comparison between migrated receiver-function cross-sections along the EASI seismic network (cross-section D–D^′ in Fig. 7) calculated using the EPcrust velocity model (top panel) and iasp91 velocity model (lower panel). The solid, dashed, and dotted lines show the Moho depths estimated by Grad and Tiira (2009) that is a compilation of Moho depths, Spada et al. (2013) who used a modified iasp91 velocity model for the migration calculations, and Hetényi et al. (2018b) who utilized EPcrust for migration, respectively. Stripes between the panels highlight the location of the orogenic root, where Hetényi et al. (2018b) interpreted the Moho at the bottom of a broad velocity-gradient imaged by teleseismic waves from below (represented by dots), while other studies (e.g., Spada et al., 2013) interpreted the Moho using controlled-source waves reflecting from the top of this gradient.

Download

5.3 Limitations

The semi-automatic workflow used here includes strict quality criteria that discard almost $\mathrm{3}/\mathrm{5}$ of the available seismic waveforms. This order of magnitude of discarded waveforms is comparable to other receiver-function studies (e.g., Hetényi et al., 2018b; Scarponi et al., 2021; Kalmár et al., 2021). However, by using these strict quality criteria, (1) we ensure that our final results are less likely to be contaminated by low-quality RFs traces, and (2) we are able to automatically process a large number of seismic waveforms and simplify manual inspections that would be extremely time-consuming. We note that these quality-control criteria and the frequency band used are tuned for identifying crustal structures and that other applications of RFs may require adjusted selections.

Despite the high density of seismic waveform data available, we emphasize that our Moho depth estimates have limitations and uncertainties that are inherently difficult to assess. For example, the Moho depth estimates can contain errors associated (1) with the quality and consistency of the manual Moho picks and (2) with the inherited uncertainties from the EPcrust velocity model used for migration (e.g., see differences in Moho depths in Fig. 9). Furthermore, we have solely focused on direct Ps conversions and have not investigated either multiples or illumination directions for effects of dip or anisotropy (e.g., Bianchi and Bokelmann, 2014; Bianchi et al., 2021). Ultimately, a major source of depth uncertainty is the 3D variability of the real seismic wave speed structure with respect to the three-layer, 0.5^∘ distance resolution EPcrust model. This and lower success in the RF quality control makes picks and interpretations particularly vulnerable in areas of sedimentary basins, and we would recommend closer (re)analyses for studies focusing on those areas. Extending the analysis by including the crustal multiples (e.g., PpSs, PpPs) is a potential future direction that would allow an estimate of Moho depth without assuming a fixed $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio over local to regional scales. For fully 3D Vs structure inversion based on RFs, we refer to Colavitti and Hetényi (2022).