Postglacial strain rate - stress paradox, example of the Western Alps active faults

. The understanding of the origins of seismicity in intraplate regions is crucial to better characterize seismic hazards. In formerly glaciated regions such as Fennoscandia North America or the Western Alps, stress perturbations from Glacial Isostatic Adjustment (GIA) have been proposed as a major cause of large earthquakes. In this study, we focus on the Western Alps case using numerical modeling of lithosphere response to the Last Glacial Maximum icecap. We show that the ﬂexural response to GIA induces present-day stress perturbations of ca. 1-2 MPa, associated with horizontal extension rates up to 5 ca. 2 . 5 × 10 − 9 yr -1 . The latter is in good agreement with extension rates of ca. 2 × 10 − 9 yr -1 derived from high-resolution geodetic (GNSS) data and with the overall seismicity deformation pattern. In the majority of simulations, stress perturbations induced by GIA promote fault reactivation in the internal massifs and in the foreland regions (i.e., positive Coulomb Failure Stress perturbation), but with predicted rakes systematically incompatible with those from earthquake focal mechanisms. Thus, although GIA explains a major part of the GNSS strain rates, it tends to inhibit the observed seismicity in the Western Alps. 10 A direct corollary of this result is that, in cases of signiﬁcant GIA effect, GNSS strain rate measurements cannot be directly integrated in seismic hazard computations, but instead require detailed modeling of the GIA transient impact.

1˘2 × 10 −9 yr -1 . As shown by the original velocities (Fig. 3a), these signals are at the limit of GNSS detection and resolution and can only be extracted from the background noise using spatial filtering / smoothing technics. The results shown in Fig. 3 are based on a spatial Gaussian filtering with a half-width of 90 km (cf., Masson, 2019), which allows highlighting coherent strain rate signals at wavelength of ca. 150-200 km that we consider appropriate for comparison with ongoing GIA deformation due to the Alpine Icecap (cf. Fig. 3 and next section). Details of the Gaussian filtering method and additional strain rate maps at 70 different filtering half-widths can be found in the Supplementary Material A2 to provide an overview of the strain rate patterns and their sensibility to the filtering wavelength.
This overall deformation pattern (slow horizontal extension / extension rates coupled with relatively fast uplift rates) is at the core of the current debate on the processes responsible for the ongoing Western Alps deformation. The role of plate tectonics being, at best, very small, alternative processes include mantle and slab-tear dynamics (Sternai et al., 2019), surface erosion 75 (Vernant et al., 2013) and Glacial Isostatic Adjustment Mey et al., 2016). Numerical modeling shows that the isostatic response to erosion can generate the observed extensional strain pattern (Vernant et al., 2013), but the associated uplift rates are significantly smaller than the GNSS velocities (Sternai et al., 2019). On the contrary, GIA from the Last Glacial Maximum can explain the present-day uplift rates Mey et al., 2016) but a detailed comparison with horizontal deformation is lacking. In the following section, we test the compatibility of GIA models with present-day 80 deformation of the Western Alps

GIA models
In order to compare GIA effects with GNSS and seismicity observations, we model the response of the lithosphere to ice loading and unloading using a simple approach of an elastic thin plate over a viscous substratum. The elastic plate flexure wm in response to the ice load is computed with the gFlex code (Wickert, 2016) using the Kirchhoff-Love thin-plate analytic 85 solution (assuming small inline deformation, zero vertical shear, and plate rotation lower than 10 • ): with kei the zero-order Kelvin function, α the flexural parameter (function of the elastic plate thickness T e and elastic parameters) and w 0m the maximum flexure at the center of the load (cf. Turcotte and Schubert, 2014;Wickert, 2016). The ice load model at LGM is from Mey et al. (2016): ca. 600 km long by 200 km wide, locally up to 2 km thick (cf. Fig. 3).
The postglacial rebound response is computed assuming an instantaneous ice melt after the Last Glacial Maximum (15000 yr BP) and a simple exponential decay associated with a linear viscosity of relaxation time τ g (Turcotte and Schubert, 2014): This simple 2D planar model is warranted by the relatively small spatial scale of the Alpine ice load. It comprises only two parameters: The flexural parameter α, which depends mainly on the elastic plate thickness T e , and the characteristic time 95 τ g . T e controls the spatial wavelength of the area affected by subsidence, the location of the forebulge and the maximum subsidence value (a higher T e results in a wider a subsidence area with a lower amplitude). The characteristic time τ g controls the relaxation speed, i.e., the time to come back to the initial pre-loading state (a higher τ g results in a longer the postglacial rebound).
In order to test the variability of GIA predictions to these two parameters, we compute models with 5 ≤ T e ≤ 100 km and 100 2000 ≤ τ g ≤ 20000 yr. Assuming that GNSS uplift rates are mostly due to GIA, comparisons between observations and model predictions indicate a best-fit parameter set ca. T e = 10-20 km and τ g = 4500-5500 yr (F-test statistics on residuals). Because GNSS data may be impacted by other processes (cf., above), we allow for a larger range of parameter values 5 ≤ T e ≤ 40 km and 3000 ≤ τ g ≤ 7000 yr. These values indicate a fast and short wavelength response (low mantle velocity coupled with a thin elastic plate), in agreement with more complex models (e.g" spatially variable plate thickness) tuned to the Western Alps 105 Mey et al., 2016).
Stress, strain and strain rates in the elastic plate are computed from the second spatial derivative of the flexural response w(x, y, t) (eq. 2) based on the Kirchhoff-Love plate theory (no vertical shear stress, only in-plane compression / tension stress).
Although the actual depth and structure of the lithosphere elastic core is unknown, standard rheological models (e.g., Burov and Diament, 1995) and the shallow depth of earthquakes in the Western Alps (less than 15km; Cushing et al., 2008;Delacou 110 et al., 2004) indicate that the seismogenic depth range is likely associated with the upper section of the model elastic plate.
In order to consider the maximum GIA stress perturbations and strain rate signals, we consider hereafter the stress, strain and strain rate values computed at the top of the elastic plate. More complex situation with a decoupled two-layer elastic plate system would need to be tested with alternative models (e.g., finite element approach, Cattin et al., 2001), but we consider such complexity unlikely in the context of a thin elastic plate such as the Western Alps. the modeled values range are about 0.9˘2.1 × 10 −9 yr -1 while the GNSS data indicates (2.5 ± 1) × 10 −9 yr -1 . In the French Alps, the model rates range between 0.4 and 0.8 × 10 −9 yr -1 in agreement with the (0.9 ± 0.7) × 10 −9 yr -1 given by the GNSS measurements. Based on these results, GIA accounts for a minimum of 30% and up to 100% of the GNSS extension rates.
In contrast, the GIA models are unable to reproduce the GNSS (0.8 ± 1) × 10 −9 yr -1 E-W extension rates in Italy or the (1.5 ± 1) × 10 −9 yr -1 strike slip strain rates along the southern France-Italy border. This area being located south of the Adria /

130
Eurasia rotation pole, we propose that these strain rates may be in part related to Adria micro-plate kinematics, consistent with the 0.3 − 0.8 mm.yr -1 of E-W extension in southern Western Alps (e.g., Nocquet, 2012).
The inner Western Alps extension is surrounded by radial shortening in the Rhine Graben, Jura Mountains and Rhone Valley regions. In the Rhine graben and Jura, the GIA modeled shortening rate directions are consistent with the GNSS strain rates, but the amplitudes are about three times smaller. This suggests that the shortening rates in these regions cannot be solely attributed 135 to Alpine GIA; processes such as Scandinavian postglacial rebound (Nocquet et al., 2005) or a buoyant mantle plume located beneath Eifel volcanic area (Kreemer et al., 2020) may contribute to the shortening signal.

GIA stress perturbations on active faults
As expected from the conceptual flexural model (Fig. 1), the GIA stress perturbations correspond to horizontal compressive stress: uniaxial NW-SE compression in southern Switzerland, near-isotropic compression on the French Alps, and E-W uniaxial 140 compression in the Southern Alps (French-Italian border). As previously mentioned, strain rate and stress present opposite signs and these maximum GIA compression (Fig. 4a) are located at the maximum GIA extension rate (Fig. 3c).
The impact of GIA stress perturbations on faults is estimated using the Coulomb failure criterion based on the fault azimuth and dip angle. The Coulomb Failure Stress ∆ CF S is computed based on the fault shear and normal stress (e.g., King et al., 1994): with µ the effective coefficient of friction (between 0.1 and 0.6), τ the shear stress and σ n the normal stress. The fault-slip style associated with this stress perturbation is given by the rake r: with τ d and τ s the fault shear stress along dip and azimuth directions, respectively. In all models, these Coulomb Failure at a given stress level if ∆ CF S < 0, or GIA inhibiting fault motion if ∆ CF S > 0. Figure 4b presents the computed rakes for three sets of typical faults in the Western Alps (Fig. 4a): the Vuache Fault, the Belledonne Fault, and a series of potential faults representative of normal-faulting earthquakes in the Internal Massifs above the Penninic Front (hereafter IMNF, Inner Massifs Normal Faults; Sue et al., 2003Sue et al., , 1999. The fault geometries are extracted 155 from the BDFA and neotectonic databases (Grellet, 1993;Jomard et al., 2017). Due to the uncertainties on the dip angles, several values are tested (30-80 • ). In Figure 4b, the GIA model rakes are shown as boxplots that represent the variability in GIA model parameters (T e and τ g ) and fault geometries (azimuth and dip). The predicted rakes compared to those derived from earthquake focal mechanisms.

160
The GIA Coulomb Failure Stress perturbation depends on the fault geometry and the assumed effective friction (eq. 3).
Models with standard friction (µ = 0.6) and near-vertical fault dips results in negative ∆ CF S in dictating that GIA tends to inhibit fault rupture in these cases. In contrast, models with low effective friction (µ = 0.1) or low fault dips result in mostly positive ∆ CF S between 0.1 and 2 MPa. These models are associated with fault-slip rakes corresponding to reverse faulting (rakes r = 50-110, Fig. 4b). This faulting style is opposite to that observed in the focal mechanisms along the IMNF (Delacou  Thouvenot et al., 2003). The Vuache Fault orientation is strongly oblique to the GIA stress field (Fig. 4a), with the northern part of the fault located in a transition zone where the GIA stress is almost zero. Predicted rakes vary between 60 and 150 (reverse faulting), inconsistent with the observed left-lateral faulting style (Rabin et al., 2018). In these three cases, GIA ∆ CF S perturbations can reach up to 1-2 MPa, 170 a level commonly considered sufficient to trigger preloaded fault rupture (e.g., King et al., 1994;Steer et al., 2014;Stein, 1999), but the modeled rakes are systematically different from the observed ones. Thus, GIA in the Western Alps cannot explain the kinematics of the aforementioned faults.
The effect of GIA is more difficult to identify for structures parallel to the main stress direction. Along the Durance Fault, in the southern Western Alps, the GIA extensive NE-SW stress perturbations decrease from 0.35-0.7 MPa in the area north of 175 the fault to 0.15-0.25 MPa around the southern segment (Fig. 5a). Due to the bends in the Durance Fault trace, the range of predicted rakes is very large for this fault (Fig. 5b). ∆ CF S are systematically positive (with low effective friction µ = 0.1).
For the northern half, GIA modeled rakes vary from 0 • to 180 • indicating that, depending on the fault segment orientation, the kinematics can be reverse, left-lateral, or right-lateral. For the southern half, GIA stress perturbations predict a right-lateral strike-slip motion. The focal mechanisms along the Durance Fault system suggest mostly left-lateral displacements (Cushing 180 et al., 2008). Thus, rupture on some segments may be promoted by the GIA stress perturbations, whereas others are clearly inhibited (Fig. 5b).

185
These deformation rates are compatible with GIA signal from the Last Glacial Maximum provided the Western Alps behave mechanically as a thin elastic plate (T e = 10-35 km) over a low viscosity mantle ( τ g = 4500-8000 yr) (e.g., Chéry et al., 2016;Mey et al., 2016, cf. Fig .3). Although the actual contribution of GIA relative to other processes (e.g., slab breakoff, erosion) is debated (Sternai et al., 2019), the good agreement between simple GIA model predictions and the observed GNSS uplift rates and horizontal strain rates suggests that GIA is a major contributor to those observed deformation rates, at least in the region 190 below the LGM icecap. The compatibility between the GNSS strain rates and the seismicity deformation style has led several authors to consider GNSS data as a constraint on seismicity budget (e.g., Sue et al., 2007;Walpersdorf et al., 2018Walpersdorf et al., , 2015. Our study shows an apparent strain rate -stress paradox, wherein GIA can explain the observed most of GNSS strain rates but not the seismicity. This apparent paradox is easily resolved by considering that the observed (and modeled) horizontal extension rates are merely a diminution over time of the compression / shortening induced at surface of the lithosphere by its 200 downward flexure under the LGM ice load. This situation is similar to that observed in the forearcs of subduction zones where a stress orientations and styles may differ from the strain rate patterns, the latter being transient signals due to the interseismic locking of the subduction fault (Wang, 2000).
A similar situation has been studied by Craig et al. (2016) for Fennoscandia and its early-postglacial reverse-faulting earthquakes. They argue that these reverse earthquakes occurred during a period of GIA extension strain rates, and thus were the 205 expression of compressive stress stored in the lithosphere over a long (tectonic) time. This analysis relies on the assumption that faults are sensitive to the stress rate of change since the LGM (reduction of compressive stress = extensional stress rate), when in fact the full GIA stress perturbation remains compressive throughout the whole glacial cycle (up to present-day), thus promoting reverse faulting earthquakes as proposed by previous studies (Hampel et al., 2009;Steffen et al., 2014;Wu et al., 1999).

210
Coming back to the Western Alps case, the opposite style between the GIA strain rates and stress perturbations results in a particular situation. GIA is likely the cause of a large part (over 50%) of the observed geodetic strain rates, but it tends to inhibit the current seismicity and its arc-perpendicular extension. Thus, although it may be tempting to do so, GNSS strain rates cannot be used in seismicity rates or seismic hazard analyses in this context. Our results point out the necessity for an additional mechanism (or several additional mechanisms) that can be the cause of the arc-perpendicular extension pattern expressed in the extension stress pattern can be modulated by the GIA stress perturbations, with a reduction of the deformation and earthquake activity during the glacial period and an increase in activity after the glacial maximum (cf. Hampel et al., 2009;Steffen et al., 2014). The integration of GIA stress perturbations in seismicity and seismic hazard studies would thus require mechanical models and cannot be simply based on direct GNSS strain rate integration.

220
Code and data availability. Maps have been made with GMT6. Figure 4b and 5b were performed with python 3.6 (Python Software Foundation), as well as models presented in the article. GPS velocities field has been processed by C. Masson using CCRS-PPP (details in Masson (2019)).Flexure is computed from gFlex code (Wickert, 2016) and glex_load code (https://github.com/juliettegroset/gflex_load.git). GPS velocities field has been processed by C. Masson using CCRS-PPP (details in Masson (2019)).

A1 Horizontal Velocity field
The velocity field (Fig. A1) is from Masson (2019). The analysis has been made on 1443 permanent stations spread on western Europe in Precise Point Positioning (PPP) from CSRS-PPP software (Kouba and Héroux, 2001). The detail of the method is described in Masson (2019).

A2 Gaussian filtering method 230
The smoothed velocity field is computed on a regular grid (0.5 × 0.5 degree). The velocity on each grid point is computed from velocity vectors of all GNSS stations, weighted according to their distance with a Gaussian function (Mazzotti et al., 2011).
The velocity vector at each point of the grid V g is defined by: With i the component (North, East, Up), N the number of stations, V n the velocity vector at the station n and σ the uncer-235 tainty of the velocity component. G n is defined by the Gaussian function parameters r g (the Gaussian half-width), and ∆ n the distance between stations and grid point: The velocity Gaussian function grid is differentiated to obtain the strain rate tensor at each point. The effect of Gaussian half-width r g on the noise in the final solution is shown in fig. A2.