Active tectonics of the Circum-Pannonian region in the light of updated GNSS network data

The Pannonian basin is an extensional back-arc basin that has undergone neotectonic inversion and is currently shortening. The understanding and quantification of present-day deformation processes during this inversion are still incomplete. To this end, we investigate the active deformation of the Circum-Pannonian region via the interpolation of GNSS-derived velocity field and the derivation of the strain rate fields. For the interpolation of the velocity field, we use ordinary kriging, a strochastic interpolation method. Our results show that estimating a strain rate field that is virtually free of short-wavelength noise requires the scaling of the velocity uncertainties, i.e. assuming a minimum standard deviation of 1 mm/yr in our case. The deformation of the Circum-Pannonian region is defined by the 2–3 mm/yr, NNE-directed motion of the Dinarides, and by the 0.5–1.5 mm/yr, WSW to SSW directed motion of the eastern areas (European foreland, East Carpathians, South Carpathians, Transylvanian basin). These opposite-sense motions define a large-scale, on average NE-SW shortening and transpression-type deformation in the Dinarides as well as in the Pannonian basin, while the East and South Carpathians undergo regional N–S extension. Neotectonic structures generally show good agreement with the strain rate field, for example in the Dinarides, Eastern Alps, or in the western Pannonian basin. However, the presence of fault-parallel shortening or biaxial shortening along sinistral neotectonic structures in the central and eastern Pannonian basin show some discrepancy between current geodetic and observed neotectonic deformation. The vertical velocity field shows dominantly 100 and 1000 km wavelength signals, the former is probably related to the response of the Pannonian lithosphere-asthenosphere system to neotectonic basin inversion, while latter can possibly be explained by far-field subsidence patterns induced by the mantle response to melting of the Fennoscandian ice sheet during the current interglacial period.


Tectonic setting
The Pannonian basin is a major back-arc basin surrounded by the Alpine, Carpathian, and Dinaric mountain chains (Fig. 1). The Miocene formation of the Pannonian basin was facilitated by the extension of a thickened orogenic area: the Cretaceous-Paleogene nappe stacks of the Adria-derived ALCAPA (Alps-Carpathians-Pannonia) and the Europederived Tisza-Dacia "mega-units" (e.g. Horváth et al. 2006). Extension in the incipient basin was coupled with large-scale translations and rotations (CCW in AlCaPa, CW in the Tisza-Dacia) related to the lateral extrusion of orogenic basement units from the surrounding mountain belts (Horváth et al. 2006;Márton and Fodor 2003;Ratschbacher et al. 1989). These lateral motions were largely accommodated by the Mid-Hungarian Line (Csontos and Nagymarosy 1998), Mid-Hungarian Shear Zone (e.g. Fodor et al. 2005;Palotai and Csontos 2010), or Mid-Hungarian mobile belt (Koroknai et al. 2020;Lőrincz et al. 2002); a major strike-slip deformation zone that has been changing its dominant kinematics between dextral and sinistral according to the distribution of extension in the ALCAPA and Tisza-Dacia basement units (Balázs et al. 2017 and references therein).
Miocene extension across the Pannonian basin has reached 150-200 km Ustaszewski et al. 2014), resulting in a thin lithosphere (Hetényi and Bus 2007;Kalmár et al. 2021;Timkó et al. 2019) and related high surface heat flow and geothermal Fig. 1 Tectonic framework of the Circum-Pannonian region showing geographical locations, regional tectonic units that constitute the basement of the Pannonian basin (after Schmid et al. 2008) and active faults (compiled from Bada et al. 2007; Koroknai et al. 2020). Vi = Vienna, Ko = Kosice, Bu = Budapest, Za = Zagreb, De = Debrecen, Cl = Cluj-Napoca, Be = Belgrade. Purple dotted line delineates the Pannonian basin 1 3 potential (Békési et al. 2018;Horváth et al. 2015;Lenkey et al. 2002). Syn-rift extension in the Pannonian basin ceased around 9 Ma (Balázs et al. 2016;Matenco and Radivojević, 2012), and was followed by basin-wide neotectonic inversion (contraction), which is still active Horváth and Cloetingh 1996;Koroknai et al. 2020). Neotectonic inversion is driven by the northward motion and CCW rotation of the Adriatic microplate, with respect to Eurasia, while the halted rollback in the East Carpathians constrains the Pannonian lithosphere towards the east (Bada et al. 2007;Horváth 1993). The combination of these two geodynamic factors have resulted in the relatively slow neotectonic contraction of the Pannonian basin, characterised by decreasing rates of north-eastward tectonic motion towards the NE, with respect to fixed Eurasia (Grenerczy et al. 2005). This contraction is largely reflected in the reverse and oblique reactivation of pre-existing normal faults, related folding and the formation of transpressional fault zones Koroknai et al. 2020).

Research goals
A remaining challenge in the Pannonian basin is the quantification and analysis of presentday deformation processes. The scarce GNSS network of the region in the 2000's resulted in the interpolation of the velocity field based on a relatively crude dataset (Grenerczy et al. 2005), which did not enable the reliable calculation of strain rate components, given the small gradients in surface velocities in the area. More recent GNSS-based studies successfully demonstrated the potential of geodetic strain analysis from a global (Kreemer et al. 2014) or European scale (Piña-Valdés et al. 2022) perspective. However, a targeted analysis of the Circum-Pannonian region is still lacking. For these reasons, the current understanding of present-day deformation processes in the Pannonian basin largely relies on structural observations from field, seismic and well surveys Koroknai et al. 2020), which are yet to be compared to a detailed geodetic strain rate analysis.
Here we supplement these studies by analysing the tectonic signal in the measurements of a substantially updated (densified) GNSS network. Our goal is to study to what degree present-day surface velocities correspond to interpretations of neotectonic deformation from non-geodetic observations. We analyse the spatial scales at which the GNSS-derived velocities can provide reliable deformation rates. We interpolate the geodetic velocity field using local ordinary kriging and derive strain rate fields of the Circum-Pannonian region and discuss their relevance for regional tectonic processes.

GNSS dataset
The continuously operating reference station (CORS) networks are the main tools providing quantitative information for ongoing tectonic motions. In Europe, the EUREF Permanent Network Densification project (EPND) aims to integrate all national CORS networks on the product level (Kenyeres et al. 2019) and to deliver dense, high quality station velocity information. The GNSS data is routinely processed by national Analysis Centres (AC) strictly following standardized EPN GNSS processing guidelines and the results are delivered in the standard SINEX (Software Independent Exchange) format. The daily/weekly position solutions coming from each ACs are intercompared, cleaned and then merged on the weekly SINEX level generating a unique solution series. Afterwards a large adjustment process is launched including all merged SINEX files and delivering a multi-year 1 3 position and velocity solution. The EPND products are a compound of 26 Analysis Center products based on Bernese and GAMIT processing. All combination processes are using the CATREF software (Altamimi et al. 2007). The CATREF software is an internationally recognized tool to combine space geodetic analysis results. It is used for generating the consecutive global ITRS (International Terrestrial Reference System) realizations. The EPND product series are using input weekly position solutions starting from GPS week 1500 (October 2008) and these are updated in each year. The actual solution series used in this study involves weekly SINEX data up to GPS week 2150 (March 2021), corresponding to time series length of maximum 13 years. The EPND database involves some 4000 stations that ever appeared in any of the AC inputs, however the velocity product published has only some 2400 entries as stations with time series shorter than 3 years and stations with apparent velocity issues are removed from the public database. The raw velocity product is filtered by an unsupervised classification filtering approach (Magyar et al. 2022). The automatic procedure marks/removes stations with time series shorter than 3 years and then searches for apparent velocity outliers in a given neighborhood around each station. This remaining velocity field is then also reviewed by a final manual filtering. Further details can be found in the dedicated EPND website (https:// epnd. sgo-penc. hu). The velocity uncertainty estimation in the prime CATREF solution is based on the least squares approach, however as it provides too optimistic uncertainties, we use the HEC-TOR software (Bos et al. 2013) to derive more realistic uncertainty estimates based on a standard power-law with white noise approximation. As the velocity estimation is coming from a single, consistent, large scale adjustment process we retained CATREF's velocity estimates, which we then complemented with HECTOR uncertainties, estimated station by station and component by component. We discard horizontal velocities with uncertainty standard deviations larger than 2.5 mm/yr.
In the present analysis we have used a subset of the actual EPND solution, D2150 supplemented with an independent velocity resource computed by the Nevada Geodetic Laboratory . NGL velocities are only used where no EPND data was available and for time series longer than 3 years. The NGL velocities are stemming from a global daily PPP (Precise Point Positioning) GNSS data analysis. Both velocity sources are expressed in the same reference frame, IGS14, which we then have transformed into ETRF2000. We have compared the agreement of the two independent velocity sources at the 446 points for which we have NGL and EPND solutions. The histogram (Supplementary Fig. A1) clearly shows, that there is no systematic bias between the two solutions, the ETRF2000 realization of both solutions is fully compatible. The horizontal velocity differences are small, with a mean of 0.01 / − 0.06 and standard deviations of 0.42 / 0.53 mm/ yr for north and east components, while there is slightly less agreement for the vertical component where the differences have a mean of 0.0 and standard deviation of 1.29 mm/yr, respectively (Suuplementary Fig. A1).
The velocity uncertainties derived from PPP are considered in general higher than the values from a DD (double difference) solutions we use in EPND. We compared the velocity uncertainty estimates at the common stations over the study area. Although the peak frequency of the computed NGL/EPND uncertainty ratio was around 1 (Supplementary Fig. A2) for all 3 components, the median uncertainty ratio is 1.71 / 1.47 / 2.29 for north, east and up, respectively, indicating that the NGL uncertainties are in most cases larger than those from EPND. Figure 2 shows the station distribution and their velocity estimates. The horizontal velocities are expressed in the ETRF2000 geodetic reference frame, where we have subtracted the Eurasian Euler pole estimate from ITRF2014 (Altamimi et al. 2016) to put the Horizontal velocities are expressed in the Eurasian-fixed reference frame ETRF2000. Error ellipses denote 1-sigma standard deviation uncertainties. The dotted box represents the area of interest for which we estimate strain rates velocities in an Eurasian-fixed reference frame. The published EPND velocities reflect the steady interseismic velocity field, in the timeframe of 2008-2021. In our study area no significant earthquakes occurred that may have affected the interseismic velocity estimates.
A first glance at Fig. 2a reveals that the GNSS velocity vectors with respect to Eurasia are small in the Pannonian Basin (as well as the areas to the north and northwest), generally below 1 mm/yr, without a clear common azimuth. Relatively high (up to 5 mm/yr) NNE to NE pointing velocities prevail in the northeastern Dinarides, with a strong decay when moving towards the Pannonian Basin. East of the Carpathians a westward motion is dominant, with a transition towards almost purely southward motions west and especially south of the Southern Carpathians. Vertical velocities show general subsidence in the region with respect to the ITRF reference frame (Altamimi et al. 2016), which is fastest in the SE-Pannonian region in our study area (2-3 mm/yr, Fig. 2b).

Velocity field interpolation
To estimate a continuous strain rate field, we interpolate the GNSS-derived point velocities using local ordinary kriging. Ordinary kriging is a stochastic weighted averaging method, that uses statistical information of the data to estimate quantities at unobserved locations (Wackernagel 2003). This statistical information is described by a correlogram that specifies the variance and the covariance/correlation of the velocities as a function of distance. Generally, nearby velocities are similar, and correlation between velocities reduces with increasing distance. To determine local correlograms (one for each interpolation point), we first empirically determine the local variance, where we determine correlation between all observations as a function of distance. We apply weighting of local data using gaussian kernels (Machuca-Mory and Deutsch 2013). The kernels allow for assigning larger weights to nearby observations such that the empirical statistics can be determined at each interpolation point. We use a half-width at half maximum of 250 km for the gaussian kernels, which implies that the local statistics are mostly determined by velocities within a 250 km radius of each interpolation point. To each set of local empirical variance and correlations we fit an exponential correlogram ρ, describing the correlation as function of distance h. Because the velocity field is non-stationary, as shows from the high velocities and gradients in the Dinaric mountains, and small velocities in the Pannonian basin, the local correlograms allow for spatial variability in the correlation and variance of the velocity field. We apply kriging separately to the east and north components of the velocity field. The interpolation estimate for unobserved location 0 is based on kriging weights w i , using the velocity z at observation point x i : Weights w can be obtained by solving the ordinary kriging equation (Montero et al. 2015):ẑ with covariance matrices C ij and Lagrange multiplier μ, which ensures that the total weight w sums to 1. We prescribe covariance is using an exponential correlogram as: with local variance C(0), correlation ρ as function of distance h between observation site i and j, velocity uncertainty variance σ 2 . In the exponential correlogram that we fit to the empirically determined local velocity correlation as function of distance, the range parameter r describes the decay of correlation. Relative sill value s in the range [0 1] and describes the degree of decorrelation at short distances. The relative sill can be described as a nugget variance σ 0 2 that describes the part of the variance that is local only (i.e. velocity errors, or local effects that affect the velocity only at a particular site): Supplementary Figs. B1 and B2 in the appendix show the local variance, range and nugget variance parameters that we have estimated from the empirical covariance of the GNSS velocities, and that we have subsequently used for the local ordinary kriging.

Strain rate determination
We calculate strain rates in a local 2D plane from the continuous description of the interpolated velocity field, by determining the horizontal part of the velocity gradient tensor L (x: east, y: north): The flat plane approximation is valid since contributions from radial deformation (at 1-2 mm/yr level) to horizontal strain rates are 2 orders of magnitude smaller compared to the contribution from gradients in east and north velocity directions (e.g. Savage et al. (2001)). We determine the gradients in the above equation by taking the velocity component difference at longitude and latitude : and where and are small steps (equal to half the grid width and height) in longitude and latitude, and ∆x and ∆y are distances on the sphere.
The strain rate tensor, or rate of deformation tensor, D can be derived from L as: and the spin tensor W is: where ω is the average angular velocity (counter-clockwise).
The strain rate tensor D is dependent on the choice of base (reference system orientation). Expressed as principal strain rates, we derive a base invariant description of the rate of deformation, independent on the choice of reference system. We denote the principal strain rate with the largest and smallest magnitude as ̇1 and ̇2 , respectively. Based on the relative magnitude of ̇1 and ̇2 we can qualitatively describe the strain rate (Kreemer et al. 2014), and we follow Broerse et al. (2021, their Eq. 43) to characterize strain rates as shortening, shear, or extension. Here, shortening and extension are distinguished as cases where one of the two principal strain rates is dominant in magnitude; the sign of the ̇1 is negative for shortening, and positive for extension. For shear both ̇1 and ̇2 have similar magnitudes, while its sign is opposite. Based on strain rates only we cannot discriminate between pure shear (as in extrusion, without rotation) or simple shear (including rotation, such as strikeslip). Furthermore, we use the second invariant of the principal strain rates to describe the general magnitude of the deformation rate.
We also determine the dilatation rate, which shows the rate of area decrease (negative) or increase (positive):

Consideration of velocity uncertainties
The covariance matrices C ij in the kriging equation include the estimated velocity uncertainties of each single GNSS-site. As is visible in Fig. 2a, in the center of the Pannonian basin, estimated velocities are generally higher than the 1−σ uncertainties, while from a first glance, a coherent motion that is shared by multiple stations is lacking. This hints at uncertainties that underestimate the degree to which local velocities represent regional velocities (at the scale of the typical inter-site distance). As strain rates are based on spatial velocity differences, uncorrelated noise in the velocity field is amplified, and errors in single GNSS-velocity estimates can lead to spurious strain rate patterns in between observation points. The spurious strain patterns manifest as strain rate sign changes in between neighbouring GNSS sites, leading to quasi-checkerboarding patterns. We therefore state that reliable strain rate patterns (an area with a single sign in strain rate) should at least be supported by three subsequent observations. As the velocity field is relatively homogeneously sampled, we can use the mean distance between neighbouring sites h. Here we define neighbours as the natural (direct) neighbours in a Delaunay triangulation. For strain rate patterns to be supported by at least three observations, the minimum wavelength T of strain rate patterns should be at least 4 times h (i.e., a half wavelength will span three sites), thus T/2 > 2 h. We vary the minimum velocity uncertainty σ min , as the velocity uncertainty largely determines (by the so-called nugget effect) whether only direct neighbouring observations receive a large weight, or whether observations lying behind the direct neighbours also receive a significant interpolation weight. In Fig. 3 we show that for increasing σ min the strength of small scale strain patterns diminishes. We test the effect on the dilatation rate, which is an approximately linear function of the principal strain rates. In the lower right panel of Fig. 3 we compare the, Fig. 3 Effect of minimum velocity uncertainty on strain rate patterns. The dilatation rate has spurious patterns in between observation sites (black dots) when there is no imposed minimum velocity uncertainty (upper left panel). These spurious patterns decrease in amplitude when the minimum velocity uncertainties increase. The lower right panel depicts the radially averaged power spectrum density (PSD) of the dilatation rate. Wavelength T/2 = h represents patterns with spatial scales equal to the mean distance between neighbouring GNSS sites. Wavelength T/2 = 2 h represent signals that follow from at least three subsequent site velocities, and which should be seen as the minimum resolvable wavelength for strain rate patterns when estimated velocities have significant uncertainties radially averaged, power spectrum density of the dilatation rate to the spatial scale of the observations. For σ min = 1 mm/yr the strain rate wavelengths below the critical threshold of T/2 = 2 h have significantly reduced, while the long wavelengths still retain much of their power. At even higher minimum uncertainty levels, the short wavelength signal is even further attenuated, but at the expense of the long wavelength strain rates. The general patterns of extension and shortening remain in place when we increase the minimum uncertainties from σ min = 1 mm/yr to σ min = 2 mm/yr, which implies that the general large scale deformation patterns with half wavelengths larger than twice the mean inter site distances are robust features. As the solution with σ min = 1 mm/yr is has dominantly signals above these length scales, we can confidently rely on the main features in the strain rate field at this minimum observation uncertainty level. In the remainder of this manuscript we thus assume that the minimum uncertainty of GNSS-derived horizontal velocities of individual sites is 1 mm/yr. As we do not derive spatial gradients from the vertical velocities, we do not apply an uncertainty threshold. We have determined residual velocities (GNSS-derived minus interpolated velocities) for the applied solution, shown in the Appendix for horizontal (Supplementary Figs. B3, B4) and vertical (Supplementary Figs. B5, B6) components. The mean of the normalized residuals (the ratio of residual/observation standard deviation) of the solution is very close to 0, showing that no significant biases arise from the interpolation method. Furthermore, the root mean square (RMS) of the normalized residuals is around 0.75 for the horizontal components, and 1 for the vertical velocity, which implies that the residuals are on average within 1 standard deviation of the input velocity uncertainties (for horizontal velocities at minimum 1 mm/yr).

Interpolated velocities
Horizontal velocities determined by the GNSS measurements with respect to the stable Eurasia plate show a characteristic NNE directed motion in the SW part of the study area (Dinarides and its transition towards the Pannonian basin) as well as a characteristic SSW directed motion in the SE part of the study area (Transylvanian basin and south-southeast Carpathians (Figs. 2a,4a). The interpolation of the velocities reinforces these observations, showing an up-to 3 mm/yr, NNE directed motion of the Dinarides, and an up-to 1.5 mm/ yr, SSW directed motion of the Romanian Carpathians (Fig. 3a). The Pannonian basin and the Western Carpathians in Slovakia are characterised by very low, 0.1-0.5 mm/yr velocities, with the values being somewhat higher in the SW part of the Pannonian basin. The effect of assuming a minimum velocity uncertainty σ min of 1 mm/yr is that the interpolated velocities in the Pannonian Basin are slightly smaller than the GNSS velocity vectors. At the same time the azimuth of the interpolated velocities is much more gradually varying compared to those from the GNSS velocities. The direction of this slow motion is E-ENE in the western and central Pannonian basin, WNW in the Western Carpathians, and SSW in the eastern Pannonian basin (Fig. 4a). Residuals between input velocities and interpolated velocities are unbiased and within the used observation uncertainties, as Supplementary Figs. B3 an B4 show.
Vertical velocities indicate a wide-spread subsidence in the Circum-Pannonian region, with the eastern and southern parts of the Pannonian basin experiencing the most rapid subsidence of 1.5-2 mm/yr (Fig. 4b). The margins of the basin, such as Fig. 4 Interpolated horizontal a and vertical b velocities. a Colors denote the horizontal velocity magnitude, circles show the GNSS-derived velocity magnitudes, all with respect to stable Eurasia. The blue vectors represent the GNSS-derived velocity vectors with 1−σ uncertainty ellipses. Black vectors represent the interpolated velocity vectors. We use a minimum uncertainty standard deviation (σ min ) of 1 mm/yr for each individual site velocity for the stochastic interpolation (local ordinary kriging). The parameters of the estimated local correlograms that we use for the interpolation can be found in the appendix. b Colours denote the vertical velocity magnitude for the input original data (circles) as well as for interpolated areas. Black lines denote active faults from the GEM Global active faults database and country borders the North Hungarian Range, the Apuseni Mountains, or the junction area between the Dinarides and the Eastern Alps show near 0 vertical velocities (either very slow subsidence or very slow uplift), hence representing a substantial relative uplift with respect to the large depocenters (Fig. 4b). Supplementary Figs. B5 and B6 show that all interpolated vertical velocities are within the uncertainties of the GNSS vertical velocities. Only in the Dinarides we see that the interpolated values are somewhat lower than the observations (but still well within the uncertainties), owing to the relatively large reported uncertainties in this region. Figure 5 shows the strain rate estimates of our reference model. The second invariant of the strain rate tensor shows the magnitude of strain rate over the area (Fig. 5a). The Croatian Dinarides deform the fastest with up to 20 nstrain (nanostrain)/yr, while the rest of the area typically deforms at rates of 1-10 nstrain/yr. The North Hungarian Range, the Western Carpathians and the European foreland show the lowest strain rates with values generally close to 1 nstrain/yr (Fig. 5a). The highest strain rates within the Pannonian basin are found in the SW (along the Slovenian-Hungarian and Croatian-Hungarian borders), reaching 6-10 nstrain/yr.

Strain rate analysis
The rotation rates are ranging between ca. 6 and -6 nrad (nanoradian)/yr, where positive values are assigned to CCW, and negative values to CW rotations (Fig. 5b). The NW of the area (Eastern Alps, Western Carpathians, European foreland) is characterised by slow (ca. 1 nrad/yr) but consistent CCW rotations. The Romanian Carpathians also rotate CCW albeit slightly faster (1-4 nrad/yr). In contrast, large parts of the Dinarides in Croatia, Bosnia and Herzegovina and Serbia together with most of the the Apuseni mountains and the Transylvanian basin are characterized by consistent CW rotations. The Hungarian part of the Pannonian basin is characterised by very slow rotations that vary between 1 and − 1 nrad/yr. The SW corner of the area (External Dinarides) show relatively fast CCW rotation.
Dilatation rates provide an indication for ongoing area extension (positive values) or contraction (negative values). Significant positive dilatation (ca. 5 nstrain/yr) characterises the Romanian Carpathians and the Transylvanian basin (Fig. 5c). In addition, the SW corner of the study area also shows positive dilatations due to the slower NNE motion of the coastal areas compared to the more internal parts of the Dinarides (Figs. 2a, 4a, 5c). The Western Carpathians and the North Hungarian Mountains show slower positive dilatation (ca. 1 nstrain/yr). In contrast, the Pannonian basin shows negative dilatation values of 1-5 nstrain/yr, while the mountain ranges in Slovenia, central and NE-Croatia and Bosnia show a faster, 5-10 nstrain/yr negative dilatation rate. Figure 5d shows the type of strains and the orientation of principal axes (shortening and extension) over the area. The junction area of the Eastern Alps, Dinarides and the Pannonian basin is dominated by NNE-SSW to NE-SW directed shortening and transpression (shortening and subordinate perpendicular extension), which is fading away in intensity inside the basin, as highlighted by the decreasing size of the principal strain rate symbols. Still, the majority of the Pannonian basin is characterised by shortening and transpression (transpression in the sense of deformation style, not stress regime), with the shortening direction changing from NNE-SSW in the western areas

Large-scale geodynamic drivers and regional deformation pattern
The new results show that the Pannonian basin is dominantly characterized by shortening and transpression (Figs. 5d, 6). This arises from the 2-3 mm/yr NNE motion of the Dinarides combined with the 1-1.5 mm/yr WSW to SSW motion of the eastern areas in Ukraine and Romania (Fig. 4a). These motions compress the Pannonian basin (Bada et al. 2007), which in response undergoes contraction (Fig. 5d). A critical condition for creating such a stress and strain rate field is the continuous NE directed push from the Adria microcontinent (CCW rotation on a larger scale, which is outside the boundaries of our interpolated area, Bada et al. 2007;Grenerczy et al. 2005;Márton and Fodor 2003), dictated by the overall convergence of Africa and Eurasia. Due to the Adria push, the area adjacent to the remaining part of the Adria microplate (the Dinarides) is characterized by the highest velocities and strain rates in the Circum-Pannonian region (Figs. 4a, 5a). Decreasing strain rates towards the inside of the Pannonian basin show that the effect of Adria push gradually decreases as more and more shortening is accommodated in the basin (Figs. 5a, d): in the central parts of the Pannonian basin and at the North Hungarian Range strain rates are already close to negligible. These results generally support and supplement earlier findings regarding the overall motion of the region (Grenerczy et al. 2005) as well as the driving mechanisms of neotectonic compression and shortening (Bada et al. 2007;Horváth 1993).

Fig. 6
Principal strain rate axes of the Circum-Pannonian region plotted over the tectonic framework and the neotectonic fault network. Fault kinematic interpretations rely on structural, reflection seismic, and focal mechanism solution data Koroknai et al. 2020;Wéber 2016). Converging arrows (red) indicate the shortening direction, diverging arrows (black) indicate the extension direction. The size of the symbol is proportional to the strain rate magnitude (corresponding to the magnitude shown on Fig. 5a; i.e. larger symbols mean higher strain rates). The dashed line denotes the Pannonian Basin. (Color figure online) In addition, our results show consistent SW-directed motion in the region of the Southern Carpathians (Fig. 4a) as well as further towards the south (Fig. 2a), where the southward component increases towards the Aegean Sea. This motion is related to the southward rollback of the Hellenic slab, which provides the space for southward motions and extension in the Aegean domain since Eocene times (e.g. Brun et al. 2016 and references therein).

Dinarides and SW-Pannonian basin
Our results show dominantly shortening and subordinate perpendicular extension (transpressive deformation) in the north-western Dinarides (N-S to NNE-SSW shortening axis, Fig. 6) and in the SW Pannonian basin (NE-SW shortening axis). The shear component (transpression) is more pronounced in the Dinarides, which fits well to the reported activity of major strike-slip faults that formed as response to the indentation of Adria in the Southern Alps (Picha 2002;Ustaszewski et al. 2014;van Unen et al. 2019;Xiong et al. 2022). The shear component diminishes towards the Pannonian basin, which is characterised by almost pure shortening. The Zala basin in the SW Pannonian (Fig. 6) exhibits several major neotectonic fold structures (e.g. Budafa anticline). The fold axes are generally oriented parallel with major Miocene normal faults (E-W to ENE-WSW), some of which were reactivated as reverse faults during the neotectonic inversion Horváth 1995;Koroknai et al. 2020). The orientation of the folds suggest an N-S to NNW-SSE shortening (and compression) direction in the Zala basin (Koroknai et al. 2020), which deviates ~ 25° from the NNE-SSW to NE-SW shortening directions derived from the GNSS data (Fig. 6). In our view, the formation of the folds was constrained by the inherited, reactivated fault structures. Consequently, the observed orientation was not only defined by the acting stress field, but also by structural inheritance, resulting in an orientation slightly oblique to the overall compression and shortening directions.

Mid-Hungarian mobile belt
The Mid-Hungarian mobile belt (in the sense of Koroknai et al. 2020) is interpreted to be a ca. 75 km wide, ENE-WSW trending, currently sinistral strike-slip shear zone (Fig. 6). Our results highlight that the area of the shear zone is characterized by moderate to very low strain rates, decreasing from 10 nstrain/yr to 1-2 nstrain/yr from the SW to the NE (Fig. 5a). This shows that deformation in this formerly prominent shear zone is currently not very significant, and it is not evenly distributed.
The principal strain rate axes show that the shortening direction in the shear zone is NE-SW in the western part of the shear zone, while ENE-WSW in the eastern parts (Fig. 6). This means that the shortening direction turns subparallel with fault structures towards the E, not typical for strike-slip deformation zones. While the NE-SW shortening direction in the western part is roughly consistent with an ENE-WSW sinistral strike-slip zone, it is critical to observe that this area is characterized by almost pure shortening, with a very much subordinate extension component (Fig. 6). Active deformation is hence dominantly NE-SW to ENE-WSW shortening in the shear zone, which is curiously not clearly reflected in the neotectonic fault and fold patterns (Koroknai et al. 2020). A possible explanation for this could be long-wavelength, lithospheric-scale folding around NW-SE axes, which would not be observed on ordinary industrial reflection seismic lines (Dombrádi et al. 2010;Horváth and Cloetingh 1996). The possibly critical role of lithospheric-scale folding in inducing differential vertical motions during the inversion of the Pannonian basin has been demonstrated (Dombrádi et al. 2010), and we speculate that it could also explain the lack of structural geological evidence for the overall NE-SW shortening that characterises most of the Pannonian basin. In addition, we stress that neotectonic activity (i.e. activity in the last 6-8 Myrs) that is resolved by structural analyses such as Koroknai et al. (2020) or Fodor et al. (2005) may differ from the pattern of currently active tectonics due to the vast difference in time scale. In summary, while the presence of neotectonic sinistral faults has been demonstrated in the area of the shear zone (e.g. Koroknai et al. 2020), our results imply that overall the Mid-Hungarian mobile belt is currently not a prominent sinistral shear zone, but rather an area that accommodates decreasing amounts of shortening towards the NE, with a minor perpendicular extension component.

Derecske fault zone
One of the most prominent neotectonic features in Hungary is the Derecske fault zone in the East Pannonian basin (Fig. 5). It is originally a normal fault (linked to the formation of the so called Derecske trough), which was reactivated with sinistral kinematics (e.g. Koroknai et al. 2020). Our results support that the Derecske fault zone is currently characterized by shear-type deformation, i.e. the shortening and perpendicular extension components being almost equal in the area (Figs. 5d, 6). Considering the reported structure and kinematics of the Derecske fault zone, it is tempting to interpret this as a result of strike-slip displacements. However, the shortening direction (ENE-WSW) derived from the GNSS data is subparallel with the fault zone (Fig. 6), which is not straightforward to reconcile with a sinistral deformation zone, where shortening should be oriented NNE-SSW in an ideal case. One also needs to keep in mind that the interpretation of geodetic strain rates may be problematic for such a small area (a single, around 30 km long fault zone), since the smoothing of the velocity field might erase short-wavelength motions that could affect the local principal strain rate directions. However, the observation of parallel shortening and fault orientations also stand for the Mid-Hungarian mobile belt, which is a significantly larger area. Hence, we suggest that the question of what could accommodate the present-day, regional, faultparallel shortening in these regions might be worth addressing in the future.

Strain rate and seismicity
The analysis of the central Pannonian region (the area of the Mid-Hungarian mobile belt) has shown that geodetically derived strain rates cannot fully be explained by seismic deformation (Bus et al. 2009), supporting the view that aseismic creep might be an important deformation mechanism in the Pannonian basin (Gerner et al. 1999). The case for this hypothesis is further strengthened by the generally low seismic moment release in the regions of deep sub-basins with high heat flow (Lenkey et al. 2002). Such an area is the Hungarian side of the SW Pannonian basin, where the highest geodetic strain rates within Hungary (Fig. 5a) are coupled to relatively low seismicity (e.g. Bondár et al. 2018). The opposite is true for the North Hungarian Range, where very low to almost negligible geodetic strain rates (Fig. 5a) are coupled to relatively higher seismicity (Bondár et al. 2018), probably due to the difference in the rheological stratification of the lithosphere, defined by rock compositions and temperature (Lenkey et al. 2002).

Vertical velocities
Geological observations, such as eroded and tilted sediments, and river incision data clearly demonstrate that the margins of the Pannonian basin, such as the Transdanubian Range, the Western Carpathians, the East Carpathians, and the Apuseni Mountains have been uplifting during the entire Quaternary, while deep sub-basins such as the W-Pannonian basin and most of the E-SE Pannonian basin have experienced further subsidence Horváth and Cloetingh 1996;Ruszkiczay-Rüdiger et al. 2020;Ruszkiczay-Rüdiger et al. 2005;Sanders et al. 1999;Zuchiewicz 2009, for relations between references and geographical locations see Fig. 7). This pattern of vertical motions does not correlate well with the GNSS-derived vertical velocities of this study (Fig. 4b), or with that of Piña-Valdés et al. (2022), which show subsidence for almost the entire Circum-Pannonia region. However, it is clear that the margins of the Pannonian basin are still relatively uplifting with respect to the deep basins, which means a near 0 vertical motion of the margins compared to the 1-2 mm/ yr subsidence of the basins in the ITRF reference frame (Fig. 4b, Altamimi et al. 2016). Piña-Valdés et al. (2022) suggested that the wide-spread subsidence in Central-Eastern Europe is not related to horizontal strain rates, which would actually result in uplift for most of the Pannonian region, due to the ongoing shortening. Hence, the wide-spread subsidence is probably related to deep-seated, sub-lithospheric processes (Piña-Valdés et al. 2022). Given that the GNSS-derived relative vertical motions correspond well with the geological observations (Fig. 7), we suggest that a regional, larger wavelength source is  Sanders et al. (1999) responsible for the currently observed, wide-spread subsidence. We also emphasize that the current subsidence rates with respect to the ITRF reference frame are not compatible with clear geomorphological-geological evidences for Quaternary uplift of the basin margins (e.g. Ruszkiczay-Rüdiger et al. 2020).
These considerations point towards the superposition of smaller-wavelength (100 km) and larger-wavelength (1000 km) sources of vertical motions, with the large-wavelength component not being constant throughout the Quaternary. In agreement with previous studies, we interpret the smaller-wavelength source (causing the relative uplift of the basin margins) to be related to the compression of the thinned Pannonian lithosphere, causing lithospheric-scale bending (Dombrádi et al. 2010;Horváth and Cloetingh 1996). This process and vertical motion pattern might further be influenced by the compressed asthenosphere below the Pannonian basin, surrounded by a thick orogenic rim, resulting in the escape of the asthenosphere perpendicular to the generally NE-SW compression direction (Liptai et al. 2022). However, a quantitative effect on the topography of this deep-seated process is currently lacking.
Regarding the observed larger-wavelength subsidence (Fig. 4b), we speculate that it is related to the Europe-scale buoyancy-driven flexure and mantle flow induced by the melting of the ice in N-Europe since the last glacial maximum (~ 21 kyr ago). This process, Glacial Isostatic Adjustment (GIA), causes significant current day uplift rates due to isostatic compensation in formerly glaciated areas, such as in northern Europe. The same process also leads to a collapse of the former so-called peripheral bulge, a previously uplifted area that surrounded the ice sheet. During the collapse of the peripheral bulge the asthenospheric mantle flows towards the uplifting regions (Whitehouse 2018). Model studies predict that this collapse of the peripheral bulge may cover large parts of Europe (Marotta et al. 2004), leading to a rim of subsidence stretching from southern England, to the Mediterranean (Simon et al. 2018) and from the eastern Mediterranean towards eastern Europe, as shown by Vishwakarma et al. (2022), based on models from Caron et al. (2018). While these studies are partly based on vertical velocities derived from GNSS data, no data is used at the typical distances between the Pannonian basin and the extent of the former Fennoscandian ice sheet. Hence, a broad region of subsidence is based on model predictions of mantle flow, rather than a fit to data. This broad region of subsidence includes the Circum-Pannonian region (Marotta et al. 2004;Caron et al. 2018;Vishwakarma et al. 2022), and is estimated to be < 1 mm/yr. Hammond et al. (2021), Serpelloni et al. (2022) as well as Valdés et al. (2022) estimate a similar pattern and magnitude of subsidence from interpolated GNSS data, if we disregard a number of smaller regions of uplift that may be affected by mantle upwelling or isostatic adjustment. If this large-scale subsidence is indeed a response to deglaciation, it is a periodical phenomenon on the time scale below that of the Quaternary, meaning that relative vertical motions capture the average Quaternary displacements better than the currently observed absolute vertical motions. Adjusting the colormap of Fig. 4b so it becomes centred around the mean vertical displacement (− 0.9196 mm/yr) highlights relative vertical motions, which correspond quite well to Quaternary geological observation in Circum-Pannonian region (Fig. 7), such as the uplift of the Transdanubian range (Ruszkiczay-Rüdiger et al. 2020, the Western Carpathians (Jacko et al. 2021;Zuchiewicz 2009), or the Romanian mountain ranges (e.g. Matenco et al. 2007;Sanders et al. 1999).

Comparison to GNSS-based tectonic studies
A notable difference and improvement compared to previous GNSS studies focusing on the Circum-Pannonian region (Grenerczy et al. 2005) is demonstrating the slow WSW directed motion of the NE European foreland (in the north-eastern corner of Fig. 3), which contributes to the squeezing of the Pannonian basin. Lacking data in that region, the velocity interpolation of Grenerczy et al. (2005) showed opposite, slow NE-directed movement, which would decrease the shortening rates inside the Carpathian chain. Furthermore, the quantification of strain rate fields in our study is a major step forward since this last GNSSbased study focusing on the Pannonian region.
A more suitable base for comparison is provided by recent studies of Piña-Valdés et al. (2022) and Serpelloni et al. (2022). These studies incorporate GNSS data at a larger scale and while they include the Pannonian region their focus is on larger areas; the former is a Europe-wide study, while latter is focusing on the Mediterranean region. Because of different interpolation approaches, inverse distance smoothing and wavelet interpolation, respectively, some differences exist with respect to our results. Furthermore, in case of Serpelloni et al. (2022), where the Circum-Pannonian region is not densely covered, the strain rate field is qualitatively different, where their solution suggests dominant extension (positive dilatation rates) in the Northern Dinarides and large parts of the Pannonian basin (e.g. the area of the Mid-Hungarian mobile belt). This is the opposite of our results and those of Piña-Valdés et al. (2022), and also Kreemer et al. (2014), which are based on a denser data coverage for this region, and show a general sense of shortening, in line with published neotectonic observations and models (e.g. Bada et al. 2007;Koroknai et al. 2020).
Our approach of setting a threshold value for the minimum horizontal velocity uncertainty makes use of the phenomenon called nugget-effect in kriging literature (e.g. Montero et al. 2015), and that causes the interpolation to become increasingly dependent on observations farther away from the interpolation point. Kriging is a weighted average interpolation and by virtue of spatially averaging over a large number of observations it is possible to estimate velocity values that are at or below the level of the (assumed) uncertainties of the individual point velocities (at the expense of small spatial details). The benefit of using a minimum velocity uncertainty is that it relates the degree of smoothing to the level of observation errors and/or local small-scale effects of the GNSS-derived velocities. The effect of a minimum uncertainty is ultimately a smoothing of the strain rate field; we show the dependency of strain rate spatial smoothness on the minimum velocity uncertainty in Fig. 3. Locally, the amount of smoothing depends on the ratio of observation uncertainty (the sum of observation uncertainty variance and nugget variance) to observation variance (C(0)). This implies that the (high) northward velocities in the Dinarides are less smoothed compared to the northward velocities in the rest of the domain, and the eastward velocities, as the local observation variance for the north component in the Dinarides is higher than elsewhere, see the upper panels of Supplementary Fig. B1. Hence, the smoothness of the interpolated velocity field does not only depend on the spatial density of the network, such as in a wavelet approach as used in Serpelloni et al. (2022), but also on the local signal to noise of observations.
In case of Piña-Valdés et al. (2022), main differences compared to our study are the dilatation rate magnitudes at the North Hungarian Range and its continuation in Slovakia, as well as the SW Pannonian basin. In former case, our data and results show a consistently slow dilatational (< 2 nstrain/yr) deformation, while the results of Piña-Valdés et al. (2022) show a much faster deformation reaching ca. 5 nstrain/yr. In the SW Pannonian case, our data (Fig. 5c) show that the fast deformation in the Dinarides (up to 20 nstrain/yr) quickly fade away towards the Hungarian side of the SW Pannonian basin, reaching only 3-6 nstrain/yr (shortening). In contrast, according to Piña-Valdés et al. (2022), the Hungarian parts of the SW Pannonian basin is characterised by 8-10 nstrain/yr shortening, while similar values are also shown for the Dinarides. Given that the used GNSS datasets are very similar in the two studies, these differences in strain rate magnitudes are most likely rooted in the different methods for velocity field smoothing or interpolation. Our methods demonstrate that lower uncertainties of the data (i.e. when using the velocity uncertainties based on GNSS time series) result in higher deformation rates, however, these solutions become increasingly dominant in short-wavelength noise. This noise occurs spatial scales of the typical distance between GNSS-sites and can result from underestimation of GNSS-velocity errors. Especially in low-strain regions these errors, which are probably not related to tectonic motions (Fig. 2a), may become dominant (see also Araszkiewicz et al. 2016). In our stochastic kriging approach, we can mitigate these short-wavelength patterns by assuming a minimum uncertainty for all velocities. This effectively increases the interpolation weight of observations that lie behind the nearest neighbouring observation sites, leading to a smoother velocity field in low-strain regions. By increasing the velocity uncertainties, we interpret the short-wavelength velocity variations from site to site as resulting from local velocity variations that do not represent the velocity field at the inter-site scale. Returning to the dilatational patterns in the North Hungarian Range and the SW Pannonian basin of Piña-Valdés et al. (2022), we think that these may have resulted from amplification of velocity field noise. For this reason, we argue that higher uncertainty levels and the resulting smoother and lower strain rate magnitudes represent regional tectonic motions more realistically.

Conclusions
We investigated the active deformation of the Circum-Pannonian region through the interpolation of GNSS-derived velocity field and the derivation of the strain rate fields. We discussed the velocity and strain rate results in terms of tectonic processes behind the horizontal and vertical velocity fields and horizontal strain rate fields. Our main conclusions are the following: • To estimate a strain rate field that is virtually free of short-wavelength noise, we need to assume a minimum velocity standard deviation of 1 mm/yr. • The active deformation of the Circum-Pannonian region is largely dictated by the 2-3 mm/yr NNE-directed motion of the Dinarides, and the 0.5-1.5 mm/yr, WSW to SSW directed motion of the eastern areas (European foreland, East Carpathians, South Carpathians, Transylvanian basin). • The opposite-sense motions define a large-scale, on average NE-SW shortening and transpression-type deformation in the Dinarides and the Pannonian basin, while the East and South Carpathians undergo N-S extension. • Neotectonic structures generally show good agreement with the GNSS-derived strain rate field. However, pure shortening and fault-parallel shortening along major sinistral neotectonic structures in the Pannonian basin raise questions for the parity of neotectonic (i.e. activity in the last 6-8 Myrs) and currently active deformation processes. Also, aseismic creep and lithospheric-scale deformation processes, that are difficult to resolve with standard structural geological tools, might also be important parts of the present-day strain rate field. • The vertical velocity field shows 100 and 1000 km wavelength signals, the former probably related to the response of the Pannonian lithosphere-asthenosphere system to neotectonic basin inversion, while latter possibly explained by regional subsidence patterns