Crustal Thinning From Orogen to Back-Arc Basin: The Structure of the Pannonian Basin Region Revealed by P -to- S Converted Seismic Waves

We present the results of P -to- S receiver function analysis to improve the 3D image of the sedimentary layer, the upper crust, and lower crust in the Pannonian Basin area. The Pannonian Basin hosts deep sedimentary depocentres superimposed on a complex basement structure and it is surrounded by mountain belts. We processed waveforms from 221 three-component broadband seismological stations. As a result of the dense station coverage, we were able to achieve so far unprecedented spatial resolution in determining the velocity structure of the crust. We applied a three-fold quality control process; the first two being applied to the observed waveforms and the third to the calculated radial receiver functions. This work is the first comprehensive receiver function study of the entire region. To prepare the inversions, we performed station-wise H -Vp/Vs grid search, as well as Common Conversion Point migration. Our main focus was then the S -wave velocity structure of the area, which we determined by the Neighborhood Algorithm inversion method at each station, where data were sub-divided into back-azimuthal bundles based on similar Ps delay times. The 1D, nonlinear inversions provided the depth of the discontinuities, shear-wave velocities and Vp/Vs ratios of each layer per bundle, and we calculated uncertainty values for each of these parameters. We then developed a 3D interpolation method based on natural neighbor interpolation to obtain the 3D crustal structure from the local inversion results. We present the sedimentary thickness map, the first Conrad depth map and an improved, detailed Moho map, as well as the first upper and lower crustal thickness maps obtained from receiver function analysis. The velocity jump across the Conrad discontinuity is estimated at less than 0.2 km/s over most of the investigated area. We also compare the new Moho map from our approach to simple grid search results and prior knowledge from other techniques. Our Moho depth map presents local variations in the investigated area: the crust-mantle boundary is at 20–26 km beneath the sedimentary basins, while it is situated deeper below the Apuseni Mountains, Transdanubian and North Hungarian Ranges (28–33 km), and it is the deepest beneath the Eastern Alps and the Southern Carpathians (40–45 km). These values reflect well the Neogene evolution of the region, such as


Geological and Geodynamic Setting
Our study area, the Pannonian Basin of Central Europe (Figure 1) is a continental back-arc basin system surrounded by the Alpine, Carpathian, and Dinaridic orogens. Its Miocene extensional basin formation and crustal thinning followed a pre-Neogene orogenic evolution that resulted from the opening and subsequent closure of the Neotethys and Alpine Tethys oceanic realms (Csontos & Vörös, 2004). By the end of Paleogene times, collision and successive events of shortening created a thickened orogenic area composed of two tectonic megaunits with different paleogeographic positions (Balla, 1986). The Alps-Carpathians-Pannonia (AlCaPa) in the NW is an Adriatic-derived orogenic unit that was extruded from the Eastern Alps, while the Tisza-Dacia in the SE derives from the Dinaridic continental collision zone (Figure 1; Schmid et al., 2020). The Neogene formation of the Pannonian Basin, coupled with extrusion of the pre-Neogene basement units from the adjacent orogens (Horváth et al., 2006), has created large amounts of translations and opposite sense rotations, that is, counterclockwise in AlCaPa and clockwise in Tisza-Dacia accompanying their extension (Márton & Fodor, 2003). These units became juxtaposed in the Middle Miocene along a wide transcurrent fault system, that is, the Mid-Hungarian (Fault) Zone (Csontos & Nagymarosy, 1998).
Neogene back-arc extension of the Pannonian Basin was coeval with the contraction recorded at the exterior of the Carpathians orogen (Horváth et al., 2006;Matenco & Radivojević, 2012). Available estimates on KALMÁR ET AL.  the amount of Miocene extension are in the order of 150-250 km (Balázs et al., 2018;Horváth et al., 2015;Ustaszewski et al., 2014) that resulted in high surface heat flow (Lenkey et al., 2002) and thin crust and lithosphere (Horváth et al., 2015;Posgay et al., 1991;Tari et al., 1999) and connected to successive volcanic events (Harangi & Lenkey, 2007;Kovács et al., 2007). From ca. 20 Ma, syn-rift extension created metamorphic core complexes and similar extensional domes at the W and SE basin margins exhuming deep crustal rocks (Fodor et al., 2008;Tari et al., 1999) and the gradual opening of deep depocentres filled by sediments reaching more than 7 km in thickness in the most extended areas, such as in the Danube Basin in the NW and in the SE part of the Great Hungarian Plain (Figure 1; Balázs et al., 2018).
The syn-rift extension and half-graben formation in the Pannonian Basin ceased between 13 and 9 Ma (e.g., Balázs et al., 2016;Matenco & Radivojević, 2012). The subsequent basin evolution was controlled by postrift thermal cooling and the basin-wide inversion during the northward drift and CCW rotation of the Adriatic indentation (Bada et al., 2007;Ustaszewski et al., 2014). Such inverted structures developed mainly in the SW Pannonian Basin and along the Mid-Hungarian (Fault) Zone (Figure 1). This deformation is still presently active creating intraplate seismicity connected to inverted normal faults and fault propagation folds and sinistral and dextral strike-slip faults in the center and along the Dinaridic margin of the basin, respectively (Bada et al., 2007;Tóth et al., 2002).

Previous Crustal Studies
In the previous century, few studies focused on the determination of the depth of crust-mantle boundary in the Pannonian Basin (e.g., Gálfi & Stegena, 1957, 1960Mituch, 1964). In the 21st century, further studies dealt with Moho discontinuity (e.g., Grad et al., 2009;Horváth et al., 2015), but these suffered from large data gaps due to the scarcity of instrumental coverage and had to rely on background models of crustal thickness. Geochemical studies based on middle and lower crustal xenoliths from the western Pannonian Basin inferred an initially overthickened and dry crust prior to Miocene extension (e.g., Németh et al., 2015;Török et al., 2013). Local and unified crustal structure from receiver function studies were so far not implemented in the Pannonian Basin.
The first digital active seismological studies of the crustal structure in the Pannonian Basin used 2D seismic reflection and refraction profiles such as the CELEBRATION 2000 (Guterch et al., 2003) and ALP 2002 (Bielik et al., 2018;Brückl et al., 2007). The first passive seismological study was conducted by travel time tomography calculations (Bus, 2001;Wéber, 2002). The first S-wave velocity and crustal thickness investigation based on receiver function inversion was published a few years later (Bus, 2003;Hetényi & Bus, 2007). Several studies presented global Moho depth and crustal velocity structure throughout Europe, including the Pannonian Basin (e.g., Grad et al., 2009;Tesauro et al., 2008), but generally at lower spatial resolution.
To overcome the difficulty of the poor station coverage, passive seismological campaigns were conducted by a consortium led by the University of Leeds, UK. The first campaign, that is, the Carpathian Basin Project (CBP) took place between 2005 and 2007 in the western part of the Pannonian Basin and the second campaign, that is, South Carpathian Project (SCP) took place between 2009 and 2011 in the eastern part of the Pannonian Basin and the Southern Carpathians. The CBP project resulted in a teleseismic tomography study (Dando et al., 2011) and detailed receiver function analysis beneath 46 broadband stations (Hetényi et al., 2015). The SCP data set were exploited for body wave tomography (Ren et al., 2012), ambient noise tomography including CPB stations too (Ren et al., 2013) and Moho depth determination from receiver function analysis (Kalmár et al., 2019). These studies all agreed on a thinned crust beneath the Pannonian basin area, and slow velocities of the sedimentary layers, yet a full seismological coverage of the area from the Eastern Alps to the Southern Carpathians was still lacking.
The next seismological breakthrough was brought about by a number of recently installed permanent stations and the AlpArray project starting at the end of 2015 (Gráczer et al., 2018;Hetényi, Molinari, et al., 2018). As a result, a number of new studies have recently been published at higher resolution, also on smaller areas of the region, such as the Moho depth determination from receiver function analysis in western part of the Hungary (Kalmár et al., 2018), S-wave velocity model based on ambient-noise tomography in the Vienna Basin region (Schippkus et al., 2018) and between the Eastern Alps and the Pannonian Basin (Szanyi et al., 2021), P-wave velocity imaging from traveltime tomography in the Pannonian Basin (Timkó et al., 2019) and crustal thickness determination based on receiver functions under the Dinarides and neighboring areas (Stipčević et al., 2020).

Research Goals
The recent developments in permanent networks and temporary (CBP, SCP, AlpArray) deployment campaigns that took place in the Pannonian Basin and its surroundings made it possible to carry out a comprehensive receiver function study on over 17 years of data. Using so far unprecedented station density of 50 km station spacing allows us achieving a much finer spatial resolution in the determination of the velocity structure of the crust.
Our semi-automated, three-tier quality control procedure enables the processing of large quantities of waveforms and screen out those with low signal-to-noise ratio (SNR) or lacking signals altogether, thus producing a high-quality set of receiver functions.
To demonstrate the consistency of our results regarding the crustal structure, we apply the H-Vp/Vs grid search method (Zhu & Kanamori, 2000), Common Conversion Point (CCP) migration (Zhu, 2000), and then mainly the neighborhood inversion algorithm (Sambridge, 1999) on receiver functions. We also develop a new 3D visualization method to obtain the 3D crustal structure of the Pannonian Basin. This method is able to describe velocity conditions and discontinuity of the layer (e.g., low-velocity sediments, Conrad, Moho). The upper and lower crustal thicknesses are very important factors in the Basin evolution procedure, but so far these have not been defined beneath our investigated area. Furthermore, we can map the transition zone between Pannonian Basin and surrounding mountains, which can provide information on the thinning of the crust beneath these areas. These results will help the better understanding of the geodynamical settings and tectonic evolution of the region.
Our uniform waveform pre-processing and quality control procedures applied before the receiver function analysis are the most important contributors to ensure robust, reliable results. Our results reveal a number of hitherto unknown geological and geophysical information of the Pannonian Basin and surrounding areas, thus allowing for more accurate and crustal tectonic interpretation of the region. Figure 2 shows the teleseismic earthquakes between 28° and 95° epicentral distances and magnitudes larger than 5.5 that we considered for the receiver function analysis. The data set consists of 3,098 events from the USGS catalog (https://earthquake.usgs.gov/earthquakes) between January 1, 2002 and March 31, 2019. We downloaded the broadband three-component waveforms of these events recorded at the stations in  Academy of Sciences of Ukraine (UA), the National Institute for Earth Physics (RO; 10.7914/SN/RO), Seismological Survey of Serbia (SJ), the Department of Geophysics, University of Zagreb and the Croatian Seismological Survey (CR) and the Slovenian Environment Agency (SL; 10.7914/SN/SL). For the permanent stations, we used all available data since they entered into operation up until March 31, 2019, and for the GeoRisk network, up until August 31, 2017. Most of these networks are members of the Central and Eastern European Earthquake Research Network, CE 3 RN (Lenhardt et al., 2021).

Quality Control and Receiver Function Analysis
To perform receiver function analysis we downloaded altogether 454,089 broadband three-component waveforms from the KRSO (permanent and CBP stations), ORFEUS-EIDA (AlpArray temporary network) and IRIS (SCP stations) data centers using a 900-s time window, 300 s before and 600 s after the theoretical first P-arrival time predicted by ak135 (Kennett et al., 1995).
The original raw three-component waveforms (ZNE data) were band-pass filtered (Butterworth-type) between 0.1 and 1 Hz before applying the quality control procedures. In order to select the optimal frequency bands for our methods, we tested various frequency ranges (0.1-0.5, 0.1-1, and 0.1-2 Hz). These ranges contain all frequencies which may bring useful information from teleseismic earthquakes. We found that the 0.1-2 Hz range contains too much high-frequency noise to reliably determine receiver functions for most part of the study area. Nevertheless, the other two frequency ranges contain enough information for the crustal discontinuities to describe major changes in the crust. The electronic supplement shows examples for three stations in different geological setting to demonstrate the suitability of our choice of frequency band.
We band-pass filtered the raw waveforms between 0.1 and 1 Hz in the H-Vp/Vs grid search and CCP migration method. We selected a somewhat narrower frequency band (0.1-0.5 Hz) in the nonlinear Neighborhood Algorithm inversion for relatively smoother convergence in the waveform fitting while still retaining information from the crustal structures (Hetényi, Plomerová, et al., 2018;Wu et al., 2018).
In order to automatically select the most suitable waveforms with the highest SNR, we applied our threestep quality control process (Kalmár et al., 2019).
In the first quality control step, QC1 we ran the STA/LTA detector (Trnkoczy, 2012) with a detection threshold of 3.5, to ensure that we indeed have a P arrival with reasonable SNR in the selected time window. QC1 screened out some 47% of the waveforms, including those that did not have all three components available.
We performed the second step, QC2 on the remaining 240,828 waveforms by calculating various ratios of shorter time windows signals as proposed by Hetényi et al. (2015) and Hetényi, Plomerová, et al. (2018): (a) for each event, only traces with an root mean square (rms) comprised between 0.4 and 2.5 times the median of all traces for that event were kept; (b) the maximum of the peak (−5 to 20 s with respect to P arrival) was required to be 6 times larger than the maximum of the background (from −30 to −5 with respect to P arrival); (c) the same maximum of the peak was required to be 9 times more than the rms of the background.
This step screened out further 29% of the waveforms that passed QC1 leaving us with 171,255 waveforms from 221 stations.
Having finished the first two quality control steps, we rotated the accepted ZNE waveforms into the ZRT coordinate system according to the theoretical back azimuth. We calculated the receiver functions using the iterative time domain deconvolution (Ligorría & Ammon, 1999) with 150 iterations. Finally, we performed the third quality control step, QC3 on the radial receiver functions based on the value and time delay of its dominant amplitude. We used a 30-s time window at these radial receiver functions, 5 s before and 25 s after first P-arrival time. Our QC3 rejected receiver functions as poor quality with positive highest amplitude of more than 0.1 and less than 0.8 between −0.2 and 2 s delay time. On stations without a sedimentary layer (e.g., GCIS on Figure 3), the maximum amplitude of the radial receiver function is located around 0 s. However, in the presence of sediments (e.g., BEHE on Figure 3) it can extend, but it is definitely located within 2 s from the direct P arrival in the whole study area. Therefore the selection of proper amplitude and time range is important, because this can handle effect of the low-velocity sedimentary layer (Kalmár et al., 2019;Yeck et al., 2013). These criteria left us with 30,216 high quality radial receiver functions.
The automated quality control steps allowed us to process a large number of waveforms with minimal manual intervention. Despite the quality control procedures screening out more than 80% of the original data set, we left with a considerable amount of radial receiver functions with indisputable quality that would also ensure the quality of the inversion results. Figure 3 illustrates the quality of stacked receiver functions (without move-out correction) along A-A′ profile before and after applying QC3. The noisiest receiver functions are located in the sedimentary part of the Pannonian Basin and typically recorded on temporary stations (e.g., CBP2K and CBP4K stations on Figure 3). We present more details about QC in the electronic supplement.

H-Vp/Vs Grid Search Method
To obtain a first rough estimate of the crust-mantle boundary, we performed the H-Vp/Vs grid search method (Zhu & Kanamori, 2000) on the radial receiver functions that passed quality control. Because of the diverse geological characteristics, we set the initial Vp velocity for each station separately based on previous studies (Behm et al., 2007;Grad et al., 2006;Janik et al., 2009;Tasarova et al., 2009). We set the limits of the Moho depth search range, H between 20 and 40 km in the basin area and between 20 and 50 km in the mountains. We set the Vp/Vs ratio search range, between 1.5 and 2. The selected H and Vp/Vs intervals are physically reasonable values for the region, and the optimal values obtained from the H-Vp/Vs grid search remained well within these limits in every single station.
Furthermore, we set the weights of Ps, PpPs, and PpSs + PsPs multiples separately for each station (Kalmár et al., 2019) in order to take into account the effect of low-velocity sedimentary layers on the receiver func-KALMÁR ET AL.
10.1029/2020JB021309 6 of 24 tions. Finally, we performed the H-Vp/Vs grid search method uniformly at the 221 other investigated stations (Note that 54 SCP stations were already investigated in Kalmár et al., 2019). Figure 4 shows the H-Vp/Vs results for a station on hard rock and another on thick sediment. The effect of the thick, low-velocity sedimentary layer manifests itself as a low amplitude in the P arrival and large-amplitude Ps conversion at the bottom of the sedimentary layer that are difficult to separate between 0 and 2 s increasing with epicentral distance (Yeck et al., 2013). Furthermore, in many cases, the highest amplitude includes the effect of both the sediment and the first P arrival in this frequency range (Hetényi et al., 2006;Zheng et al., 2005). We list the H-Vp/Vs analysis parameters and all results in the electronic supplement. These results serve as a first order reference for all stations, but cannot be considered as best estimates because of the assumptions made on the average Vp values, and because of a single Vp/Vs value for the entire crust.

CCP Migration
In order to assess the Moho depth results of the H-Vp/Vs grid search method, we imaged this discontinuity with CCP migration method (Zhu, 2000) using a recent 1D local velocity model (Gráczer & Wéber, 2012). In the migration, one of the most important parameters was the sedimentary basin depth correction, because the selected velocity model contains only a thin sedimentary layer (Kalmár et al., 2019). Therefore, we used a Neogene basement depth map compiled recently from reflection seismic profiles and well data (Balázs et al., 2018). Owing to the sedimentary depth correction, the P-to-S conversions are located at the appropriate depth, as could be observed in the deeper areas of the sedimentary basin (center of the Pannonian Basin and Vienna Basin). The migrated sections are also corrected for the station elevation above sea level.

Inversion and New Interpolation Method
The main focus of our study was the nonlinear receiver function inversion at each station, in order to get the most accurate image of the depth and velocity conditions of the three interfaces in the crust from teleseismic converted waves. The nonuniqueness of the receiver function inversion was described by Ammon et al. (1990). Therefore we applied the widely used Neighborhood Algorithm (NA) method (Sambridge, 1999) that gives an ensemble of acceptable solutions. NA inversion was performed individually for each of the 221 stations. This method preferentially samples the well-fitting parts of the parameter space, and has proven to be efficient at multi-layer receiver function problems. The implemented NA software uses frequency domain deconvolution for the synthetic receiver function calculation from the each sample model with the Thomson-Haskell matrix method (Haskell, 1953;Thomson, 1950). In contrast, we used iterative time domain deconvolution on the observed receiver functions, this is leading to minor differences in waveforms. For the entire investigated area, we considered a three-layer crust (sedimentary layer, upper crust, and lower crust) above a half-space mantle. We tried to define the simplest model for the area as this is well justifiable from the rich prior knowledge of the area, as well as the clearly perceivable effect of sediments on the receiver functions. The inversion procedure searched the thickness and Vp/Vs ratio of the crustal layers, and the S-wave velocity of the layers at their top and bottom. The depth-steps of the S-wave velocity discretization are not varied in the NA software's forward model, and we only retrieve layer's velocities at their top and bottom.
Initial search ranges of Moho depth and Vp/Vs ratio have been set using our H-Vp/Vs grid search results. Other searchable parameter bounds, such as the base of the sediments, depth of the Conrad surface and Vp/Vs ratio, and the S-wave velocity of the layers at their upper and lower boundaries, were selected from previous studies (Janik et al., 2011;Ren et al., 2013;Simonova & Bielik, 2016;Šumanovac, 2010). The full set of search parameters are listed in the electronic supplement. Based on the above, the result models contain two (upper crust and lower crust, in this case the sediment thickness is zero) or three layers (in areas covered by sediment). The limits of the gray area represent the parameter search ranges in the model figures (Figures 6d, 7c, 7f and 7i). We focused on the resolution of Conrad and Moho depths, but we used either 0 or minimum 1 km thickness for the sedimentary/slow layer (0 in the mountainous areas, which are certainly not covered by sediments based on the geological map) in the NA inversion.
In all cases, inversions have been run for 400 iterations with 50 models tested in each step. In the subsequent iteration step, the best 30 results defined the parameter space to be resampled. We tested these general inversion parameters and we have found these values able to constrain the investigated discontinuities and velocities in the crust. The reliability of the results is represented as the misfit value of the receiver function fitting on the 30-s time window. This time range includes every necessary peak (P, Ps, PpPs) of the investigated receiver functions (Singh et al., 2016). Moreover, we investigated the properties of the best KALMÁR ET AL.  2,000 models (10% of all models), and found the median model is a very good representation of these set of acceptable models (Figures 6 and 7).
In order to take into account possible topography in the major discontinuities such as a dipping Moho, we plotted the receiver functions sorted by back azimuth (binned every 30° with an overlap of 5°) at each station. In the case when the P and Ps peaks arrived at almost the same time from every back azimuth direction, we used a single receiver function stack without move-out correction in the inversion, illustrated in Figure 6. However, at several stations, the peaks had quite variable delay times, therefore we decided to identify similar receiver functions and group them by back azimuth ranges for separate inversions at a same station (Figure 7). To identify the groups, we also used tangential receiver functions. This resulted in two and three stacks of receiver functions from various back azimuth directions. Thus, the effect of anisotropy and dipping boundaries were partly covered by back-azimuthal grouping and separate inversion for those bundles (the overall analysis of receiver functions at all stations for dip and anisotropy related signals revealed a few local occurrences, but not a general pattern). Furthermore, we have taken into account gaps in the data. This was relevant when the data set coverage was not uniform in back azimuth directions beneath the station. For each station we list the back azimuth range of inversion groups, data gaps, nonlinear inversion results and misfits in the electronic supplement.
The uncertainty of the median model was determined with Kernel Density Estimation (KDE) method fitted on the best 2,000 models (Botev et al., 2010;Bowman & Azzalini, 1997). The KDE works well in complex and random parameter space, dynamically estimates the probability density function to different parameters (Xu et al., 2015). With this approach we were able to determine meaningful uncertainty estimates even KALMÁR ET AL.  though the distribution of the inverted parameters showed a large variation (Figure 8). The investigated distributions are about 60% not perfect bell-shaped (Gaussian) and the remaining 40% are skewed (e.g., logarithmic) in this study. However, the KDE calculation underestimates the variance of a sample drawn from a Gaussian distribution (Xu et al., 2015). In order to compensate this, we determined the ratio of the KDE variances to that of a theoretical Gaussian distribution in a synthetic test of 2,000 models that resulted in a ratio of 1.56. We calibrated the KDE variance estimates by scaling them up with this ratio (yellow lines on Figure 8), so that the KDE estimate will match the variance of a normally distributed population. The scaling somewhat increased the KDE variance estimates, but they are still between the one and two-sigma range (Duputel et al., 2012). The electronic supplement lists the uncertainty values for all inverted parameters. Following the NA procedure described above, we obtained one or more 1D velocity model(s) at every station. In order to derive a 3D model from the 1D velocity models, we have developed a new geometric representation of the obtained models combined with the widely used natural neighbor interpolation method. This model is not directly 3D but is constructed by a much larger number of 1D models than classical interpolation of a single 1D model per station. Our aim is to step forward from previous practices by approaching a very simplified migration procedure, and thus obtain a new, 3D image of structures with less interpolation effects. However, it is important to mention that the method can be used with high reliability if the BAZ coverage is adequate in most of the stations, otherwise we cannot estimate the direction and the angle of the dip of the discontinuities. Furthermore, using isotropic medium for this imaging algorithm is certainly an approximation. We call this approach the Natural Neighbor Cone Interpolation (NNCI).
Natural neighbor because we use this popular interpolation method (3D out of 1D velocity models) at each discontinuity (Amidror, 2002). The method is based on Voronoi tessellation of a discrete set of spatial points (Sibson, 1981). The natural neighbors of any node are those in the neighboring Voronoi cells, those to which the node is connected by the sides of the Delaunay triangle (Belikov et al., 1997). In our special case, we determined pre-grid with 0.01° longitude and latitude resolution.

10.1029/2020JB021309
10 of 24 Figure 7. Receiver function waveforms and inversion results at the JAVC station. We grouped the data by back azimuth at this station. We determined two groups, the yellow square represents the first (120°-30°), the green shows the second (30°-120°) based on radial and tangential receiver function of the station. (a) Radial and tangential receiver functions sorted by back azimuth, binned every 30° with an overlap of 5°. Time is relative to direct P wave. The smearing of the Ps peak begins after 120° back azimuth. (b) Receiver function waveform fit of all data, observed (black), median model (red), and 2,000 (10%) best models (green) receiver function. (c) 1D velocity-model inversion results of all data. Inversion results for (d-f) the first group and (g-i) the second group. Figure  descriptions are the same as in Figure 6.
Cone because it is the 3D geometry that describes the volume covered by teleseismic rays beneath a station ( Figure 9). This new visualization algorithm of the obtained models can handle data gap in back azimuth coverage and thus interpolates only data that are present, from its proper position. The key parameters in this novel algorithm are the radius of each cone at the respective discontinuity. These radii are determined by the highest angle of incident S waves (32°) and the first Fresnel zone so that the raypath lies entirely within this calculated zone ( Figure 9). Furthermore, the calculation reduces the shadow zone beneath and between the stations and leaves less room to cover by the subsequent interpolation. Along the perimeter of the cones, the circles were divided every 30° in back azimuth, this determined the location of each of our fixed points (with the results of the inversion to assign).
The NNCI has enabled uniform management of the NA inversion results and spanned a proper base for spatial interpolation to produce detailed maps of the major discontinuities. We attached three maps showing the piercing points at the three-layer boundaries in the electronic supplement.

Results of Inversion and New Interpolation
The nonlinear inversion results and spatial redistribution using the NNCI scheme allowed us to create detailed maps of the major velocity discontinuities in the crust as well as thickness maps of the sediments, upper crust, and lower crust. These maps are the first to be based on uniform data processing at a dense network covering the entire Pannonian Basin. Previous receiver function studies (Hetényi et al., 2015;Kalmár et al., 2019;Stipčević et al., 2020) covered only relatively small regions of the Pannonian Basin and its surroundings.
Note that for the sake of better visual presentation the color scale is different for each map. We used the NA inversion results and the NNCI method to produce on all presented maps below.

Basement Map and Sedimentary Layer
The interpolated sediment thickness map is shown in Figure 10. The figure also shows the main strike-slip and normal faults with high offsets (Horváth et al., 2015) and the stations contributing to this study. Our interpolated sediment thickness map is in good agreement with previous compilations based on reflection seismic profiles and well data (e.g., Balázs et al., 2016Balázs et al., , 2018Horváth et al., 2015). Remarkable differences between the two maps are located in the deep basin area, where we image young and unconsolidated loose sediments. In contrast, the maps based on well data and seismic stratigraphic interpretation contain also strongly compacted sediments with higher Vp and Vs velocity in the deeper region (7-8 km) of the depocentres. In our results, the thickest unconsolidated sedimentary layers (5-6 km) are in the southeastern part of the Pannonian Basin and in the Danube Basin. The southwestern part of the Pannonian Basin and Vienna Basin show somewhat shallower (2-3 km) basement depth. As expected, the Eastern Alps and the Southern Carpathians lack any sedimentary cover. The typical numerical uncertainty of the sediment thickness from the KDE estimate was ±0.1 km for stations where a sedimentary layer was included, however this value is clearly less than what is resolvable at the highest applied frequencies (on the order of a kilometer).  . Methodology of Natural Neighbor Cone Interpolation, illustrated with JAVC station (not to scale). We determined three discontinuities (disco.) within the cone (Sedimentary (Sed.); Conrad (Con.); Moho). The radius of the respective cone (R Sed.;Con.;Moho ) is determined by the Fresnel zone (Fresnel Sed.;Con.;Moho ) and highest angle of incident S waves (ROC Sed.;Con.;Moho ). Along the perimeter of the cones, the circles were divided every 30° in back azimuth. The filled green (30°-120° BAZ direction) and yellow (120°-30° BAZ direction) circles represent models of groups of observed data (two groups in this case). If there were no data from a given direction, the circle is left blank. The approach was the same at every discontinuity and at every station. the typical uncertainty was ±0.2 km/s. The Vp/Vs ratio is between 1.60 and 1.80 in the sedimentary layer. Usually, the 1.60 values are located in the thin sedimentary regions and the values around 1.80 characterize thick sediments. The typical uncertainty was ±0.04 in this case, as it was a fairly well-resolved parameter. Figure 12 shows the depth of the Conrad discontinuity. The Conrad discontinuity is the shallowest (12-16 km) beneath the sedimentary basins and deepest (24-28 km) under the Eastern Alps and Southern Carpathians. The typical uncertainty of the depth of this discontinuity was ±0.6 km, this value contains error of the sediment thickness too. Figure 13 shows the S-wave velocities at the top and the bottom of the upper crust, and the Vp/Vs ratio in the upper crust. Apart from some local low-velocity spots (2.10-2.50 km/s), the S-wave velocities are quite homogeneous at the top of the upper crust (3.25-3.50 km/s) with ±0.08 km/s typical uncertainty. The same is true for the bottom of the upper crust where the S-wave velocity varies between 3.50 and 3.75 km/s. The average Vp/Vs ratio of the upper crust varies between 1.60 and 1.80. The typical uncertainty was very small (±0.02) for this parameter. Figure 14 shows the detailed depth of the crust-mantle boundary. The Moho discontinuity is the shallowest beneath the sedimentary basins (20-26 km), goes deeper below the Apuseni Mountains, North Hungarian Range, and Transdanubian Range belts (28-33 km), and it is the deepest beneath the Eastern Alps and the Southern Carpathians (40-45 km). The Moho depth is well resolved as its typical uncertainty was ±1.2 km, also contains the uncertainties of the upper two layers.

Moho Depth and Lower Crust
KALMÁR ET AL.

Upper and Lower Crust Thickness Maps
As the Conrad discontinuity does not appear to be a major discontinuity from the results presented above, we have further investigated differences between the upper and the lower crust. We calculated the S-wave velocity gradient of the upper (GradV 1 ) and the lower (GradV 2 ) crust by Equation 1 and presented them in map view (Figures 16c and 16d).
where Vs22 is the S-wave velocity at bottom of the upper crust, Vs12 is the S-wave velocity at top of the upper crust, Vs32 is the S-wave velocity at bottom of the lower crust, Vs31 is the S-wave velocity at top of the lower crust, Z Sed is the depth of the basement, Z Conrad is the depth of the Conrad discontinuity, Z Moho is the depth of the Moho discontinuity.
The velocity gradient across the upper crust shows near-zero or slightly positive values in most cases. In contrast, the lower crust has clearly positive gradients across the entire region. This gradient might be linked to the gradually increasing metamorphic degree with increasing depth of the lower crustal granulites (Török et al., 2013;Weiss et al., 1999). Since the velocity jump across the Conrad discontinuity is less than 0.2 km/s over most of the investigated area, we conclude that it is better characterized as an S-wave velocity gradient change between the upper and the lower crust. With this approach, we are also able to present here the first upper crust and lower crust thickness map of the Pannonian Basin and its wider region (Figures 16a and 16b). The upper crust is thinnest (11-14 km) beneath basins, while the thickest parts (20-25 km) are located under the mountains. In contrast, the lower crust is thinner than the upper crust, not only beneath the basin but also in a wider region: between 6 and 11 km under the Pannonian Basin and beyond, and relatively thicker (14-19 km) only when reaching the high mountains (Eastern Alps, Southern Carpathians) and the Bohemian Massif.

Vp/Vs Ratio Comparison
We Vp/Vs ratio in each layer separately. To be able to make comparisons, we calculated the mean Vp/Vs ratio in the crust using Equation 2   of H-Vp/Vs grid search method were very good initial parameters, but the inversion method refined these and we could produce more interpretable and finer map of Vp/Vs ratio values. Figure 18 shows the comparison of the Moho depth from the two most recent publications Horváth et al., 2015)  previous studies lacked uniform data coverage and did not have the spatial resolution of this study. Therefore, these maps show high homogeneity in the center of the Pannonian Basin (22-27 km) and in the mountain regions (35-40 km). The map of Horváth et al. (2015) presents thicker (basins area) and deeper (high mountains) Moho depth with more accurate values, but these results rely on their gravity background model. The Moho depth map from our H-Vp/Vs grid search method with natural neighbor interpolation presents local features and these values as first order reference, were good initial parameters in the inversion method. The inversion method with NNCI provides a much more detailed Moho map of the Pannonian Basin and its surroundings. in the Tisza-Dacia unit, where Miocene extension ultimately localized at the end of syn-rift times (Figure 1). In general, AlCaPa shows somewhat higher crustal thickness values than Tisza-Dacia or the Mid-Hungarian (Fault) Zone, but no sharp boundary can be observed between these units. This infers that Miocene extension and associated upper crustal brittle faults and lower crustal ductile deformation resulted in a similar mechanical equilibrium (Hetényi et al., 2015). The crust-mantle boundary is also better resolved in the Transdanubian Range (27.5-30 km), North Hungarian Range (27.5-33) and Apuseni Mountains area (30-35 km). Thinnest lower crust is found beneath the main sedimentary depocentres (5-10 km; Figure 16b). Furthermore, the eastern part of the Transdanubian Range that lacks upper crustal faulting has similar low lower crustal values suggesting lower crustal thinning during the extension of the Pannonian Basin.

Moho Depth Comparison
In general, the western margin of the Pannonian Basin toward its Alpine and Dinaridic margin shows lower values of crustal thickness gradients than along the eastern margins toward the Apuseni Mountains or Southern Carpathians (see approximate dips on Figure 14). The smoother transition in the west is in well agreement with previous inferences on the gradual extrusion of the initially overthickened Pannonian basement units by detachments and rotating low-angle normal faults associated by the lateral flow of the weak lower crust (Fodor et al., 2008;Horváth et al., 2015;Tari et al., 1999). On the contrary, the eastern basin margin was rather affected by a series of high angle brittle normal faults leading to a sharper Moho variation. Similar asymmetrical patterns of upper and lower crustal thinning in the depocentres and neighboring basement highs is inferred by reflection seismic sections and thermo-mechanical models (Balázs et al., 2018). Finally, we found that, the Eastern Alps and the Southern Carpathians exhibit the thickest crust were between 40 and 45 km, putting the Moho somewhat deeper than previous studies. The highest crustal thickness gradient is observed between the most extended Great Hungarian Plain part of the eastern Pannonian Basin and the Eastern and Southern Carpathians where Miocene collision resulted in a thick orogenic crust.

Conclusions
We performed the first comprehensive receiver function analysis in the Pannonian Basin and surrounding regions using the most recent data set (221 stations) available. We mapped the thickness of major crustal layers and determined their S-wave velocity and Vp/Vs ratios. We also estimated the uncertainty of each parameter based on inversion results. Our study is based on a relatively long time-span ( broadband waveforms with uniform automatic waveform processing and quality control procedures. The receiver functions obtained from a network of closely spaced stations allowed us to infer a 3D structural and shear-wave velocity model of the region. We have developed a new interpolation and visualization algorithm (NNCI), in order to image seismic features (including dip estimates) as accurately as possible. We plan to investigate anisotropy in more detail using the tangential component of the receiver functions (e.g., Bianchi et al., 2010;Link et al., 2020) in the future. However, to do that with high reliability, it is necessary to have uniform azimuthal coverage to extract anisotropy information from the receiver functions and at the moment we have sufficient coverage only at the permanent stations with a long operational history.
KALMÁR ET AL. The Conrad depth, upper crust, and lower crust thickness maps are the first for the Pannonian Basin region and will provide reliable constraints for geodynamical numeric modeling studies. Our results help the better understanding of the Pannonian Basin evolution history and other geological aspects. The obtained sedimentary layer thickness and Moho depth maps show good correlation with previous ones, and they are richer in details. We propose that the Conrad discontinuity is a change in velocity gradient between upper and lower crust, rather than a large velocity jump at the lower-upper crust boundary. The Moho depth map presents local variations with more finely and better resolved values than previous investigations. The Pannonian Basin is thinner (20-26 km) beneath basins area and thicker (40-45 km) under the high mountains at the periphery of the investigated region. The results are obtained from broadband data of seismological stations with uniform and well-proven processing procedures. Overall, the dense seismic network with large amount of quality-controlled data processed here allowed to infer a 3D structural and shear-wave velocity model of the region, which is a valuable new information for any region.

Data Availability Statement
Further information about the data availability see under the Data tab.