Lithospheric Structure of the Circum‐Pannonian Region Imaged by S‐To‐P Receiver Functions

The lithosphere‐asthenosphere boundary and mid‐lithospheric discontinuities are primary attributes of the upper mantle. The Pannonian region is an extensional sedimentary basin enclosed by collisional orogens. Here, we estimate the negative phase depth of S‐to‐P receiver functions to image the lithospheric thickness and other discontinuities with high resolution, based on the recent dense seismological broadband networks. The lithosphere‐asthenosphere boundary is relatively shallow (<90 km) in the Pannonian Basin system, and deeper (∼90–140 km) in the surrounding orogens, where average surface heat flow values are higher (120 mW/m2) and lower (50–70 mW/m2), respectively. The 1D and 2D common conversion point migration with 3D velocity model provide comparable but different resolution images beneath the wider region of the Pannonian Basin. We obtained deeper values in the Western (∼120 km) and Southern‐Carpathians orogens (∼135 km). Furthermore, we provide new information on the lithospheric thickness and its seismic properties in the eastern part of the study region (e.g., Apuseni Mountains (∼95 km), Eastern‐Carpathians (∼120 km), Moesian Platform (∼90 km) and Transylvanian Basin (∼85 km). The shallower negative phase depth can be interpreted as the lithosphere‐asthenosphere boundary beneath the Pannonian Basin system in agreement with its high heat flow values. In contrast, the deeper negative phase depth estimates in the colder surroundings can be interpreted as intra‐ or mid‐lithospheric discontinuities, when compared with local seismic tomography models. In this region, the correlation with heat flow implies that the observed negative phase depth is of thermo‐chemical or rheological nature.

The Pannonian Basin has a thin (<70 km) lithosphere over a wider area at a large distance from the orogenic arc.Published estimates of the LAB depth suffer from significant uncertainty (at least 20 km) in the study region (cf., Horváth et al., 2006;Plomerová & Babuška, 2010).Furthermore, due to technical advances of acquisition and processing, new methodologies are available today, such as S-to-P Receiver Function (SRF) analysis.In addition, in the past decade, significant improvements have been achieved in the coverage and quality of temporary and permanent seismological broadband networks of the ACPDR.This provides an outstanding opportunity to image the shallow upper mantle and crust with high resolution.In contrast to the Moho interface, the LAB is usually not as sharp and clear discontinuity from a seismological aspect.The accuracy and resolution of the LAB, however, can be greatly improved from more recent seismological data using the SRF technique (Hu et al., 2011;Kind et al., 2012).The application of the method is very common in areas where the station coverage is dense, because in this way we can provide a coherent regional image of the lithospheric thickness (e.g., Hu et al., 2011;Kind et al., 2020;Liu et al., 2020;Zhang et al., 2012).Some SRF studies have already been conducted in a wider environment of the Pannonian Basin, but these do not show a coherent picture and their resolution is typically poorer (Kind et al., 2017(Kind et al., , 2021) ) or cover only part of the region (Belinić et al., 2017;Klébesz et al., 2015).
The geophysical data which provided the first information on the lithospheric structure of the circum-Pannonian region are more than 20 years old.Determining the lithospheric thickness is fundamental because it plays a major role in many geodynamic and geochemical processes such as volcanism, tectonic movements and earthquake nucleation.The geodynamic processes also include the coeval occurrence of thin back-arc mantle lithosphere, back-arc extension, non-isostatic forearc subsidence, slab retreat and associated mantle flow patterns (Balázs et al., 2022;Horváth et al., 2015;Liptai et al., 2022;Qorbani et al., 2015Qorbani et al., , 2016)).Currently, the rather uncertain and poor resolution of the LAB determination hinders the refinements of these features and models.
Within the lithosphere, other detectable discontinuities have been reported mostly from old cratonic regions.The mid-Lithospheric Discontinuity (MLD) was discovered more than two decades ago (e.g., Thybo, 2006;Thybo & Perchuc, 1997;Yuan & Romanowicz, 2010).The MLD is a significant boundary often reported in thick and cold cratonic continental lithosphere.Despite ongoing research, the origin of the MLD and its relationship to the LAB remain a topic of current interest.The MLD is a relatively thin (from a few kilometers up to a few tens of kilometers thick), characterized by lower seismic velocity beneath the discontinuity and higher conductivity (Selway et al., 2015) and its depth is typically at ∼100 km but can vary between 60 and 160 km (Fu et al., 2022;Rader et al., 2015).The importance of MLD has been highlighted in relation with various tectonic processes such as the removal or destruction of cratonic roots (Wang & Kusky, 2019), displacement of shallow cratonic lithosphere (Wang et al., 2017), delamination of lower continental lithosphere or localized subsequent deformation events (Dunbar & Sawyer, 1988).The origin of the MLD is a timely and controversial topic which includes elastically accommodated grain boundary sliding (Karato & Park, 2018), channel flow (Yang et al., 2021), anisotropy (Selway et al., 2015;Wirth & Long, 2014); presence of amorphous carbonates (CaCO 3 ; Zhang et al., 2023), pyroxene and garnet crystallization (Liu et al., 2022), partial melting (Mierdel et al., 2007), fossilized melts and LAB (Rader et al., 2015) and break-down of hydrous minerals (Green & Liebermann, 1976;Kovács et al., 2021).
Here, we aim to provide better constraints on the depth estimates of the aforementioned upper mantle discontinuities (i.e., LAB and intra-lithospheric discontinuities, like MLD) based on the Negative Phase Depth (NPD) of the SRF (Abt et al., 2010) and compared them to previous geophysical studies.This will allow us to refine the evolution of the Pannonian Basin, such as how the extrusion and thinning of the originally thicker lithospheric units happened.Was it driven by subduction roll-back and related passive asthenospheric flow (Horváth et al., 2006) or Alpine collision or Dinaric slab roll-back driven active mantle flow (Faccenna et al., 2014;Horvath & Faccenna, 2011;Kovács et al., 2012).This natural laboratory allows us to study a deeper insight into the origin and relation between the shallower LAB and the deeper intra-lithospheric discontinuities and their role in recent geodynamic processes.A new integrated SRF analysis and interpretation could shed new light on less known geodynamic, geochemical, petrological and volcanological processes in the ACPDR and may have wide implications for similar continental back-arc basins undergoing inversion worldwide.

Stations
Waveforms for all available broadband three component seismic stations located in the ACPDR were collected from 14° to 30° East and from 43.5° to 51° North (Figure 1).The waveforms were recorded at 389 (155 permanent and 234 temporary) seismological stations.The permanent stations in the study area belong to several countries: Hungary, Austria, Czech Republic, Slovakia, Ukraine, Romania, Serbia, Croatia and Slovenia.Most of these networks are members of the Central and Eastern European Earthquake Research Network, CE3RN (Lenhardt et al., 2021).

Teleseismic Events
Figure 2 shows the teleseismic earthquakes between 55° and 85° epicentral distances and magnitudes larger than 5.5 that we selected for the SRF analysis.This epicentral distance is mainly to avoid post-critical incidence and to keep the rays within the mantle.The data set consists of 2,525 events from the USGS catalog between 01.01.2002 and 31.01.2022.We downloaded the broadband three component waveforms of these events recorded at the stations shown in Figure 1.Based on these, we supplemented the data set with ca. 2 years of data from PACASE, CBP, and SCP temporary projects as well as more than 3 years of data from the AlpArray temporary network (e.g., Figure 2b).The permanent stations used all available data since they entered into operations up until 31.01.2022(e.g., Figure 2a).The exact stations operating time range and events number can be viewed in Supporting Information S2.

Quality Control and SRF Analysis
The SRF analysis is an effective tool to study seismic discontinuities in the upper mantle (Abt et al., 2010;Bock, 1991;Fischer et al., 2010;Kumar et al., 2012;Sacks et al., 1979;Snoke et al., 1977;Yuan & Romanowicz, 2010).The method is based on the S-to-P conversion of the seismic waves, which can be determined by the deconvolution procedure (Hu et al., 2011;Kind et al., 2012;Yuan et al., 2006).However, the large amount of data becomes unmanageable without quality control (QC) and filters, so we implemented a strict automatic procedure to select the best possible waveforms.The three-component seismic data were resampled (0.05 s), demeaned, detrended, and tapered.Data were then band-pass filtered in the frequency band of 0.04-1 Hz.To perform SRF analysis we downloaded 630,210 broadband three-component waveforms.The first quality control (QC1) was performed in the time window −200 before and +200 s after the predicted S-wave arrival (with IASP91 global velocity model) on the filtered ZNE components (Figure 3).This QC step is a complex signal-to-noise ratio investigation shows the number of quality-controlled SRFs as a function of BAZ at station A262A (temporary).The location of the stations is shown in Figure 1.(Colavitti et al., 2022;Hetényi et al., 2018b;Kalmár et al., 2021) adapted to S waves.We defined 3 different independent criteria, calculating the various ratios of the signals: 1.For each event, only traces with a root mean square (rmsall) between 0.5 and 3 times the median of all traces of each component for that event were kept; 2. The maximum amplitude of the S peak (maxpk) was required to be 7 times more than the rms of the background (rmsbg); 3. The absolute magnitude of the S peak amplitude (lgthS) within maxpk was required to be 4 times larger than the maximum amplitude of the pre-S wave window (maxbg); If all three (ZNE) components meet the moderately strong criteria (1,2,3), we accept these waveforms (Figure 3).These steps within QC1 removed ∼50% of the data from the data set, leaving us with 317,628 waveforms.
For the SRF calculation we used a shorter time window, 120 s centered around the first S-wave arrival.This window should theoretically capture target discontinuity conversion signals from the lithosphere-asthenosphere boundary or shallower discontinuities (e.g., first negative and positive phase, Yuan et al., 2006).We rotated the accepted ZNE waveforms into the ZRT coordinate system according to the theoretical back azimuth.After that, we rotated the ZRT waveforms into the local ray coordinate system LQT (separating the P, SV, and SH components) based on incidence angles.The L component points in the propagation direction of the incident SV wave, thus primarily containing energy of the converted Sp phase and nearly zero direct SV wave energy when passing through homogenous horizontal layers.The Q component, on the other hand, is perpendicular to the L component and contains significant direct SV wave energy (Kind et al., 2012).We present an example in Supporting Information S1, based on the waveforms in Figure 3.As an additional step, we tested several more rotation methods based on Kennett and Engdahl (1991), Reading et al. (2003), and Svenningsen and Jacobsen (2007).A small change in the rotation matrix can be observed in the different methods.However, we preferred the most complex and well-resolved rotation based on Kind et al. (2012).
We calculated the raw SRFs using the iterative time-domain deconvolution (Ligorría & Ammon, 1999) with 300 iterations.The Gaussian width was 0.3 during the deconvolution.In this case, we deconvolved the L component with the Q component, yielding raw SRF opposite in time and amplitude to the more commonly used P-receiver functions (Kumar et al., 2010;Yuan et al., 2007).To make it more comparable to the actual Earth's structure, amplitude and time axes of the calculated SRF were reversed yielding, positive Sp arrival time comparable with the P receiver function studies.We used a 35 s time window in the final SRFs for the second quality control (QC2).The resulting SRFs were subjected to the QC2 (Figure 4).In this case, we highlighted the importance of positive (SpP) and negative (SpN) Sp conversions, which contain important information about the lithospheric thickness.Poor quality SRFs were rejected if the rms of the S signal (between 0 s and the end of the SpN) is less than three times the rms of the background noise (between end of the SpN and 30 s).These QC criteria yielded 65,490 high-quality SRFs spread unevenly among stations (see Supporting Information S2).

1D Migration
In order to interpret the resulting SRFs and place constraints on the lithospheric thickness and other possible intra-lithospheric discontinuities, we employed and compared the results from two approaches to convert from time-domain to depth domain.First approach was to use 1D migration of the stacked SRFs with the IASP91 (Kennett, 1991) global velocity model.
This migration method is used on the post stack SRFs.It is therefore essential to perform the moveout correction before the stacking procedure because the arrival time of the converted waves is a function of the incident angle, which is controlled by horizontal slowness.Horizontal slowness is determined by hypocentral depth and the epicentral distance of the earthquake.In this study, as we investigated only the S phase, the horizontal slowness is between 9 and 12 s/° corresponding to 55° and 85° epicentral distance.Based on these criteria, we used a fixed reference epicentral distance (70°) during the moveout correction.After that, the stacking of the SRFs is an important step because only a small percentage of the internal energy is transformed.Stacking boosts coherent small energy arrivals, highlighting, phases of interest (SpN in this study) from the obtained SRF.The stacked SRF representing the lithospheric structure around the seismological stations.
In the last step, we performed time-to-depth conversion with the selected velocity model (Figure 5), which serves as a depth value below the station.We tested the ak135 (Kennett et al., 1995) and other global velocity models, due to a lack of reliable local velocity models.We preferred IASP91 to maintain consistency with the reference models of previous studies (Belinić et al., 2017;Kind et al., 2021).The median value of the SpN peak was the accepted depth value for 1D migration (Abt et al., 2010;Hu et al., 2011;Selway et al., 2015).
We applied the 1D migration procedure and estimated the Negative Phase Depth (NPD) at 389 stations) (see Figure 6; Supporting Information S2).We introduced the generic term NPD here to refer to the depth in the lithosphere  where negative amplitude occurs.This choice avoids attaching direct geodynamic or geochemic connotations.
We show that the NPD can be interpreted in various and sometimes contrasting ways.These results provide a good estimate of the thickness of the lithosphere in the studied area, but the determinations are impeded by large uncertainties, especially at the stations where we managed to generate a small number of SRFs (see Supporting Information S2 and Section 4.3), which is why the obtained results must be treated and interpreted with caution.Bias may stem from not accounting for 3D path effects or uneven azimuthal distribution of sources for the incoming waves, which may result in pulse smearing (Figure 4c).The 1D migration provides a first order depth interpretation of SRFs in this study.

2D Common Conversion Point Migration
The 2D Common Conversion Point (CCP) migration (Zhu, 2000) is the best way to convert time domain SRFs to depth along the ray path of the Sp phase.The CCP method on SRFs is widespread and has been reliably used in many studies (Akanbi & Li, 2016;Kind et al., 2012Kind et al., , 2015Kind et al., , 2017Kind et al., , 2020Kind et al., , 2021;;Knapmeyer-Endrun et al., 2017;Kumar et al., 2012;Wang et al., 2016).We imaged the NPD with this method using the IASP91 global velocity model and 3D model of Piromallo and Morelli (2003) and Zhu et al. (2015).3D models provide more detailed results and slightly deeper NPD values than 1D models (∼9 km, Figure 7; Supporting Information S2).Also, after smoothing, the negative phase depth is localized in a narrower range in the case of 3D than 1D models.The 3D models estimated based on different seismic tomography studies are similar.However, we prefer the Zhu et al. (2015) model due to its higher resolution and more accurate depiction of sedimentary basins compared to the model of Piromallo and Morelli (2003).Migration was performed on selected cross-sections.We set bin sizes of 2 and 1 km on the horizontal and vertical directions, respectively.Since the station coverage is not completely uniform, we used a 50 km smoothing factor between the stations (Supporting Information S1 shows our parameter tests of the smoothing factor).This smoothing factor is 3 times lower than that used by Kind et al. (2017), which investigated the Central European upper mantle discontinuities (e.g., Lehmann and 410 km, etc.).This type of CCP migration is performed before SRF stacking; therefore, all seismic rays traveling within the cross-section contribute to the final image.
Based on these, we have the possibility to calculate the real piercing points (Crotwell et al., 1999) of the SRFs (the point where the Sp conversion takes place in the vicinity of the station).The piercing point locations for each SRF were determined based on the depth values obtained from the 1D migration.Afterward, using the CCP migration, it was possible to determine the new NPD for each SRF, which was used to produce the final map (Figure 8).The station coverage is better in the western part of the study area than in the eastern part, yet our study still presents the most comprehensive SRF analysis and lithosphere thickness to date in the complex circum-Pannonian region.
Figure 8 shows all the piercing points (65,490), which we used to determine the NPD map, the average NPDs around the station environment from the CCP method (picked by hand).This dense data coverage provides clearly more reliable and new results from hitherto unresolved areas (e.g., Carpathians; Moesian Platform) than the results obtained from 1D post-stack migration (see in Supporting Information S1 and Supporting Information S2).

Uncertainty Calculation
Quantifying uncertainty is a crucial component when interpreting SRFs as it reveals the limitations and interpretability of the results.During bootstrap  resampling, we randomly sampled the same number of SRFs with replacement from the original data set at each iteration before migration (Efron & Tibshirani, 1986).We optimized the procedure used by Abt et al. ( 2010) for our SRF data set (Figure 5).To estimate the uncertainty of the NPD depth, we calculated a standard deviation on the NPD from a total of 100 iterations.Since the emphasis in this study is on the depth of the lithospheric discontinuity (especially NPD), we use the ±2σ (95%) standard deviation as the uncertainty measurement for the NPD.
It is important to note that the actual distribution of the data may deviate from a Gaussian distribution, which may affect the magnitude of the uncertainty.We present the stacked SRF results at LTVH station in Figure 5.
For the negative phase pick, the uncertainty in depth is taken as the maximum (ε + 2σ) and minimum (ε − 2σ) distances from the pick on the calculated SRF.The negative and positive errors are equal if we have a Gaussian distribution.However, if the negative phase distribution is not perfectly normal, it will be larger in one direction than in the other (e.g., LTVH).The values of the calculated uncertainty plotted on Figures 6 and 8 can be found in Supporting Information S2 for all stations.For CCP sections, we reported the average uncertainty per station, not the error of all SRFs individually.

Negative Phase Depth Maps
The NPD results obtained from the 1D and 2D CCP migration approaches and the previously available LAB map (Horváth et al., 2006) for the Pannonian Basin and its surroundings are shown in Figures 9a-9c.As mentioned above, this map uses a wide range of various geophysical input data but was hampered by a sparse seismic station coverage (Figure 9a).The western and outside of Hungary part of the LAB map between 14°E and 19°E is based on Babuška and Plomerová (1988, 1992, 1993) tomographic study.Unfortunately, the geophysical source data used for the previous LAB map is not open access, but it seems to be more accurate within Hungary (e.g., Great Hungarian Plain) than beneath the adjacent regions (e.g., Transylvanian Basin; Carpathians; Eastern Alps).The LAB depth increases linearly beyond the central Pannonian Basin toward the orogenic rims.Moreover, this map bears significant uncertainties, where the average errors can reach up to ±10 km which must be considered during the geological and geodynamic interpretations.
The NPD map shown in Figure 9b was constructed based on 1D migration data from 389 seismological stations using the natural neighbor interpolation scheme (Sibson, 1981).The data points do not cover the study area homogeneously, thus there are certain regions where our determinations are less reliable (e.g., Transylvanian Basin, Eastern-Carpathians).However, our NPD estimates bear significantly smaller uncertainty (±6 km, Supporting Information S2).The NPD obtained from the 1D migration method closely resembles previous maps (e.g., Tari et al., 1999) beneath the study area, showing similar, relatively shallow NPD under the Great Hungarian Plain (∼55 km) and Vienna Basin (∼70 km).Our map shows improved new information on the NPD for the Transylvanian Basin (∼85 km), the area of the Western (∼120 km), Southern (∼135 km) and Eastern Carpathians (∼120 km), and the Moesian Platform (∼90 km).The NPD for the Bohemian Massif (∼110 km) correlates well with the previous LAB map, which can be attributed to the fact that in these areas the map strongly relies on the results of the seismic tomography model of Babuška and Plomerová (1988, 1992, 1993).
We also constructed another map of the NPD results obtained from the 2D CCP migration (Figure 9c).This map was produced with the natural neighborhood cone interpolation algorithm (Kalmár et al., 2021) and contains 65,490 measured data points (Figure 8).The piercing points show relatively good coverage throughout the entire area, the western part is especially well covered.The northern (e.g., Western-Carpathian) and eastern areas (e.g., Eastern-Carpathian) can be refined with the inclusion of additional data.This map shows generally shallower NPD values than the 1D map (Figure 9b).On all three maps, the NPD correlates well with the LAB depth beneath the Great Hungarian Plain.This area is characterized by a thin lithosphere, associated with high surface heat flow (>90 mW/m 2 , Lenkey et al., 2021).In the areas where the lithosphere is thicker and the surface heat flow is lower (<70 mW/m 2 ), NPD values obtained from 2D CCP migration are less likely consistent with the LAB, but is probably indicative of another shallower intra-lithospheric discontinuity.Similar sharp "anomalous" NPD can also be observed in cross-sections (Figures 7 and 10).Similar interpretations can be proposed for other SRF studies covering the region (Geissler et al., 2010;Kind et al., 2017).Long period magnetotellurics studies (Jones et al., 2010;Korja, 2007) also identified conductive shallow layers within the lithosphere which do not seem to coincide with the LAB in thicker and older lithosphere.
Overall, the calculated NPD results are suitable for making maps in both cases of 1D and 2D migration, and they must be interpreted together in order to get the most accurate lithosphere thickness and properties.From the results, it can be seen that the NPD values obtained from the 2D SRF CCP migration closely approximate LAB (Rychert et al., 2020) in areas characterized by higher surface heat flow (>90 mW/m 2 ).

Comparison of Results With Tomography and Previous Studies
In Figure 10 we present the 2D CCP migrated cross-section from the Bohemian Massif through the Pannonian Basin to the Moesian Platform (B-B′ in Figure 1), across the entire study area from the NW to SE.This 2D CCP profile perfectly outlines the lithospheric structure diversity of the area, crossing the orogens (e.g., Apuseni  (Babuška & Plomerová, 1988, 1992, 1993;Horváth et al., 2006;Tari et al., 1999 Mountains and the Southern-Carpathians), the basin areas (e.g., Vienna Basin and Great Hungarian Plain) and the controversial Vrancea slab.Furthermore, for a more complex interpretation, in Figure 10  The NPD depths estimated under the Pannonian Basin correlate well with negative velocity anomalies in tomographic images, whereas such correlation is less obvious under the surrounding orogenic rims and massifs (e.g., Bohemian Massif and Southern-Carpathians).This suggests that the 2D CCP migration NPD results in the areas with higher heat flow (e.g., Pannonian Basin) represent the LAB depth.In the high topography areas, NPD signals probably indicate the presence of another lithospheric discontinuity on the migrated SRFs.Beneath the Bohemian Massif, the Moesian Platform and the East European Platform, discontinuities at ∼100 km depth could be identified as MLD.A local asthenospheric upwelling adjacent to the Vrancea slab was inferred in many previous studies (e.g., Martin et al., 2006;Petrescu et al., 2020;Ren et al., 2013;Russo et al., 2005).Figure 10 shows a slight disturbance in the estimated NPD depth and local uplift.The NPD is slightly elevated in the NW forefront of the Vrancea seismogenic body, which is in line with previous studies suggesting asthenospheric upwelling, volcanism and a relatively elevated LAB in the SE Carpathians (e.g., Harangi et al., 2013;Tiliță et al., 2018;Tondi et al., 2009).Korja (2007) investigated the LAB depth beneath Europe based on electrical conductivity estimated with magnetotelluric methods.Similar studies have been conducted in the Pannonian region for a long time (Ádám & Wesztergom, 2001;Ádám et al., 1996), but delimiting the exact LAB depth is not trivial because the values obtained here change in a relatively wide zone (∼45-90 km) into the extensional Pannonian Basin with a thin lithosphere.Plomerová and Babuška (2010) conducted the first global tomographic study to investigate the LAB in the Pannonian region.They found that in areas with thick sedimentary cover, such as the Great Hungarian Plain and Vienna Basin, the LAB is located at a depth of approximately 70 km.This estimate is still considered to be a very accurate approximation of the LAB depth today.For the surrounding orogens, this European-scale model provides LAB values of ∼125 km.Such inferences (Babuška & Plomerová, 1988) were integrated into later reviews on the lithospheric structure of the Pannonian Basin.
The area was partially involved in the first European SRF study of Geissler et al. (2010), which estimated the LAB depth at 4 stations in the Pannonian Basin: PSZ (∼101 km), BUD (∼74 km), SOP (∼89 km), and PKSM (∼103 km).The 1D migration results from Geissler et al. (2010) show a deeper LAB than this study.We obtained ∼64 km for PSZ, ∼77 km for SOP, ∼64 km for BUD and ∼74 km for PKSM (in this study MORH station) with the 1D migration.These differences can be attributed to technological development (e.g., QCs, rotation) and the increase in quantity and quality of the available data.
Previous local SRF studies partially overlapping with our study region include the Northern Hungarian region (e.g., PSZ station), where Klébesz et al. (2015) determined a LAB depth of 65 km (±10 km).Our estimate of the LAB is also ∼64 km, but with a much smaller uncertainty (±4 km; Supporting Information S2) due to the larger, quality-controlled data set and considerably improved data processing routine.Belinić et al. (2017) focused on the Dinarides and its surroundings and similarly used 1D migration with the IASP91 velocity model.The LAB depths determined in the southwestern part of Hungary (∼65 km), in the foreland of the Eastern Alps (∼100 km) and on the northern edge of the Dinarides (∼80 km) agree well with our NPD estimates from the 1D migration routine (Figure 9b).
A somewhat more comprehensive study by Kind et al. (2017) provides a few 2D CCP migrated cross-sections in the Pannonian Basin region.These NPD values are similar to the ones we estimated, in the Pannonian Basin (shallow lithosphere) and also in the Southern and Eastern-Carpathians (thicker lithosphere), where the depth values most likely do not indicate the LAB.Therefore, the geological and petrological interpretation of these phenomena is essential for understanding lithospheric structure and evolution.

Geological Implications: Interpretation of the NPDs as Signals From the LAB or Other Lithospheric Discontinuities
The lithosphere-asthenosphere boundary is often defined based on geophysical properties, such as contrasts of elastic, thermal, electrical, or seismic properties in the uppermost mantle (e.g., Artemieva, 2009Artemieva, , 2023)).While traditionally the ca.1,300°C isotherm was used to define the base of the rheologically strong lithosphere, later studies proposed the ca.1,100-1,200°C isotherm based on the presence and boundary of the underlying melt-bearing mantle (e.g., Kovács et al., 2021;Rychert et al., 2020).A shear wave velocity drop of 2%-10% is usually associated with the presence of partial melts (Saha et al., 2018;Thybo, 2006), changes in deformation creep mechanisms (diffusion and dislocation creep, Aulbach et al., 2017;Fei et al., 2016), a high water content layer produced by fluid metasomatism (Aulbach et al., 2017).Fluids can migrate from the underlying subduction or thermo-chemical plume sources and can accumulate in a mantle layer with low temperature and pressure, increasing water capacity at a depth shallower than 100 km (e.g., Fu et al., 2022;Karato et al., 2015).It has also been proposed that carbonated silicate melts and pargasite dehydration might be connected to intra-lithospheric discontinuities (cf.Kovács et al., 2021;Niu & Green, 2018).Elastically accommodated grain boundary sliding (EAGS) is also proposed to account for the lower seismic velocities associated with the LAB (Karato, 2012).Jackson and Faul (2010) experimentally identified an EAGS driven relaxation peak (at 1 Hz) at ∼1300 K for dry and melt barren upper mantle material with 1 mm grain size.This implies that seismic velocity drop could also be explained without the presence of partial melt at temperatures close to what is usually assumed in the vicinity of the LAB.In line with these inferences, another recent study (Samae et al., 2021) proposed that grain boundary sliding could be enhanced at the glass transition temperature of amorphous olivine, which also contributes to dropping the viscosity and seismic velocity at the LAB.The glass transition temperature of amorphous olivine agrees well with the temperature 1300-1400 K typically present at geophysical anomalies associated with the LAB (Pollack & Chapman, 1977;Samae et al., 2021).
Other models attributed these variations in the geophysical properties at the LAB to solid-state processes such as grain size variations (Austin & Evans, 2007), varying redox conditions (Cline Ii et al., 2018), or anisotropy (Rychert et al., 2007).Our new results show that the NPD determinations accurately resemble the LAB depth in the central part of the circum-Pannonian region where the basin underwent large amounts of extension in the Miocene leading to the present-day thin crust and lithosphere and high surface heat flow (Balázs et al., 2022;Lenkey, 1999).
Beyond the Pannonian Basin, our NPD estimates generally do not coincide with previous constraints on the depth of the LAB.Instead, in most cases NPDs indicate a shallower layer within the high seismic velocity continental lithosphere at ∼100 km depth rather than the LAB, which should be located at deeper levels in the upper mantle (Figure 10).We argue that NPDs indicate intra-lithospheric discontinuities (cf., MLD) in areas with thicker (>100 km) and colder lithosphere (Bohemian Massif, East-European Platform, Moesian Platform).This interpretation also explains that the NPD depression increases the rims of the Pannonian Basin where the heat flow is lower, and consequently, the depth of the 1,100°C is deeper.If the LAB and MLD were mainly gradual thermal boundaries, the corresponding negative phase signal would be too low in some regions to be detectable on SRFs (Eaton et al., 2009).Accordingly, the LAB and MLD would rather correspond to chemical discontinuities or phase transitions showing larger impedance contrasts where S-P energy conversion is more efficient (Sodoudi et al., 2013).
The dichotomy between intra-lithospheric discontinuities and LAB in the study area is overprinted by more subtle variations in the depth of the NPD.The most important changes occur at the transition between the Pannonian Basin and the orogenic rims coincident with the transition from the intra-lithospheric discontinuity (orogenic rim) to the LAB (Pannonian Basin).This transition may be gradual or sharp and is associated with subvertical negative velocity anomalies in the nearby upper mantle (Figures 7 and 10).One of the sharpest changes can be detected at the transition from the Dinarides to the Pannonian Basin (Figure 7), where the NPD is abruptly elevated from ∼100 km to ∼60-70 km.This sudden offset may be explained by the delamination of the lower part of the Dinaric continental lithosphere (Balling et al., 2021;Belinić et al., 2017).In contrast, the transition of the intra-lithospheric discontinuity beneath the Western Carpathians toward the stable European Platform appears to be smoother (Figure 7).This may be because lower convergence values may have taken place along the Western Carpathians.Instead, there were larger horizontal strike-slip movements between the Alcapa tectonic unit and the European lithosphere, also supported by the lack of well-developed orogenic roots under the Western Carpathians (Grad et al., 2006;Hrubcová & Środa, 2015;Hrubcová et al., 2010;Szafián & Horváth, 2006).

Conclusion
In this study, we obtained new information on the lithospheric structure of the broader Pannonian Basin region including the surrounding orogens with a resolution that was not available before.This is due to the new developments in methodology as well as substantial increase in the number of broadband seismic stations in the region.Precise waveform processing and the SRF calculation procedure are essential to produce reliable results.We applied two independent quality control procedures and tested multiple SRF rotation and calculation methods.The SRFs were interpreted in two different ways, which resulted in a new NPD map of the study area.The obtained results were analyzed from both geophysical and geological perspectives and compared to the previous studies.Our results show that the LAB is the shallowest under the central part of the Pannonian Basin, where the heat flow is the highest.The thickest lithosphere is imaged at the periphery of the study area.We propose that in these areas with lower heat flow values, NPD is not consistent with the LAB interpretation.We identify this interface as another intra-lithospheric discontinuity.The geological and petrological interpretation provides a new picture of the lithosphere-asthenosphere system within the circum-Pannonian region, which can be generalized to areas affected by similar geodynamic processes.

Figure 1 .
Figure 1.Distribution of three component broadband seismometers of the circum-Pannonian region in Central-Eastern Europe (inset).Colored triangles represent different networks of seismic stations (bottom left corner).Labeled stations are shown or discussed in the manuscript.Profiles A-A′ and B-B′ are the location of the migrated cross-sections.The scale presents the elevation of the investigated area in the bottom right corner.Abbreviations: Apuseni Mountains (AM); Bohemian Massif (BHM); Dinarides (DI); Eastern Alps (EA); Eastern Carpathians (EC); East European Platform (EEP); Great Hungarian Plain (GHP); Moesian Platform (MP); Southern Carpathians (SC); Transylvanian Basin (TB); Vienna Basin (VB); Vrancea seismogenic zone (Vr); Western Carpathians (WC).

Figure 2 .
Figure 2. Epicenter distribution of the 2,525 teleseismic events used for SRF calculation in this study.The circum-Pannonian region is shown by the green rectangle.The events are colored according to hypocentral depths.The red rose diagram (a) shows the number of quality-controlled SRFs as a function of BAZ at the PSZ station (permanent).The blue rose diagram (b) shows the number of quality-controlled SRFs as a function of BAZ at station A262A (temporary).The location of the stations is shown in Figure 1.

Figure 3 .
Figure 3. Theory (a) and accepted example of the QC1 procedure at PSZ station.The QC1 was performed on the ZNE components and the waveform was accepted if all three components met the established criteria.The presented earthquake parameters are shown in the lower left corner.The meaning of the abbreviations can be found in the main text.

Figure 4 .
Figure 4. Example and theory (c) of the QC2 procedure at final stack SRF of the PSZ station.Note that the time axis has been reversed.(a) Figure presents the stack SRF at station PSZ without any QCs.(b) Figure presents the stack SRF at PSZ with only QC1.(c) Figure presents the stacked SRF at PSZ with QC1 and QC2.The meaning of the abbreviations can be found in the main text.

Figure 5 .
Figure 5. Example of the 1D migration at LTVH station.The depth of the negative phase is 45 km; this station is located in southwestern Hungary (in Great Hungarian Plain), where the heat flow is the highest (ca.110 mW/m 2 based on Lenkey et al., 2002, 2021) and the lithosphere is the thinnest beneath the investigated area.The figure also presents the theory of the uncertainty calculation, which is based on Abt et al. (2010).

Figure 6 .
Figure 6.Map of the estimated depth of the negative phase from 1D migration with IASP91 global 1D velocity model beneath the ACPDR.The up and down arrows inside the circle show the upper (ε + 2σ) and lower (ε − 2σ) limits of the uncertainty values in the case of 1D migration.In this case, the colors of the arrows represent the uncertainties, the size is uniform everywhere.The thin lithosphere is located beneath the central Pannonian Basin and extends toward south and the thicker part located beneath the surrounding orogens (e.g., Carpathians, Eastern Alps, etc.).The values shown on the figure can also be found numerically in Supporting Information S2.

Figure 8 .
Figure 8. Piercing points distribution (gold dots) and average NPD for the stations (colored circles) from 2D CCP migration method with 3D velocity model of Zhu et al. (2015).The up and down arrows inside the circle show the upper (ε + 2σ) and lower (ε − 2σ) limits of the uncertainty values in the case of 2D CCP migration.In this case, the colors of the arrows represent the uncertainties, the size is uniform everywhere.Due to denser station distribution, the western part of the study region (between 14°E and 20°E) shows a better station coverage compared to the eastern part.The values shown on the figure can also be found numerically in Supporting Information S2.

Figure 7 .
Figure 7. Migrated SRF cross-section along the A-A′ line (shown in Figure 1) with (b and e), IASP91 global velocity model; (c and f), 3D model of Piromallo and Morelli.(2003); (d and g), 3D model ofZhu et al. (2015).The top panel shows the topography above the profile and the location of the seismic stations (colors of the stations are the same as presented in Figure1).The left side (b, c, d) shows the raw CCP migrated results with three different velocity models as mentioned above.The right side (e, f, g) shows the smoothed (50 km factor) and interpreted CCP migration results with three different velocity models, as mentioned above.The orange circles (LAB) and stars (intra-lithospheric discontinuity) show the NPD beneath the profile, transparency of the symbols supporting easier interpretation.Abbreviations are the same as in Figure1.

Figure 9 .
Figure 9.Comparison of the lithospheric thickness maps from (a) Previous studies(Babuška & Plomerová, 1988, 1992, 1993;Horváth et al., 2006;Tari et al., 1999) (b) 1D migration results of the SRFs; (c) 2D CCP migration results of the SRFs with 3D velocity model of Zhu et al. (2015).The color scale is the same on each map, the black interpolated area is off the scale in the previous LAB depth map.The black continuous line (on panel b and c) indicates the well-resolved area based on fix points of the interpolations.The Supporting Information S1 contains the surface heat flow map of the study area (based on Fuchs et al. (2021)).
Figure 9.Comparison of the lithospheric thickness maps from (a) Previous studies(Babuška & Plomerová, 1988, 1992, 1993;Horváth et al., 2006;Tari et al., 1999) (b) 1D migration results of the SRFs; (c) 2D CCP migration results of the SRFs with 3D velocity model of Zhu et al. (2015).The color scale is the same on each map, the black interpolated area is off the scale in the previous LAB depth map.The black continuous line (on panel b and c) indicates the well-resolved area based on fix points of the interpolations.The Supporting Information S1 contains the surface heat flow map of the study area (based on Fuchs et al. (2021)).
we also show the surface heat flow distribution based on Fuchs et al. (2021), the P wave tomography anomalies of Ren et al. (2012) and the S wave tomography results of Timkó et al. (2022).

Figure 10 .
Figure 10.Migrated SRF cross-section and velocity perturbation anomalies of the P and S wave tomographic models along the B-B′ line.The top panel shows the topography above the profile and the surface heat flow curve and data points of the profile from the global heat flow database (Fuchs et al., 2021).The second panel shows the interpreted CCP migration results and the location of the seismic stations.The third panel presents the P wave tomography model of Ren et al. (2012).The fourth panel shows the S wave tomographic model of Timkó et al. (2022) based on the surface wave dispersion curve measurements.The orange circles (LAB) and stars (intra-lithospheric discontinuities) show the average NPD beneath the profile.Above the last profile, you can see the joint interpretation of the anomalies.The abbreviations of structural elements are the same as in Figure 1.