# Bayesian Estimation of Surface Strain Rates in the Peri-Adriatic, Balkans and Aegean region from GNSS Velocities

Marianne Métois, Cécile Lasserre, Aimine Méridi, Maxime Henriquet, Thomas Bodin

► To cite this version:

Marianne Métois, Cécile Lasserre, Aimine Méridi, Maxime Henriquet, Thomas Bodin. Bayesian Estimation of Surface Strain Rates in the Peri-Adriatic, Balkans and Aegean region from GNSS Velocities. Tektonika, 2025, 3.1, pp.82-98. 10.55575/tektonika2025.3.1.99. hal-05094138

HAL Id: hal-05094138

https://hal.science/hal-05094138v1

Submitted on 2 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Bayesian Estimation of Surface Strain Rates in the Peri-Adriatic, Balkans and Aegean region from GNSS Velocities

Marianne Métois ID*¹, Cécile Lasserre ID¹, Aimine Méridi ID¹, Maxime Henriquet ID², Thomas Bodin ID³

¹Université Claude Bernard Lyon 1, ENS de Lyon, Université Jean Monnet, CNRS, LGL-TPE, UMR5276 | ²Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Geoazur, Valbonne, France | ³Institute of Marine Sciences (ICM-CSIC), Barcelona, Spain

**Abstract** The calculation of spatial derivatives of the lithosphere surface velocity field has become a routine task in geosciences. It provides a continuous description of strain rates characteristics (dilation, rotation, strain rate intensity, etc) and facilitates the comparison between observed strain and theoretical strain computed from rigid or viscous lithospheric deformation models. However, this derivation is highly non-unique and standard methods may fail at correctly computing the strain rate tensor. This is especially the case in regions with strong spatial heterogeneities in strain rates like the Balkans, significantly limiting the understanding of their geodynamics. We tackle this issue in the peri-Adriatic and western Aegean region, using the bayesian transdimensional approach implemented in the B-Strain code, successfully applied in California. We show that our results, providing reliable probabilistic estimates of strain rates all over the region, are less prone to interpolation artifact than those produced from standard techniques, opening new perspectives for modeling and understanding the forces at play in this seismically active region of Europe.

Executive Editor:
Robin Lacassin
Associate Editor:
Jack Williams
Technical Editor:
Mohamed Gouiza

Reviewers:
Corné Kreemer
Kathryn Materna
Enrico Serpelloni

Submitted:
2 December 2024

Accepted:
14 April 2025

Published:
29 May 2025

## 1 Introduction

The tectonics of the Eastern Mediterranean is dominated by the convergence between the Eurasian, Nubian and Anatolian plates (Figure 1) and characterized by the highest seismic hazard in Europe. Some of the boundaries between the colliding plates are well defined, such as the Hellenic subduction zone (where Nubia/Anatolia converge at a rate of ~4 cm/yr; *Floyd et al.*, 2010; *Nocquet*, 2012) or the North Anatolian Fault (NAF, accommodating the Eurasia/Anatolia relative motion at a rate of 2.5 cm/yr; *Nocquet*, 2012), associated with clear seismic style (thrusting and dextral slip, respectively, Figure 1). In contrast, the junction between Nubia and Eurasia near the northwestern end of the Hellenic subduction trench remains poorly documented, while the connection zone between the western end of the NAF and normal faults in the Corinth Gulf appears broad with more distributed deformation (*Armijo et al.*, 1996; *Jolivet et al.*, 2013). GNSS measurements carried out in the last decades over the Italian Peninsula show that the Adriatic promontory has a slightly different motion from that of the Nubian plate. It is well described by the motion of two rigid blocks, Adria and Apulia (which boundaries remain unclear, Ap and Ad in Figure

1a), pushing as indenters into the Alps and the Dinarides with a velocity of 2 to 4 mm/yr towards the north-east (*D'Agostino et al.*, 2008; *Serpelloni et al.*, 2005; *Nocquet*, 2012).

The Balkan region located between the Alps and the highly straining Aegean domain has long been considered a stable part of the Eurasian plate (e.g. *Jolivet et al.*, 2013), experiencing negligible straining. As a result, most rigid block models proposed to explain GNSS surface velocities encompass the inner Balkans, north of the Rhodopes and Hellenides ranges (R and H in Figure 1a), within a single large and stable Eurasian block (*Reilinger et al.*, 2010; *Floyd et al.*, 2010; *Vernant et al.*, 2014; *Briole et al.*, 2021; *Ergintav et al.*, 2023). However, recent GNSS measurements published by *Metois et al.* (2015); *d'Agostino et al.* (2020); *Piña-Valdés et al.* (2022); *Serpelloni et al.* (2022) show that the entire peninsula is subject to significant deformation. Non-zero and spatially varying horizontal velocities can indeed be resolved resulting in a clockwise rotation towards the Aegean domain, extending up to central Serbia (see Figure 1b in an Eurasia-fixed reference frame). In this context, Albania appears as a NS-trending shear zone acting as the pivot region of this rotation, affecting the Dinarides to the north, the Albanides (both separated

<footnote>* marianne.metois@univ-lyon1.fr</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1** – Seismotectonic context of the studied Balkans and Aegean region. a - Tectonic context: active faults from the European Fault-Source Model EFSM20 project (Basili et al., 2022) are plotted as black curves, faults in bold are those discussed in this study. Rigid plates are Nubia (Nu), Eurasia (Eu), Thyrrenian (Ty), Apulia (Ap), Adria (Ad), Anatolia (An). Dotted and plain gray lines stand for major plate boundaries: North Anatolian Fault (NAF), Kefalonia Fault (KF), Eastern Alpine Front (EAF), Hellenic subduction backstop (HB), Calabrian subduction (CB). Other major structures are indicated: Corinth Gulf (cg), Scutari Peck lineament (sp), Pannonian Basin (pb), Moesian platform (mp), together with major mountain belts as the Alpine Arc (AA), Dinarides (D), Balkans range (BR), Rhodopes (R), Albanides (A), and Hellenides (H). Profile lines studied in detail in Figures 6 and 7 are plotted as dashed and dotted gray lines. b - Current day kinematics. White and purple arrows: cleaned combined GNSS velocity field from Piña-Valdés et al. (2022) in the Eurasia-fixed reference frame. c - Focal mechanisms and instrumental seismicity from Global CMT (Harvard) and USGS catalogues (Mw2.5+; USGS, 2025) over the 1997-2024 and 1990-2024 time spans, respectively. Current political boundaries are plotted as thin black curves.
by the NE-SW-trending Scutari-Peck lineament), and the Hellenides to the south (Figure 1a). The style, amount, and shape of the distributed deformation over the Dinarides and Albanides are still poorly constrained due to the lack of dense GNSS interseismic measurements and the low strain rates expected there.

Despite these limitations, several authors have proposed to map and analyze the strain rate tensor, $\dot{\epsilon}$, over this broad area based on combinations of GNSS velocity fields measured at a set of discrete locations (e.g. Metois et al., 2015; d'Agostino et al., 2020; Briole et al., 2021; Piña-Valdés et al., 2022; Serpelloni et al., 2022). The goal is to decipher the balance of forces acting on the lithosphere and to better estimate the seismic hazard in the area. The inversion of $\dot{\epsilon}$ from discrete measurements of horizontal velocities ($V_{east}$ and $V_{north}$) is known to be highly non unique: most commonly used inversion techniques are prone to biases, and generate artifacts (Baxter et al., 2011). The solution often depends on user-defined parameters, which prevents uncertainty assessment, thus hindering further interpretation, in particular where strain rates are low (Pagani et al., 2021).

In this study, we invert for the strain rate tensor over Italy, the Balkans and continental Greece. We use (i) the recently published combined GNSS velocity field by Piña-Valdés et al. (2022), that offers the best coverage to date, in particular in Croatia and Bosnia-Herzegovina, and (ii) the B-Strain code published by Pagani et al. (2021), that allows a bayesian inversion of the strain rate tensor using a transdimensional approach. We provide probabilistic continuous maps of the strain rate tensor invariants (e.g. second invariant, dilatation), vorticity, and interpolated horizontal velocities. We examine the statistical indicators of their probability density functions (PDF), and make the results widely accessible to all on an online webtool and an archived repository.

These results open the door to detailed tectonic and geodynamic analysis of the area based on a refined knowledge of surface deformation. To highlight their valuability for further geodynamical interpretations, we discuss some key characteristics of our surface strain rate maps along representative profiles.

# 2 Data and Method

## 2.1 GNSS horizontal interseismic velocity field

The interseismic velocity fields published so far in the western and central Balkans still suffer from the lack of publicly available raw GNSS data, despite the existence of private permanent GNSS networks with geodetic standards (e.g. Metois et al., 2015; d'Agostino et al., 2020; Piña-Valdés et al., 2022). Here, we chose to use the horizontal GNSS velocity field published by Piña-Valdés et al. (2022) in our area of interest, ranging from the western Alpine front to the Black sea and from the Hellenic backstop to northern Carpathians (Figure 1-b). This velocity field comes from a combination of

several independent velocity fields, rotated in a common Eurasia-fixed reference frame. It has been cleaned from low quality data (see Piña-Valdés et al., 2022, for further details on their preferred filtered and smoothed velocity field, that we use here). 1867 horizontal velocities computed over variable time-spans over the 1996-2020 period are encompassed in our study area, with associated rescaled uncertainties, of 0.2 mm/yr on average on the East and North components. For comparative tests, we also use the raw unfiltered data set (including the outliers) combined by Piña-Valdés et al. (2022); their corresponding analysis is provided in the Supporting Information.

Since vertical velocities are often associated with large uncertainties due to several processing issues and non-tectonic deformation, and since it is still unclear how to properly incorporate these velocities in pseudo-3D strain rate tensor calculation (see for instance Mazzotti et al., 2005; Piña-Valdés et al., 2022), we discard here the vertical data set provided by Piña-Valdés et al. (2022).

## 2.2 Strain rate calculation using B-Strain

Our study area is the most seismically active region in Europe (Figure 1), but the seismic style and expected strain rates are spatially very heterogeneous: the relatively slowly straining Balkans peninsula contrasts with surrounding regions with major active structures, slipping at rates ranging from several millimeters per year (in the Apennines and Corinth Gulf) to a few centimeters per year (along the NAF and the Hellenic trench). Because of these large variations in strain rates, standard interpolation methods of measured GNSS velocity fields would tend to fail to correctly recover both the spatial pattern and the amplitude of the expected strain tensor components. A well-known drawback is that the level of smoothing for interpolating velocities measured at stations needs to be defined by the user. The amplitude of the spatial derivatives (i.e. the components of the estimated strain rates) therefore depends directly on this chosen level of smoothing (see Pagani et al., 2021, and the discussion therein). If the smoothing is chosen to be large, regions of high strain rates will be overall correctly recovered and no signal will be visible in low deformation areas; whereas if the smoothing is low, the results will be highly sensitive to outliers and will appear artificially patchy. Some refined interpolation and inversion methods, such as VISR developed by Shen et al. (2015) and discussed further in this study, invert for a spatially varying smoothing parameter to get around this problem. Wavelets decompositions (Tape et al., 2009; Serpelloni et al., 2022) are other approaches to manage varying smoothing but they remain highly dependent on the network density and regularization choices. Another issue is that such standard interpolation techniques tend to generate interpolation artifacts when crossing areas with high velocity gradients, depending on the network geometry (Baxter et al., 2011).

Here, we use the B-Strain code presented in Pagani et al. (2021) and successfully tested in California
to address these methodological pitfalls and generate probabilistic maps of the strain rate tensor over our study area. B-Strain is made publicly available via a Gitlab repository and associated documentation. With B-Strain, we estimate four types of continuous fields from our discrete dataset of surface velocities measured at GNSS stations, namely:

* the horizontal velocities $v_{east}$ and $v_{north}$

* the vorticity $\omega = \frac{1}{2}(\partial_y v_x - \partial_x v_y)$, with positive values for counterclockwise rotation,

* the divergence $d$ defined as $d = \dot{e}_{xx} + \dot{e}_{yy}$ with $\dot{e}_{ij} = \frac{1}{2}(\partial_j v_i + \partial_i v_j)$, here positive when dilation dominates while negative values represent compression,

* the second invariant of $\dot{e}$ defined as $I_2 = \sqrt{\dot{e}_{xx}^2 + \dot{e}_{yy}^2 + 2\dot{e}_{xy}^2}$, which is a convenient way to estimate the strain rate intensity.

We sample the probability density functions (PDF) associated with each component of the velocity field, strain rate tensor, and its invariants, by a reversible-jump Monte-Carlo Markov chain algorithm. The continuous horizontal velocity field to invert for is parameterized with a Delaunay triangulation with a variable number of nodes that do not have to be at the observation locations (Pagani et al., 2021). In this way, the complexity of the solution is not imposed by the user but retrieved by the inversion. B-Strain generates a large ensemble of 2D models which distribution can be used to approximate the probabilistic solution. Since each individual model has a different parametrization, the statistics on the ensemble provide continuous maps. We also jointly invert for a hyperparameter that scales all data uncertainties, which is equivalent to estimating the level of data fitting required by the data. As in California where we used the MIDAS data set to invert for strain rates (Pagani et al., 2021), we parametrize B-Strain so that the formal uncertainties given by Piña-Valdés et al. (2022) are allowed to be increased by up to 10%. This a priori bound reflects our degree of confidence in the velocity field, in average.

To ensure good convergence and sampling of the model space, we run the inversion on 96 parallel processors for 1 million iterations (a 2D Delaunay model is collected at each iteration), resulting in a full day of calculation on the PSMN mesocenter facilities in Lyon. Figure 2-a shows that the inversion converges toward models containing around 380 triangles, and that the inversion process requires the uncertainties associated with horizontal velocities to be 10% higher than the initial values (thus reaching the maximum allowed value that we impose on the inverted hyperparameter). We verify that the nodes location is mostly controlled by the level of information contained in the dataset (here the velocity gradients) rather than by the stations density. Indeed, the highest node density shown in Figure 2-b is not correlated with the most densely instrumented areas (the Po plain and the Appenines in northern and central Italy, for instance). We run B-Strain on a 200×200 pixels grid, after projecting the GNSS stations coordinates

onto a cartesian grid in the UTM 34N, resulting in a pixel size of about 10 km. The full description of the parametrization used for this computation, including a priori bounds on each parameter, is available in the archive repository hosting our results (file parameters.in, see Data Availability section).

# 3 Results

## 3.1 Visualisation of the results

Because the outputs from the B-Strain inversion are Probability Density Functions (PDFs) associated with the invariants or principal directions of the strain rate tensor (divergence, second invariant), the vorticity and velocity components, at each grid point of the inverted fields, their visualization is challenging. Our results can be explored online, dynamically, using the dedicated web tool B-Strain Plotter and associated documentation. Maps of the statistical indicators associated with the PDFs (e.g. mean, median, maximum mode) can be downloaded as geotiffs, and the full PDFs are available in the associated archived repository. The PDFs can be represented as full distributions for each single pixel, or along profile lines. We developed this web tool to make our results reusable and allow their exploration and exploitation by people interested in studying strain rates over a specific area or profile. In this study, we present only a limited set of figures that will hopefully help the reader to get a general idea of the main findings associated with our computation. However, many other insights can be derived from the B-Strain outputs and will serve as a basis for future work.

We plot in Figure S1 (in Supporting Information) the residuals between the interpolated velocity field (median) and the observations, and in Figure S2 we provide full map views of the PDF’s median value for the east and north components. In Figure 3, we show maps of three different statistical modes of the PDF (median, average, and maximum mode of probability) for the divergence, second invariant and vorticity of the velocity field (see definitions in Section 2). To facilitate their interpretation, we mark the contour line where the standard deviation of the PDF associated with the horizontal velocity is higher than 3 mm/yr (green line on Figure 3), our chosen threshold to mask areas (in all inverted fields) with poorly constrained results (see also Figures 2-c and 4-right). We could have chosen to mask each inverted field independently based on the standard deviations of their own PDF. However, we prefer using a single masking strategy for all the quantities we plot, for clarity. It is worth noting that large standard deviations in horizontal velocities are observed in the Marmara sea in northwestern Turkey and Aegean sea in Greece, where GNSS observations are sparse, but also where expected gradients and velocities are large. The chosen mask may therefore be too restrictive there. Overall, maps of median, average, or maximum mode show consistent patterns in well-constrained areas, regardless of the quantity considered. However, maximum mode maps are systematically rougher, with triangular shapes due to our parametrization remaining visible in some
<figure><img src="image_4.png" /></figure>

**Figure 2** – a- Top and center: Evolution of the number of Delaunay triangles and of the scaling parameter of measurement uncertainties during the iterations. The so-called burnin period corresponds to the first set of iterations (here 400k) that are discarded in the final results. Bottom: histogram of the number of models as a function of the number of triangles used. The x-axes boundaries are the a priori bounds used as B-Strain input. b- Map of the density of nodes (i.e. triangle submits) for each pixel at the end of the convergence. The density is the highest where expected gradients of deformation are the highest (NAF, Apennines, Corinth Gulf), and not where the GNSS station density is the highest. c- Map of the standard deviation (in mm/yr) computed for the horizontal velocity based on the standard deviation of the PDF obtained for interpolated $V_{East}$ and $V_{North}$, with contours every 3 mm/yr. In the following, areas above this ~3 mm/yr value, either offshore or where data are lacking, are masked as unresolved areas.
<figure><img src="image_5.png" /><figcaption>**Figure 3** – Maps of the divergence (top, positive stands for extension), second invariant (middle) and vorticity (bottom, positive in counter-clockwise) of the interseismic velocity field in the study area, extracted from the B-Strain outputs. Three statistical indicators of the PDFs are shown, namely the median (left), average (central) and maximum probability mode (right). The green line represents the 3 mm/yr standard deviation value we use as a quality threshold (see corresponding low quality / masked areas on Figures 2 and 4).</figcaption></figure>

places (see Discussion section). This drawback has also been observed in California and recall that maps of statistical indicators are not equivalent to consistent physical models. These remnants of the Delaunay triangulation tend to fade out in median and average maps, making them more suitable for geophysical interpretation. However, we consider median maps to be more reliable than average maps, as the average of an asymmetric distribution may be significantly different from the median and maximum mode values. This is strikingly the case in the Adriatic sea (an area considered to be poorly constrained, i.e. with horizontal velocities above the quality threshold), where the average values of the second invariant are unrealistically high compared to the median values (see Figure 3 middle center). In this area, some models presenting unrealistically high values of the second invariant were retained during the inversion process, due to the lack of strong constraints on the velocity values, thus biasing the average value. These discrepancies between average and median maps arise in all areas masked based on our quality criteria.

We therefore recommend not to interpret the strain rates characteristics in these areas.

Note that the offshore areas with poor data coverage depict the highest strain rates. This is because we parametrize and invert for velocities rather than their spatial derivatives (i.e. strain). In the regions where no data is available, the posterior probability distribution on velocities is similar to the prior. The Delaunay nodes then take completely random velocity values between $V_{min}$ and $V_{max}$, without any spatial consistency, resulting in a white noise on velocity and producing extreme spatial derivatives, and thus high strain levels.

## 3.2 Strain rate maps analysis and comparison with previous work

We compare our maps of the PDFs median to the ones published by Piña-Valdés *et al.* (2022) using the exact same dataset but a different method for strain rate inversion, the VISR3D method, one of the most commonly used method based on an adaptative
<figure><img src="image_6.png" /><figcaption>Figure 4 – Maps of the divergence (top), second invariant (middle) and vorticity (bottom) of the velocity field in the study area, obtained using the VISR3D interpolation method (left; Piña-Valdés et al., 2022) or the B-Strain inversion method (center, median of the PDF from this study). Right column shows the shaded topography, active faults from SHARE project and the B-Strain outputs to allow a better location of major features. Dashed black lines: -2 nstrain/yr contour for divergence, 5 nstrain/yr for I₂, zero contour for rotational. Arrows point toward significant features discussed in the main text.</figcaption></figure>

smoothing strategy to obtain an interpolated velocity field, fully described in Shen et al. (2015). The algorithm has been modified by Piña-Valdés et al. (2022) to calculate a pseudo-3D second invariant, assuming a fully elastic lithosphere, resulting in a slightly different formula for $I_2^{3D}$:

$$I_2^{3D} = \sqrt{\dot{\epsilon}_{xx}^2 + \dot{\epsilon}_{yy}^2 + 2\dot{\epsilon}_{xy}^2 + \left(\frac{-2\mu_0(\epsilon_{xx} + \epsilon_{yy})}{\lambda + 2\mu_0}\right)^2}$$

with $\lambda$ and $\mu_0$ the Lame coefficients for the Earth's crust. Part of the differences observed between our $I_2$ estimates and $I_2^{3D}$ obtained by Piña-Valdés et al. (2022) could therefore be associated with this pseudo-3D assumption. These differences remain negligible at first order and for the $\lambda$ and $\mu_0$ values characteristic of the crust.

Indeed, results of the two methods are overall in good agreement, with a correlation close to 1 for the divergence, second invariant and vorticity (see Figures S3 and S4 in Supporting Information).

Both methods show the following consistent patterns, that also arise from other calculations made on slightly

different GNSS velocity fields (Metois et al., 2015; Serpelloni et al., 2022, see also Figure S5 in Supporting Information):

* Dilatation (positive values of divergence) dominates in western Turkey, in the Apennines down to southern Italy and northeastern Sicily, Greece, Macedonia and the inner Albanides, while compression dominates in the Dinarides (and possibly in the Adriatic although poorly constrained), the Po plain and the eastern Alpine front. A sharp transition between dilational and compressive domains (dashed contour in the right panel of Figure 4) can be followed along the Northern Apennines front (Label 1 on central panel of Figure 4) and along the coast of Albania (Label 2 on central panel of Figure 4);

* The second invariant is maximal (>100 nstrain/yr) along the North Anatolian fault in Turkey and the Corinth Gulf in Greece, where a clear decrease in both $I_2$ amplitude and width of the deformation zone from the western to the eastern part of the gulf is observed. High values are also
observed in northeastern Sicily and along the Apennines range with two localized maxima in the Northern (Norcia-L'Aquila) and Southern (Irpinia) Apennines;

* A clockwise rotation is recovered over Greece and western Macedonia, along the NAF, as well as over the so-called Apulia microplate encompassing Albania, Montenegro, the italian Puglia region and south of the Adriatic platform (see Figure 1 for microplates location). A counterclockwise rotation is observed in the Po plain, in the western Alps and in the northeastern Adriatic, resembling the previously described Adria microplate (D'Agostino et al., 2008; Serpelloni et al., 2005). A counterclockwise rotation is also observed in southwestern Anatolia and as patches in central and southwestern Italy.

However, significant differences emerge in the results of the two approaches. The B-Strain maps appear much less patchy than those derived from VISR3D (and from the wavelet-based method used by Serpelloni et al. (2022) in the area, see Figure S5 in Supporting Information), especially in areas of low strain rates (e.g. the inner Balkans and the Alpine arc), highlighting structures that were previously difficult to detect. In particular, the NS-trending limit between the compressive and extensive domains in the Balkans, although seen on both the VISR3D and B-Strain solutions as mentioned above, can be drawn with much more confidence (as outlined by the dashed contour in Figure 4 right panel and label 2 on its central panel). In the continuity to the South, offshore the western coast of Greece, B-Strain also recovers, despite the low confidence, the limit between compression and extension associated with the Hellenic trench along the entire Hellenic arc, unlike VISR3D (label 4 in Figure 4). Such a deformation contrast is expected offshore but is usually imposed in strain rates maps by using information coming from focal mechanisms or by locally imposing a specific rheology of the lithosphere (Perouse et al., 2012; Metois et al., 2015; Kreemer and Chamot-Rooke, 2004; Kreemer et al., 2004; d'Agostino et al., 2020). In the north of our study area, the northern limit of the compressive domain in the eastern Alps and the Dinarides also appears clearly both in the divergence and second invariant maps, as an arc-shaped limit (label 3 in central panel of Figure 4). North of this limit (~47°N), strain rates are uniformly low with divergence and second invariant values lower than 10 nstrain/yr. Other interesting features arise clearly in the B-Strain results from the analysis of the vorticity map. In particular we observe an abrupt change in the direction of the rotation along a ~N-S-trending structure crossing Eastern Macedonia and Eastern Greece near longitude 23°E, roughly following the Vardar suture zone (label 5 in Figure, 4 Jolivet et al., 2013). The eastern Alpine compressive front and the NW-SE, mainly right-lateral, Dinaric fault system in the continuity (Moulin et al., 2014), are also associated with an abrupt change in the sense of rotation (label 6 in Figure 4).

To further compare the results of the B-Strain code (this study) and the VISR3D method (Piña-Valdés et al., 2022), we draw several profile lines crossing major active structures (see their location in Figures 1-a and 7) presented in Figures 6 (for roughly N-S profiles), and 7 (E-W profiles). The variation of the East and North components of the interpolated velocity along these profiles is also shown on Figures S6 and S7 Supporting Information. These profiles present the full PDF, normalized by its area, and the median, average and maximum modes (Figure S8, in Supporting Information, presents an alternative view of these profiles normalized based on the PDF maximum value). Where strain rates are high, for instance across the North Anatolian Fault (NAF, profile 1 on Figure 6 left), or the Corinth Gulf (cg, profile 2 on Figure 6 center), B-Strain recovers higher amplitudes than VISR3D, with more than 100 nstrain/yr difference in the second invariant locally (see Figure S4, in Supporting Information, for map view of the residuals between the two approaches). In areas of low strain rates, the B-Strain-derived profiles are smoother with lower amplitudes in average: a large dispersion is observed in the divergence, second invariant or vorticity values obtained with VISR3D where B-Strain gives nearly zeros (see spreaded VISR3D values at zero B-Strain values in Figure S3 in Supporting Information). For instance, going from Bucarest (44.33°N,26.03°E) to Izmir (38.41°N,27.13°E) along profile 1 (see Figures 6 and 5), the median of the second invariant given by B-Strain increases from nearly zero to 20 nstrain/yr over a 150 km distance, from the Balkans range to the Thrace basin, then jumps abruptly to 220 nstrain/yr where crossing the NAF and the Marmara sea, before decreasing abruptly back to 50 nstrain/yr and remaining stable up to the Bergama graben (BG), where it increases again and reaches ~75 nstrain/yr. Along the same profile 1 line, the VISR3D method retrieves a smooth increase of the second invariant, centered on the NAF, but reaching 120 nstrain/yr only, while missing the increase in strain rate associated with the Anatolian grabens (EFZ, BG) - also compare in map view Figure 4. Similarly, the B-Strain code correctly recovers the strong variation in the divergence expected across the NAF left-lateral fault, contrary to the VISR3D retrieving a smoother increase. The two methods differ significantly in the amplitude (by nearly 100 nstrain/yr) and width (from simple to double) of the negative vorticity anomaly observed across the NAF. Overall, the results from the VISR3D method provide values that are often reasonable given the local PDF, i.e. are amongst the most probable ones. However, in some places, the values are out from the 90% probability envelope, and are therefore unlikely values.

## 3.3 Principal directions and tectonic style

We derive the principal components of the strain rate tensor $\epsilon_1$ and $\epsilon_2$ for each pixel, and represent the mode of maximal probability of their PDF in Figure 5-left on a decimated grid, since the visualization of the full PDF is particularly difficult for these quantities
<figure><img src="image_7.png" /></figure>

**Figure 5** – Left: active faults from the EFMS project (Basili et al., 2022) together with principal directions of strain rate plotted on top of the median divergence. Dashed green line stands for the masking contour based on the 3 mm/yr standard deviation of the horizontal velocity. Right: tectonic style categorized based on the $s = \frac{3}{4} + \frac{1}{\pi}tan^{-1}(\epsilon_2/\epsilon_1)$ criteria. The active faults are color-coded based on their assumed principal motion based on Basili et al. (2022) (red: reverse, green: strike slip, blue: normal). Profile lines 1 to 5 are plotted for reference (see Figures 6 and 7).

(Pagani et al., 2021). Windroses of the principal components' PDF can be vizualised for each pixel with the B-Strain Plotter webtool. We propose in Figure 5-right to plot the dominant tectonic style as expected from GNSS velocities based on the commonly used ad-hoc criteria $s = \frac{3}{4} + \frac{1}{\pi}tan^{-1}(\epsilon_2/\epsilon_1)$ ($s > 0.5$ stands for extension, $s = 0.5$ for pure strike-slip, $s < 0.5$ for compression, see Houseman and England (1986); Gordon and Houseman (2015)). These representations allow for an easier comparison with the tectonic style given by focal mechanisms and tectonic studies in the area.

Based on this criteria, most of the Italian peninsula appears to experience normal faulting with a ~N45 extension direction, roughly perpendicular to the Apennines belt. In the large extensive domain covering Greece and Macedonia, extension is dominantly N-S, with localized rotation of this maximum direction of extension in particular in the Albanides where it strikes N100, and in northwestern Greece around Ioaninna (20.9°E, 39.6°N) where it strikes N120. Compression is ~N15 to N20 across the northern Apennines front and Dinarides, while the principal compressive direction strikes roughly NS across the Eastern Alpine front. A restricted compressive zone is observed above the central Hellenic subduction zone in Crete, with a N20 dominant direction, in agreement with previous findings of Serpelloni et al. (2022) (see also Figure S5 in Supporting Information). Narrow areas are outlined as prone to dominant strike-slip motion, in particular along the NAF up to the Evia gulf in its western continuation (~24°E, 40°N) as well as in the northern Balkans (southern Serbia, Bosnia Herzegovina, and western Croatia). Finally, the Western Alps, where strain rates are very low and therefore more uncertain, are characterized by larges zones of extension and strike-slip with limited patches of compression, in

agreement with previous studies (Walpersdorf et al., 2018; Mathey et al., 2021; Grosset et al., 2023).

# 4 Discussion

## 4.1 Benefits and limitations of the B-Strain method

The application of the B-Strain inversion method to the GNSS velocity field of Piña-Valdés et al. (2022) over the Balkans confirms its efficiency for computing velocity gradients over regions of highly contrasted strain rates. The maps of the median of standard strain rates indicators presented in Figures 3, 4 and 5 show little variations in low strain stable areas, together with sharp and significant jumps where crossing major tectonic features (NAF, Corinth Gulf, etc), as expected.

The method therefore correctly retrieves smooth variations and abrupt gradients using a single parametrization. In particular, we get rid of patchy oscillations, often present in previous strain maps in the inner Balkans area (Metois et al., 2015; Piña-Valdés et al., 2022) but associated with interpolation issues that prevent the interpretation of small gradients (Baxter et al., 2011). We are therefore confident that we can use our results to go further in the tectonic interpretation of strain rates, and encourage the community to use the B-Strain Plotter tool to do so.

However, the transdimensional inversion method applied in the B-Strain code presents some limitations that must be underlined. In particular, we would like to recall that, even if we claim getting rid of the subjective choice of smoothing parameters, the a posteriori probabilistic solution depends on both the parametrization and on the a priori distribution, as developed below. Furthermore, for visualization
<figure><img src="image_8.png" /><figcaption>**Figure 6** – Normalized PDFs associated with divergence (top), second invariant (center) and vorticity (bottom) along three ~NS profile lines numbered in Figure 5 from 1 to 3. From left to right: from Budapest to Izmir crossing the North Anatolian Fault (NAF, left column), from Belgrade to southern Peloponnese crossing the Corinth Gulf (central), from Verone to Calabria along the Apennines belt (right). Y axis is in $10^{-7}$ or $10^{-8}$ strain/yr, as indicated in the upper left corner. Normalized PDF is color coded, while median, average and maximum mode evolution are represented as solid lines. Dashed blue lines stand for the same quantities extracted from the VISR3D results from Piña-Valdés *et al.* (2022). The top panel shows topography and major geographic or tectonic features. Gray shaded areas stand for badly resolved zones based on our quality criteria (standard deviation of the horizontal velocity PDF is higher than 3 mm/yr). Some of the main tectonic features crossed by the profiles are indicated as vertical dotted grey lines (Basili *et al.*, 2022; Ganas *et al.*, 2018).</figcaption></figure>

purposes, our ensemble of sampled models representing the posterior solution is depicted with histograms whose bounds and number of bins need to be chosen by the user. For example, in our study area, we expect the second invariant to range from 0 to 1000 nstrain/yr (expected values in stable Eurasian plate and highly straining NAF, respectively), and we only plot the posterior distribution within these limits. Of course, increasing the limits and the number of bins leads to larger files to store as outputs. Similarly, the spatial resolution of our outputs could be increased and may provide a more detailed view

of the strain rates' spatial variations, but it would result in large files.

In the transdimensional inversion method developed by Bodin and Sambridge (2009); Sambridge *et al.* (2013), the level of noise in the data can be treated as an unknown parameter to be inverted for. In the B-Strain code, this is done through the inversion of a hyperparameter that scales the GNSS velocities uncertainties and that is allowed to vary from zero (when uncertainties are much larger than the real errors) to
<figure><img src="image_9.png" /><figcaption>**Figure 7** – Same as for Figure 6 but for roughly EW profiles (profiles 4, 5 and 6 in Figure 5) ranging from Grosseto to Budapest crossing the Apennines and the Dinarides (left), and from Bari to the Black sea crossing the Albanides and Rhodopes ranges (right). Profile 6 (Figures 1 and 5) crosses Eolian islands and Sicily from North to South. Some of the main tectonic features crossed by the profiles are indicated as vertical dotted gray lines.</figcaption></figure>

1.1 (when uncertainties have to be increased by 10% to better resemble the expected level of noise; this appears to be the case in our study, as shown in Figure 2). However, the observed velocity field inverted here comes from a refined combination of various heterogeneous data sets (see Piña-Valdés et al., 2022, for full details) which have their own processing strategy and method to estimate velocity uncertainties. It is therefore rather unlikely that a single scaling parameter can capture alone the true level of noise and uncertainties associated to the measurements. This, amongst many other reasons, is in favor of making all raw data from GNSS networks available to ensure a self consistent processing and correct estimates of uncertainties. Many networks, in central Balkans in particular, do not freely distribute

these raw data to date.

Finally, to emphasize the resistance of our method to outliers, we run B-Strain with the exact same parametrization on the raw velocity field provided by Piña-Valdés et al. (2022), i.e. with outliers. The results provided in Figure S9 (in Supporting Information) are very consistent at the first order and for median values with the ones obtained with the cleaned velocity field. Local discrepancies arise in the amplitude and shape of highly straining zones, in particular in northern Apennines and across the Corinth Gulf when looking at divergence and vorticity. The second invariant appears much less sensitive to outliers.
## 4.2 Current strain rates and active tectonics: validation and new findings

### 4.2.1 GNSS versus seismotectonic observations

Our strain rates computation allows for a detailed confrontation of the current surface strain rate patterns observed by GNSS with the active tectonics knowledge of the region, derived from instrumental and historical seismicity and neotectonic studies.

Where significant earthquakes have been recorded since 1990 and associated with a focal mechanism, the tectonic style predicted by GNSS based on the $s = \frac{3}{4} + \frac{1}{\pi}tan^{-1}(\epsilon_2/\epsilon_1)$ criteria is consistent at the first order with the one observed through seismicity. The Figure S10, in Supporting Information, shows that earthquakes occurring in areas where $s$ is lower than 0.7 are mostly reverse, mostly normal where $s$ is higher than 1.3, and strike-slip for intermediate values of $s$. However, all seismic styles are observed whatever the value of this $s$ criteria, which has therefore to be seen as a first order indicator of the current seismic style. This can partly be explained by the structural inheritance of the lithosphere, i.e. the fact that the occurrence of a rupture depends jointly on (i) the stress state resulting from the residual crustal stress accumulated over the successive seismic cycles and lithosphere history and the current-day stress; and (ii) the orientation of pre-existing faults.

The EFSM data-base of active faults (European Fault-Source Model; Basili et al., 2022) assigns a dominant tectonic style to each documented fault, that we also compare with our ad-hoc criteria $s$ on Figure 5. The agreement is overall very good, with some notable exceptions. For instance, the complex transpressive zone in the slovenian and croatian Dinarides (~39-42°N) is categorized as a nearly purely strike-slip area based on our simplified categorization, while thrust faults exist. Similarly, although a thrust style is expected based on our strain rates in the Pannonian basin, strike-slip faults are mapped. Of course, in Romania, the deep thrust-type seismicity associated with the old slab break-off under Carpathian is not reflected by the surface strain style, dominantly extensive.

We analyze the relationship between strain rates and known active faults along six profile lines shown in Figures 6 and 7. Based on the elastic loading theory, we expect the second invariant $I_2$ to be maximal in the vicinity of active locked faults, at distances from the fault trace that depend on the geometry and indepth extent of the seismogenically locked zone (Okada, 1985). The norm of the divergence is supposed to be maximal near normal and thrust faults, and at the extremities of locked strike-slip finite faults or segments. An abrupt change in the divergence's sign should also be observed when crossing a major strike slip fault (see Pagani et al., 2021, for instance). Vorticity may vary from one side of the fault to the other if the fault accommodates high velocity gradients. Amongst the numerous active faults described over our study area, some are associated with

well resolved discontinuities in the strain rate and we list in the following some of the main interesting features contained in our results, either refining previous studies or bringing new information.

### 4.2.2 Strain rates highlight along-strike segmentation of the Apennines

Strain rates are far from uniform along and across the Apennines mountain belt as previously discussed in the literature (e.g. D'Agostino, 2014; Serpelloni et al., 2022; Nucci et al., 2023). Profile 3 (Figure 6 right panel) helps better quantifying the along-strike variations of invariants from the Po Plain to Calabria, with a 22 nstrain/yr jump in $I_2$ distributed among the Ferrara thrust system and the other North Apennines thrusts. Significant variations in divergence and vorticity are observed in between these systems, with dominant compression and counter clockwise rotation in the Po Plain, which appears consistent with the expected deformation of an accretional prism. The belt itself is dilating down to Potenza (40°N) and experiencing very low rotation rates, consistently with the dominant normal-faulting earthquakes striking Italy, with an overall decrease in $I_2$ and dilatation rates toward South. A clear decrease in $I_2$ and dilatation rates is observed in between the Northern and Southern Apennines where crossing the Tremiti and Mattinata NW-SE strike slip faults, and between the Southern Apennines and Calabria where crossing the Potenza and Pollino faults, confirming a segmentation of the belt in terms of current-day straining. The Pollino fault is also marked by a sharp jump in vorticity, suggesting that this fault zone is a major active feature.

We further explore the relationship between strain rates and active normal faults in the Apennines by looking at an EW profile crossing the belt in Northern Apennines, at the latitude of Grosseto and Budapest (47°N, profile 4 in Figure 7). The location of two major normal faults systems crossed by the profile (the Alto Tiberina and Umbria faults) correlates with the highest values of $I_2$ and divergence (~35 nstrain/yr). A sharp decrease in both quantities is observed just East of the Umbria fault, 20 km west of the main Apennine thrust's front.

### 4.2.3 Contrasted tectonic styles between the Nubian and Tyrrhenian domains in Sicily

Our strain rate maps also highlight key features affecting the island of Sicily in southern Italy, a strongly deforming region, as attested by the occurrence of major historical earthquakes along the eastern sicilian coast (e.g. Noto earthquake in 1693, Mw ~ 7.4, Messina earthquake in 1908, Mw ~ 7.1) and Calabria (e.g. Calabrian earthquakes in 1783, Mw 5 to 7). Thanks to GNSS measurements on the Eolian islands, ~ 50 km north from the sicilian northern coastline, the B-Strain derived deformation field is reliable and can be interpreted not only onshore but also offshore Sicily, in the southeastern Tyrrhenian Sea. Note that the activity of Mount Etna volcano is associated with high strain rates ($I_2$ ~ 80
nstrain/yr) in eastern Sicily that are not directly related with the regional tectonic setting.

The divergence velocity field (Figure 4 and profile 6 on Figure 7) shows clear compression (50 nstrain/yr) at ~38.6°N along the E-W trending thrust fault systems of the Ustica and Eolian islands, accommodating the slow convergence between Nubia and the Tyrrhenian Sea (Figure 1; Ventura et al., 2014; Palano et al., 2012). Southward in onshore NE Sicily, the extension dominates and reaches up to 60 nstrain/yr within the Nebrodi range and across the Messina strait, together with a net clockwise rotation that may be associated with the SE retreat of the Ionian slab (Figure 1; Ventura et al., 2014). We also note that the maximum $I_2$ (70 nstrain/yr) in that region is located in the vicinity of the so-called Cefalù-Etna seismic zone consistently with recent InSAR-based velocity fields (Henriquet et al., 2022a). Significant extensional strain continues toward the east across the Messina strait and within Calabria up to the southern Apennines. At ~37.5°N in Central Sicily, subtle compression around 10 nstrain/yr is evidenced when crossing the front of the Sicilian chain (Gela Thrust Front, profile 6 on Figure 7), consistently with previous estimates of 20 nstrain/yr of N-S compression in the Catania plain (Mattia et al., 2012). South of this front, in the Hyblean foreland a slight counterclockwise rotation is well evidenced for the first time and rises questions regarding the present-day geodynamics at the western edge of the Ionian subduction (Figure 1).

### 4.2.4 Alpine and Dinaric compression

One striking feature on the divergence map is an arc-shaped compressive limit, roughly parallel to the eastern Alps front, but located 50 km north of the front's surface trace (see Figures 4 and 5). This compressive limit -or zero divergence line- (label 3 in Figure 4) continues eastward in Slovenia and Croatia where it rotates toward a ~N120 direction, and roughly follows the topographic step bounding the Panonnian basin (Bada et al., 2007)(Drava through in profile 4 on Figure 7). To the south in Serbia and Kosovo, this nearly zero divergence contour joins the sharp compressive/dilation limit on the western front of the Albanides (label 2 in Figure 4).

The overall compressive domain including the Po plain, the eastern Alps and the Dinarides does not form a consistent block since tectonic styles strongly vary within it. In particular, the locations of the Rasa strike slip fault in Slovenia and of the N135 Dinarides thrusts along the croatian coast correlate with a very sharp change in rotation direction from anticlockwise to the west to clockwise to the East (see profile 4 in Figure 7 and Figure 4). This confirms the existence of a major plate boundary between the Adria microplate and Eurasia. Along the Grosseto-Budapest profile 5 (Figure 7), $I_2$ slightly increases from the western Dinarides to the longitude of Zagreb and Petrinja, where major earthquakes have occurred recently (Baize et al., 2022; Henriquet et al., 2022b; Herak and Herak, 2023),

while divergence is very low. This area is therefore characterized by a strike-slip tectonic style (Figure 5), consistently with the main active faults (Basili et al., 2022) and observed focal mechanisms. Interestingly, the recent Petrinja earthquake (45.4306°N, 16.2473°E) that happened at the eastern limit of this zone was nearly purely strike-slip while a more composite mechanism was expected (Figure 5). Compression increases again up to 10 nstrain/yr in between the Petrinja fault zone and the Drava through; reverse faulting is promoted there with a N15 compressive principal direction (Figure 5).

### 4.2.5 Spatial extent of the South Balkans Extension Region and major features

Our strain rates calculations allow for better defining the contours of the so-called South Balkans Extension Region (SBER; Burchfiel et al., 2000, 2006), i.e. the wide extensional area where velocities are mostly directed towards the Hellenic subduction zone (with respect to stable Eurasia) and where multiple active faults are reported (Figures 1b and 4). We do constrain the northern limit of the SBER at ~42.5°N by looking at the near zero contour of divergence (Figure 4) and at the map of the tectonic style (Figure 5). However, the SBER extends up to 44°N locally in western Serbia, while it goes only up to 42.5°N in Kosovo, with low extensional values observed down to 41.5°N in central Macedonia. To the East, the Balkans range is included in the SBER, with low dilatation rates observed in the Moesian plain. The Carpathian are also in extension with rates lower than 5 nstrain/yr.

We can have a closer look at the tectonic style and strain rates inside the SBER by looking at two profile lines ranging from the southern Peloponnese to Belgrade (nearly NS profile 2, Figure 6) and from Bari to the Black Sea (nearly EW profile 5, Figure 7). Profile 2 is dominated by the large extension rates observed across the Corinth Gulf (up to 160 nstrain/yr) concentrated along a relatively small area (narrower than 40 km) in between the Kastania and Maratias normal faults. The transition from the Peloponese area straining at around 30 nstrain/yr to the Corinth gulf is very sharp, while a smooth decrease in strain rate is observed along a 75 km long zone to the North before reaching similar levels of strain rates in the Thessalia basin. This dissymmetry could be consistent with the proposed structure of the Corinth Gulf half-graben, which extension is mainly controlled by normal faults located on its southern shoreline (Armijo et al., 1996; Rigo et al., 1996). In the transition zone between the Corinth Gulf and the Thessalia basin, several steps with localized changes in strain rates are observed when crossing the Sperchia, Mavrommati and Rimnio faults either in $I_2$, in the vorticity or in both, suggesting that these fault zones all accommodate part of the current strain accumulation.

Toward North in the Hellenides (profile 2 in Figure 6), $I_2$ remains nearly constant up to the northern end of the Florina basin (41.5°N in Macedonia) where it starts to decrease down to 10 nstrain/yr near Pristhina and nearly zero when crossing the Prokupje strike slip
fault in Serbia. Over this northern end of the SBER zone, the Florina fault zone (D'agostino et al., 2022) correlates with a reduction of the vorticity, while the northern end of the Florina basin may be associated with the eastward prolongation of the Kicevo SW-NE fault (Burchfiel et al., 2006; Dumurdzanov et al., 2005). It is to note that the PDF there appear non gaussian (km 200 to 400), which suggests difficulties in converging and properly recovering surface velocities. This could be due to significantly different velocity fields merged in this region, generating an inconsistent overall description of the current surface straining. The Tetovo normal fault, just west of Skopje, is known to be active (Burchfiel et al., 2006; D'agostino et al., 2022) and is correlated with a sharp change in vorticity, I₂ and divergence along the EW profile 5 presented in Figure 7. It marks the easternmost extent of the highly extensive Albanides (30 nstrain/yr), and the start of the transition from a clockwise rotating domain (Albanides) to a counterclockwise area to the East. In more details, vorticity is almost zero from Skopje to the Macedonia-Bulgaria boundary, just west of the Strouma valley bounded by numerous active faults (Athanassios et al., 2005; Meyer et al., 2002). This observation suggests that this zone, that correlates roughly with the Vardar suture and is the locus of a significant change in tectonic style and deformation rates, could be considered as a boundary between two tectonic domains moving independently. East of the Strouma valley, I₂ decreases under 15 nstrain/yr and reaches finally zero when crossing the Maritsa fault zone (profile 5, Naydenov et al., 2013; Dimitrov et al., 2006). The latter is associated with a ~10 nstrain/yr decrease in both I₂ and vorticity, confirming it as an active structure.

## 5 Conclusion and perspectives: toward a better integration of strain rates in regional geodynamical models and seismic hazard assessment

Reliably recovering surface strain rates has major implications on our understanding of the physical processes controlling the current deformation since it provides key constrains on various modeling strategies.

On the one hand, considering that surface deformation can be modeled at first order by the relative eulerian motion of rigid and elastic blocks bounded by active faults (e.g. McCaffrey, 2009), many authors have used complex block models to describe and recover the surface deformation around the wide aegean region (Nyst and Thatcher, 2004; Floyd et al., 2010; Reilinger et al., 2010; Vernant et al., 2014; Briole et al., 2021; Ergintav et al., 2023), and derive slip rates on main active faults. However, a single block was defined north of the Rhodopes and Hellenides belts, considering the inner Balkans as part of the stable Eurasian plate and failing in describing in details fault rates and tectonic style over this slowly straining region. Our new vorticity and divergence maps exhibit sharp gradients in tectonic style and rotation direction there, that could be used

to draw better defined block boundaries. For instance, the abrupt N-S change in rotation near the Vardar suture zone is a good candidate for a consistent block boundary. Clear changes in rotation rates also appear in central Greece and Peloponese and could serve as objective criteria to define block boundaries there. Such modelling work is beyond the scope of this study but should definitively be conducted in the future.

On the other hand, continuous models of lithospheric straining considering a thin-sheet of viscous lithosphere, submitted to forces as the ones resulting from gravitational potential energy contrasts (GPE), boundary conditions, mantle drag or slab retreat effects, highly depend on surface strain rates observations (Le Pichon and Kreemer, 2010; Perouse et al., 2012; Métois et al., 2015; England et al., 2016; d'Agostino et al., 2020; Daout et al., 2023). These modeling strategies bring constrains on the lithosphere effective viscosity, the direction and amount of the expected stress at the major plate boundaries. As an example, along the North Anatolian Fault, England et al. (2016) chose their best model and the associated power-law rheology for the lithosphere based on a direct comparison of the predicted versus observed width and peak value of the deformation zone experiencing large values of the second invariant. We show in profile 1 (Figure 6) crossing the NAF that B-Strain predicts a width half that of the VISR3D method, with doubled peak values, which could come from significantly different rheological laws than previously modeled.

Finally, we think that this careful mapping and estimate of strain rates is needed to make the best use of the long-term GNSS velocities in seismic hazard calculations. For instance, several recent studies have focused on using GNSS-derived strain rates estimates to compute geodetic moment rates, which can then be compared with effective seismic moment rates (D'Agostino, 2014; Rong et al., 2016; Jenny et al., 2004; Donniol Jouve et al., 2024; Beauval et al., 2018). Even if the Kostrov formulation that allows for geodetic moment rate calculation mainly depends on the considered area and assumed seismogenic depth, it also takes into account the maximum strain rate that has therefore to be correctly determined (Kostrov, 1974). The B-Strain method provides in our opinion reliable estimates of this maximum strain rate, together with a consistent associated uncertainty (i.e. the standard deviation of the PDF) which is often poorly propagated in current interpolation methods. Applied to the most complete existing GNSS velocity field over the peri Adriatic, Balkans and Greece area, B-Strain allows us to propose a novel mapping of current strain rates which uncertainties are realistically estimated. These results, combined with other calculations made on the area using different data sets and methods (e.g. Serpelloni et al., 2022), could therefore be used with confidence in future probabilistic seismic hazard assessment calculations in this seismically very active region of Europe.
## Acknowledgements

This work has been supported by the Programme National de Télédétection Spatiale (PNTS), grant n°PNTS-2019-7; and by the CNES grant "BALKASAR" focused on Sentinel 1-A. The computation of this work is performed on the PSMN platform of ENS de Lyon, funded by the LABEX Lyon Institute of Origins (LIO, ANR-10-LABX-2660066) of University of Lyon. The B-Strain Plotter tool has been developed thanks to the FASCIL project funded by Labex-LIO and led by S.Durand who must be thanked here. It has been completed with the help of the DyoNeo company. We would like to thank V.Guiller and C.Pagani for their initial attempts in vizualisation of B-Strain outputs, and J. Pina-Valdes for making his work easily reusable. We would like to thank A. Socquet and T. Ueda for inspiring discussions. Special thanks to SNCF TER Rhone-Alpes that significantly contributed to this work by randomly extending journey times. The authors thank the reviewers for their careful reading of the manuscript and relevant feedback which helped improve the initial manuscript.

## Author contributions

T.B. designed the inversion method. M.M. conducted the inversion over the Balkans and produced the figures. M.M and C.L designed the study. C.L., T.B., M.H., A.M. and M.M. wrote the manuscript. All authors contributed to the scientific discussion of the results.

## Data availability

The original velocity field comes from supplementary materials from (Piña-Valdés et al., 2022). The B-Strain outputs associated with this work are available on the EaSyData repository, and could be explored online on the BStrain Plotter webtool.

## Competing interests

The authors declare no competing interests.

## Peer review

This publication was peer-reviewed by Corné Kreemer, Kathryn Materna, and Enrico Serpelloni. The full peer-review report can be found here: Review Report.

## Copyright notice

© Author(s) 2025. This article is distributed under the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited, and any changes made are indicated.

## References

Armijo, R., B. Meyer, G. King, A. Rigo, and D. Papanastassiou (1996), Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean, *Geophysical Journal International*, 126(1), 11–53.

Athanassios, G., S. Stefan, D. George, D. Nikolai, S. Sotiris, T. Christina, F. Georgi, and P. Spyros (2005), Active fault segmentation in southwest Bulgaria and Coulomb stress triggering of the 1904 earthquake sequence, *Journal of Geodynamics*, 40(2-3), 316–333.

Bada, G., F. Horváth, P. Dövényi, P. Szafián, G. Windhoffer, and S. Cloetingh (2007), Present-day stress field and tectonic inversion in the Pannonian basin, *Global and Planetary Change*, 58(1-4), 165–180.

Baize, S., S. Amoroso, N. Belić, L. Benedetti, P. Boncio, M. Budić, F. R. Cinti, M. Henriquet, P. J. Rupnik, B. Kordić, et al. (2022), Environmental effects and seismogenic source characterization of the December 2020 earthquake sequence near Petrinja, Croatia, *Geophysical Journal International*.

Basili, R., L. Danciu, C. Beauval, K. Sesetyan, S. P. Vilanova, S. Adamia, P. Arroucau, J. Atanackov, S. Baize, C. Canora, R. Caputo, M. M. C. Carafa, E. M. Cushing, S. Custódio, M. B. Demircioglu Tumsa, J. C. Duarte, A. Ganas, J. García-Mayordomo, L. Gómez de la Peña, E. Gràcia, P. Jamšek Rupnik, H. Jomard, V. Kastelic, F. E. Maesano, R. Martín-Banda, S. Martínez-Loriente, M. Neres, H. Perea, B. Sket-Motnikar, M. M. Tiberti, N. Tsereteli, V. Tsironi, R. Vallone, K. Vanneste, and P. Zupančič (2022), European Fault-Source Model 2020 (EFSM20): online data on fault geometry and activity parameters, doi: 10.13127/EFSM20.

Baxter, S. C., S. Kedar, J. W. Parker, F. H. Webb, S. E. Owen, A. Sibthorpe, and D. Dong (2011), Limitations of strain estimation techniques from discrete deformation observations, *Geophysical Research Letters*, 38(1).

Beauval, C., J. Marinière, H. Yepes, L. Audin, J.-M. Nocquet, A. Alvarado, S. Baize, J. Aguilar, J.-C. Singaudo, and H. Jomard (2018), A new seismic hazard model for Ecuador, *Bulletin of the Seismological Society of America*, 108(3A), 1443–1464.

Bodin, T., and M. Sambridge (2009), Seismic tomography with the reversible jump algorithm, *Geophysical Journal International*, 178(3), 1411–1436.

Briole, P., A. Ganas, P. Elias, and D. Dimitrov (2021), The GPS velocity field of the Aegean. New observations, contribution of the earthquakes, crustal blocks model, *Geophysical Journal International*, 226(1), 468–492.

Burchfiel, C., R. W. King, A. Todosov, V. Kotzev, N. Durmurdzanov, T. Serafimovski, and B. Nurce (2006), GPS results for Macedonia and its importance for the tectonics of the Southern Balkan extensional regime, *Tectonophysics*, 413(3), 239–248.

Burchfiel, C. B., R. Nakov, T. Tzankov, and L. H. Royden (2000), Cenozoic extension in Bulgaria and northern Greece: the northern part of the Aegean extensional regime, *Geological Society, London, Special Publications*, 173(1), 325–352.

D'Agostino, N. (2014), Complete seismic release of tectonic strain and earthquake recurrence in the Apennines (Italy), *Geophysical Research Letters*, 41(4), 1155–1162.

D'Agostino, N., A. Avallone, D. Cheloni, E. D'anastasio,
S. Mantenuto, and G. Selvaggi (2008), Active tectonics of the Adriatic region from GPS and earthquake slip vectors, *Journal of Geophysical Research: Solid Earth* (1978–2012), 113(B12).

d'Agostino, N., M. Métois, R. Koci, L. Duni, N. Kuka, A. Ganas, I. Georgiev, F. Jouanne, N. Kaludjerovic, and R. Kandić (2020), Active crustal deformation and rotations in the southwestern Balkans from continuous GPS measurements, *Earth and Planetary Science Letters*, 539, 116, 246.

Daout, S., N. d'Agostino, E. Pathier, A. Socquet, J. Lavé, M.-P. Doin, M. Riesner, and L. Benedetti (2023), Along-strike variations of strain partitioning within the Apennines determined from large-scale multi-temporal InSAR analysis, *Tectonophysics*, 867, 230, 076.

Dimitrov, D., T. Camelbeek, J. Ruegg, I. Georgiev, and E. Botev (2006), Surface seismic deformations in the Plovdiv Region (Bulgaria) by space geodesy and seismology data, *SENS*, pp. 14–16.

Donniol Jouve, B., A. Socquet, C. Beauval, J. Piña Valdès, and L. Danciu (2024), Consistency between the Strain Rate Model and ESHM20 Earthquake Rate Forecast in Europe: insights for seismic hazard, *EGUsphere*, 2024, 1–32.

Dumurdzanov, N., T. Serafimovski, and B. C. Burchfiel (2005), Cenozoic tectonics of Macedonia and its relation to the South Balkan extensional regime, *Geosphere*, 1(1), 1–22.

D'agostino, N., A. Copley, J. Jackson, R. Koçi, A. Hajrullai, L. Duni, and N. Kuka (2022), Active tectonics and fault evolution in the Western Balkans, *Geophysical Journal International*, 231(3), 2102–2126.

England, P., G. Houseman, and J.-M. Nocquet (2016), Constraints from GPS measurements on the dynamics of deformation in Anatolia and the Aegean, *Journal of Geophysical Research: Solid Earth*, 121(12), 8888–8916.

Ergintav, S., M. Floyd, D. Paradissis, H. Karabulut, P. Vernant, F. Masson, I. Georgiev, A. Ö. KONCA, U. DOĞAN, R. King, et al. (2023), New geodetic constraints on the role of faults and blocks vs. distribute strain in the Nubia-Arabia-Eurasia zone of active plate interactions, *Turkish Journal of Earth Sciences*, 32(3), 248–261.

Floyd, M., H. Billiris, D. Paradissis, G. Veis, A. Avallone, P. Briole, S. McClusky, J.-M. Nocquet, K. Palamartchouk, B. Parsons, et al. (2010), A new velocity field for Greece: Implications for the kinematics and dynamics of the Aegean, *Journal of Geophysical Research: Solid Earth* (1978–2012), 115(B10).

Ganas, A., V. Tsironi, E. Kollia, M. Delagas, C. Tsimi, and A. Oikonomou (2018), Recent upgrades of the NOA database of active faults in Greece (NOAFAULTS), *Proceedings of the 19th General Assembly of WEGENER, Grenoble, France*, pp. 10–13.

Gordon, R. G., and G. A. Houseman (2015), Deformation of Indian Ocean lithosphere: Evidence for a highly nonlinear rheological law, *Journal of Geophysical Research: Solid Earth*, 120(6), 4434–4449.

Grosset, J., S. Mazzotti, and P. Vernant (2023), Glacial-isostatic-adjustment strain rate–stress paradox in the Western Alps and impact on active faults and seismicity, *Solid Earth*, 14(10), 1067–1081.

Harvard (), *Global Centroid Moment Tensor Catalog*, http://www.globalcmt.org/CMTsearch.html.

Henriquet, M., M. Peyret, S. Dominguez, G. Barreca,

C. Monaco, and S. Mazzotti (2022a), Present-day surface deformation of Sicily derived from Sentinel-1 InSAR time-Series, *Journal of Geophysical Research: Solid Earth*, 127(3), e2021JB023,071.

Henriquet, M., B. Kordic, M. Métois, C. Lasserre, S. Baize, L. Benedetti, M. Spelić, and M. Vukovski (2022b), Rapid remeasure of dense civilian networks as a game-changer tool for surface deformation monitoring: The case study of the Mw 6.4 2020 Petrinja Earthquake, Croatia, *Geophysical Research Letters*, 49(24), e2022GL100,166.

Herak, M., and D. Herak (2023), Properties of the Petrinja (Croatia) earthquake sequence of 2020–2021–Results of seismological research for the first six months of activity, *Tectonophysics*, 858, 229, 885.

Houseman, G., and P. England (1986), Finite strain calculations of continental deformation: 1. Method and general results for convergent zones, *Journal of Geophysical Research: Solid Earth*, 91(B3), 3651–3663.

Jenny, S., S. Goes, D. Giardini, and H.-G. Kahle (2004), Earthquake recurrence parameters from seismic and geodetic strain rates in the eastern Mediterranean, *Geophysical Journal International*, 157(3), 1331–1347.

Jolivet, L., C. Faccenna, B. Huet, L. Labrousse, L. Le Pourhiet, O. Lacombe, E. Lecomte, E. Burov, Y. Denèle, J.-P. Brun, et al. (2013), Aegean tectonics: Strain localisation, slab tearing and trench retreat, *Tectonophysics*, 597, 1–33.

Kostrov, V. (1974), Seismic moment and energy of earthquakes, and seismic flow of rock, *Physics of the Solid Earth*, 1, 13–21.

Kreemer, C., and N. Chamot-Rooke (2004), Contemporary kinematics of the southern Aegean and the Mediterranean Ridge, *Geophysical Journal International*, 157(3), 1377–1392.

Kreemer, C., N. Chamot-Rooke, and X. Le Pichon (2004), Constraints on the evolution and vertical coherency of deformation in the Northern Aegean from a comparison of geodetic, geologic and seismologic data, *Earth and Planetary Science Letters*, 225(3-4), 329–346.

Le Pichon, X., and C. Kreemer (2010), The Miocene-to-present kinematic evolution of the Eastern Mediterranean and Middle East and its implications for dynamics, *Annual Review of Earth and Planetary Sciences*, 38, 323–351.

Mathey, M., C. Sue, C. Pagani, S. Baize, A. Walpersdorf, T. Bodin, L. Husson, E. Hannouz, and B. Potin (2021), Present-day geodynamics of the Western Alps: new insights from earthquake mechanisms, *Solid Earth*, 12(7), 1661–1681.

Mattia, M., V. Bruno, F. Cannavò, and M. Palano (2012), Evidences of a contractional pattern along the northern rim of the Hyblean Plateau (Sicily, Italy) from GPS data, *Geologica Acta: an international earth science journal*, 10(1), 1–8.

Mazzotti, S., T. S. James, J. Henton, and J. Adams (2005), GPS crustal strain, postglacial rebound, and seismic hazard in eastern North America: The Saint Lawrence valley example, *Journal of Geophysical Research: Solid Earth*, 110(B11).

McCaffrey, R. (2009), Time-dependent inversion of three-component continuous GPS for steady and transient sources in northern Cascadia, *Geophysical Research Letters*, 36(7).

Metois, M., N. D'Agostino, A. Avallone, N. Chamot-Rooke,
A. Rabaute, L. Duni, N. Kuka, R. Koci, and I. Georgiev (2015), Insights on continental collisional processes from GPS data: Dynamics of the peri-Adriatic belts, *Journal of Geophysical Research: Solid Earth*, 120(12), 8701–8719.

Meyer, B., R. Armijo, and D. Dimitrov (2002), Active faulting in SW Bulgaria: possible surface rupture of the 1904 Struma earthquakes, *Geophysical Journal International*, 148(2), 246–255.

Moulin, A., L. Benedetti, A. Gosar, P. J. Rupnik, M. Rizza, D. Bourlès, and J.-F. Ritz (2014), Determining the present-day kinematics of the Idrija fault (Slovenia) from airborne LiDAR topography, *Tectonophysics*, 628, 188–205.

Naydenov, K., I. Peytcheva, A. von Quadt, S. Sarov, K. Kolcheva, and D. Dimov (2013), The Maritsa strike-slip shear zone between Kostenets and Krichim towns, South Bulgaria—Structural, petrographic and isotope geochronology study, *Tectonophysics*, 595, 69–89.

Nocquet, J.-M. (2012), Present-day kinematics of the Mediterranean: A comprehensive overview of GPS results, *Tectonophysics*, 579, 220–242.

Nucci, R., E. Serpelloni, L. Faenza, A. Garcia, and M. E. Belardinelli (2023), Comparative Analysis of Methods to Estimate Geodetic Strain Rates from GNSS Data in Italy, *Annals of Geophysics*.

Nyst, M., and W. Thatcher (2004), New constraints on the active tectonic deformation of the Aegean, *Journal of Geophysical Research: Solid Earth*, 109(B11).

Okada, Y. (1985), Surface deformation due to shear and tensile faults in a half-space, *Bulletin of the seismological society of America*, 75(4), 1135–1154.

Pagani, C., T. Bodin, M. Métois, and C. Lasserre (2021), Bayesian Estimation of Surface Strain Rates From Global Navigation Satellite System Measurements: Application to the Southwestern United States, *Journal of Geophysical Research: Solid Earth*, 126(6), e2021JB021,905.

Palano, M., L. Ferranti, C. Monaco, M. Mattia, M. Aloisi, V. Bruno, F. Cannavò, and G. Siligato (2012), GPS velocity and strain fields in Sicily and southern Calabria, Italy: Updated geodetic constraints on tectonic block interaction in the central Mediterranean, *Journal of Geophysical Research: Solid Earth*, 117(B7).

Perouse, E., N. Chamot-Rooke, A. Rabaute, P. Briole, F. Jouanne, I. Georgiev, and D. Dimitrov (2012), Bridging onshore and offshore present-day kinematics of central and eastern Mediterranean: Implications for crustal dynamics and mantle flow, *Geochemistry, Geophysics, Geosystems*, 13(9).

Piña-Valdés, J., A. Socquet, C. Beauval, M.-P. Doin, N. D'Agostino, and Z.-K. Shen (2022), 3D GNSS velocity field sheds light on the deformation mechanisms in Europe: Effects of the vertical crustal motion on the distribution of seismicity, *Journal of Geophysical Research: Solid Earth*, 127(6), e2021JB023,451.

Reilinger, R., S. McClusky, D. Paradissis, S. Ergintav, and P. Vernant (2010), Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone, *Tectonophysics*, 488(1), 22–30.

Rigo, A., H. Lyon-Caen, R. Armijo, A. Deschamps, D. Hatzfeld, K. Makropoulos, P. Papadimitriou, and I. Kassaras (1996), A microseismic study in the western part of the Gulf of Corinth (Greece): implications for large-scale normal faulting mechanisms, *Geophysical Journal International*, 126(3), 663–688.

Rong, Y., P. Bird, and D. Jackson (2016), Earthquake potential and magnitude limits inferred from a geodetic strain-rate model for southern Europe, *Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society*, 205(1), 509–522.

Sambridge, M., T. Bodin, K. Gallagher, and H. Tkalčić (2013), Transdimensional inference in the geosciences, *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 371(1984), 20110,547.

Serpelloni, E., M. Anzidei, P. Baldi, G. Casula, and A. Galvani (2005), Crustal velocity and strain-rate fields in Italy and surrounding regions: new results from the analysis of permanent and non-permanent GPS networks, *Geophysical Journal International*, 161(3), 861–880.

Serpelloni, E., A. Cavaliere, L. Martelli, F. Pintori, L. Anderlini, A. Borghi, D. Randazzo, S. Bruni, R. Devoti, P. Perfetti, et al. (2022), Surface velocities and strain-rates in the Euro-Mediterranean region from massive GPS data processing, *Frontiers in Earth Science*, 10, 907,897.

Shen, Z.-K., M. Wang, Y. Zeng, and F. Wang (2015), Optimal interpolation of spatially discretized geodetic data, *Bulletin of the Seismological Society of America*, 105(4), 2117–2127.

Tape, C., P. Musé, M. Simons, D. Dong, and F. Webb (2009), Multiscale estimation of GPS velocity fields, *Geophysical Journal International*, 179(2), 945–971.

USGS (2025), Earthquake Catalog.

Ventura, B. M., E. Serpelloni, A. Argnani, A. Bonforte, R. Bürgmann, M. Anzidei, P. Baldi, and G. Puglisi (2014), Fast geodetic strain-rates in eastern Sicily (southern Italy): New insights into block tectonics and seismic potential in the area of the great 1693 earthquake, *Earth and planetary science letters*, 404, 77–88.

Vernant, P., R. Reilinger, and S. McClusky (2014), Geodetic evidence for low coupling on the Hellenic subduction plate interface, *Earth and Planetary Science Letters*, 385, 122–129.

Walpersdorf, A., L. Pinget, P. Vernant, C. Sue, A. Deprez, and R. team (2018), Does long-term GPS in the Western Alps finally confirm earthquake mechanisms?, *Tectonics*, 37(10), 3721–3737.