Glacial isostatic adjustment strain rate - stress paradox in the Western Alps, impact on active faults and seismicity

. In many regions formerly glaciated during the Last Glacial Maximum (LGM), Glacial Isostatic Adjustment (GIA) explains most of the measured uplift and deformation rates. GIA is also proposed as a key process contributing to fault activity and seismicity shortly after the LGM and potentially up to present-day. Here, we study the impact of GIA on present-day fault activity and seismicity in the Western Alps. We show that, in the upper crust, GIA induces horizontal compressive stress perturbations associated with horizontal extension rates. The latter agree with the observed geodetic strain rates and with the 10 seismicity deformation patterns. Yet, in nearly all cases, the GIA stress perturbations tend either to inhibit fault slip, or to promote fault slip with the wrong mechanism compared to the seismicity deformation style. Thus, although GIA from the LGM explains a major part of the Western Alps geodetic strain rates, it does not drive nor promote the observed seismicity (which must be driven by other processes). This apparent strain rate - stress paradox results from the gradual diminution over time of the finite shortening induced in the upper crust by the Würm icecap load. A direct corollary of our results is that 15 seismicity and seismic hazard studies in the Western Alps cannot directly integrate geodetic velocities and strain rates, but instead require detailed modeling of the GIA transient impact.


Introduction -GIA, deformation and seismicity
Glacial Isostatic Adjustment (GIA) is the mechanical response of the Earth's crust and mantle to loading / unloading cycles of continental ice sheets and glaciers.The associated surface deformation is observed in geomorphological features such as raised paleo-shorelines and in geodetic measurements of present-day uplift rates or horizontal strain rates (cf.review in Peltier et al., 2022).The impact of GIA and the associated stress perturbations on fault activity and seismicity has been studied since the seminal work of Johnston (1987), Quinlan (1984) and Stephansson (1988), up to recent developments including more complex Earth and fault mechanics models (cf.Steffen et al., 2021).A common feature of these studies is the demonstration that GIA increases fault activity and seismicity shortly after the main deglaciation phase in many regions formerly glaciated during the Last Glacial Maximum (LGM) (eg., Wu et al., 1999;Stewart et al., 2000;Muir-Wood, 2000;Hetzel and Hampel, 2005;Steffen et al., 2014;Lund, 2015).In contrast, the potential effects of GIA stress perturbations on present-day faulting and seismicity in and near formerly glaciated regions remains debated (e.g., Bungum et al., 2010;Bungum and Eldholm, 2022;Brandes et al., 2015).
The first-order mechanics of GIA-related seismicity involves lithosphere flexure coupled with mantle relaxation in response to ice loading and unloading.Associated stress perturbations in the upper crust can reach a few MPa, sufficient to induce rupture on faults near failure equilibrium or to unclamp faults thus allowing the release of long-term stored tectonic stress (e.g., Craig et al., 2016;Hetzel and Hampel, 2005;Steffen et al., 2014).These effects are based on the same model (Fig. 1), wherein the ice loading results in a downward flexure and a horizontal compressive stress perturbation in the upper half of the elastic lithosphere beneath the load (resp.upward flexure / extensive stress in the forebulge regions).Following deglaciation, the lithosphere unbending is dampened by the mantle viscosity resulting in a gradual diminution of the initial bending stress.The associated strain corresponds to a maximum shortening at the peak of glaciation followed by a diminution of shortening during the postglacial phase.As a result, the present-day surface strain rate corresponds to an extension rate while the stress perturbation remains compressive.This apparent contradiction is similar to the strain rate -stress paradox observed in subduction zone forearcs in relation with transient interseismic locking of the megathrust fault (Wang, 2000).At the end of the GIA cycle, bending stress are fully released and the plate regains its initial background state of stress (plus the potential steady-state tectonic loading) (Fig. 1).
Owing to GNSS (Global Navigation Satellite Network) data, present-day horizontal extension rates due to GIA have been measured with increasing accuracies not only in regions formerly covered by Wisconsin and Weichselian ice sheets (Calais et al., 2006;Keiding et al., 2015;Tarayoun et al., 2018), but also in regions of smaller icecaps and mountains glaciers, such as the Würm glaciation in the European Alps (Masson et al., 2019;Nguyen et al., 2016;Walpersdorf et al., 2018).These horizontal strain rates have been compared with seismicity deformation patterns and rates, with various degrees of agreement in both styles and amplitudes (Keiding et al., 2015;Mazzotti et al., 2005;Sánchez et al., 2018).One of the goals of these comparisons is the integration of GNSS velocities and strain rates in seismic hazard models (Mathey et al., 2020), assuming that these short-term (5-20 yr) data can provide information on longer-term (10 2 -10 5 yr) earthquake activity.However, the apparent GIA strain rate -stress paradox puts strong doubts on the pertinence of using GNSS rates to compare with seismicity in regions affected by ongoing transient GIA deformation (even more so for regions with very low tectonic loading rates such as intraplate domains).
In this study, we compare GIA deformation and stress predictions with present-day GNSS strain rates, fault activity and seismicity in the Western Alpes (France, Italy, Switzerland), where GIA contributes to a large part of the measured geodetic deformation rates (Sternai et al., 2019;Stocchi et al., 2005).Specifically, we analyze the effect of GIA stress perturbations on a series of typical fault systems of the Western Alps, in order (1) to compare with their observed kinematics and associated earthquake focal mechanisms, and (2) to test whether GIA stresses tend to promote or inhibit the present-day seismicity.
The present-day geodynamics of the Central and Eastern Alps are primarily controlled by the counter-clockwise rotation of the Adria microplate relative to Eurasia, with a rotation pole east of the Western Alps near Turin, northwestern Italy (Fig. 2a; Battaglia et al., 2004;D'Agostino et al., 2008).However, this rotation kinematics is incompatible with the seismicity and geodetic horizontal deformation patterns in the Western Alps (Fig. 2.b).GNSS data indicate that the overall horizontal deformation across the Western Alps is smaller than their current resolution, i.e., less than 0.3 mm.yr -1 relative motion between northwestern Italy (Pô Plain) and the Rhône Valley -Eastern France (Masson et al., 2019;Sánchez et al., 2018).
Geodetic data also indicate that the Western Alps are affected by a significant regional uplift rate between 0.5 and 2.5 mm.yr -1 peaking in the inner high-altitude areas and decreasing towards 0 mm.yr -1 in the surrounding lowlands (Fig. 3a; Brockmann et al., 2012;Masson et al., 2019;Nguyen et al., 2016;Nocquet et al., 2016).The highest uplift rates roughly coincide with horizontal arc-perpendicular extension rates ca. 2 x 10 -9 yr -1 in the inner regions (Switzerland, French-Italian border, Fig. 3b), while the surrounding lowlands show horizontal arc-perpendicular shortening rates ca.(1-2) x 10 -9 yr -1 (Jura, Rhone Valley, Fig. 3b) (Masson et al., 2019;Nguyen et al., 2016;Sánchez et al., 2018;Walpersdorf et al., 2018).The southern Western Alps are associated with a mix of extension and strike-slip rates ca.(1-2) x 10 -9 yr -1 .These strain rates amplitudes are at the limit of resolution of GNSS techniques and can only be identified through spatial low-pass filtering to remove short-wavelength noise.In a detailed analysis of real and synthetic data using a spatial Gaussian filter, (Masson et al., 2019) show that a filter half-width of 70-100 km provides the best combination of spatial resolution and noise reduction for data in France and the Western Alps.This filter allows the identification of horizontal strain rate signals with a spatial coherence of ca.100-200 km and a formal 95% confidence interval ca.(0.2-1.0) x 10 -9 yr -1 , depending on the area, data quality, and network geometry.Considering the dimensions of the Western Alps, their icecap during the LGM, and the lithosphere elastic thickness (cf.Section 3), 100-200 km is a reasonable estimation of the expected wavelength of the GIA signal.Thus, in the following, we apply the same 90-km-half-width Gaussian filter on GNSS velocities and GIA model velocities for strain rate comparisons.Details of the method and additional strain rate maps at different filtering half-widths can be found in the Supplementary Material 1 for reference.This overall deformation pattern (slow horizontal extension rates coupled with relatively fast uplift rates) is at the core of the current debate on the processes responsible for the ongoing geodynamics of the Western Alps.Because the role of regional plate tectonics is very small or null, recent studies consider alternative processes including mantle and slab-tear dynamics (Sternai et al., 2019), surface erosion (Champagnac et al., 2007;Vernant et al., 2013), or GIA (Chéry et al., 2016;Mey et al., 2016).Numerical modeling shows that the isostatic response to erosion can generate the observed extension strain rate pattern, but the associated uplift rates are significantly smaller than the GNSS velocities (Sternai et al., 2019;Vernant et al., 2013).On the contrary, GIA from the Last Glacial Maximum can explain the present-day uplift rates (Chéry et al., 2016;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 vertical and horizontal deformation rates of the Western Alps.

Model setup
Due to the time scale of glacial loading / unloading cycles (10 2 -10 5 yr), GIA mechanical models consider both elastic and viscous deformation of the Earth's crust and mantle, with a large array of Earth rheology hypotheses and simplifications.
In its simplest form, the surface response to GIA can be modeled as that of a thin elastic plate overlying a linear viscous fluid (e.g., Mey et al., 2016).The more complex and most common modeling approach involves a 1D rheology stratification assuming an upper elastic layer, commonly 60-100 km thick, overlying multiple mantle layers with different Maxwell viscoelastic properties (cf.Peltier and Andrews, 1976;Spada et al., 2011).Higher levels of complexities can include asymmetrical 3D variations of rheology (elastic thickness and viscous properties, (e.g., Bagge et al., 2021;Steffen et al., 2006;Wu, 2005), inclusion of viscous bodies within the elastic layer (e.g., Klemann and Wolf, 1999;Wu and Mazzotti, 2007), and non-linear or transient viscous rheologies (e.g., Giunchi and Spada, 2000;Lau et al., 2021;van der Wal et al., 2013).
Owing to the relative simplicity of both its parameterization and computation requirements, the 1D, symmetrical, Maxwell, layered model (hereafter "1D Maxwell") is the standard approach in most GIA studies.However, studies testing more complex rheology structures have shown that their predictions can deviate substantially from those of the "1D Maxwell" model, with better fit to the surface observables in several cases (cf.Bagge et al., 2021;Lau et al., 2021;Steffen et al., 2006;van der Wal et al., 2013).Of particular note is the integration of non-linear transient mantle rheologies, which are also proposed to explain post-seismic relaxation following major earthquakes (e.g., Freed et al., 2010;Qiu et al., 2018).The integration of such time-dependent viscous rheologies in GIA models would imply significant biases in the standard "1D Maxwell" viscosity estimations (Lau et al., 2021), but would also challenge the precept of stress advection / migration (Steffen et al., 2015;Wu, 1992), which directly results from the Maxwell body approximation for the Earth's mantle.
In this study, we consider the simplest form of GIA model predictions and use a thin elastic plate overlying a viscous substratum.The elastic plate flexure (w) in response to the ice load is computed using the Kirchhoff-Love thin-plate analytic solution (Turcotte and Schubert, 2002;Wickert, 2016), with a flexural rigidity defined by the plate thickness (he) and its elastic parameters (Young modulus E = 10 11 Pa, Poisson coefficient n = 0.25).The elastic plate is laterally infinite.The flexure is computed from a superposition of analytical solutions (using gFlex code, Wickert, 2016).Bending strains and stresses in the Using this simplified approach is justified by several points: (1) We are only interested in the present-day state of stress and deformation, to compare with seismicity and geodetic data.Time variations of GIA responses since the LGM are not relevant to our study, as long as they reach similar present-day states.
(2) Alpine GIA models can only be constrained using present-day velocities and strain rates.Thus, the path leading to this present-day state (and thus the dependency on rheology assumptions) are unconstrained.
(3) The Western Alps icecap during the LGM has characteristic dimensions of a few 100s km (Fig. 3).This relatively small load dimension only excites upper mantle relaxation and is thus insensitive to deep mantle rheologies (Steffen et al., 2015).( 4) Small (few kyr) time variations in the deglaciation history only impact our results at present-day by a few percent owing to the mantle relaxation time (cf.Section 3.2).These justifications for our simple approach are illustrated in the Supplemental Material 2, in which we compare its results with those a standard "1D Maxwell" model.

GIA velocities and strain rates
The analytical model allows testing the sensitivity of predicted present-day GIA strain rates and stresses to the large range of assumed elastic plate rigidity and mantle relaxation time.Uplift rate predictions for models with a plate thickness ( 5≤ he ≤ 100 km) and a relaxation time (2000 ≤ tr ≤ 20000 yr) are compared with GNSS uplift rates (Fig. Supp.Mat.5), assuming that the latter are largely due to GIA (Mey et al., 2016;Nocquet et al., 2016;Sternai et al., 2019).Using F-test statistics on these comparisons, we estimate a best-fit set of (he, tr) values to he = 10-20 km and tr = 4500-5500 yr.However, because present-day uplift rates may be impacted by other processes, in the following analysis we allow for a slightly larger range of values: 5 ≤ he ≤ 40 km and 3000 ≤ tr ≤ 7000 yr (in the following, computations and figures are presented for 72 GIA models with (he, tr) values distributed evenly within these ranges).These parameter values indicate a fast and short-wavelength response to the Last Glacial Maximum, in agreement with the results of more complex models (Chéry et al., 2016;Mey et al., 2016).
All retained models predict similar surface horizontal strain rate patterns, with N-S to NNW-SSE extension rates in the inner Western Alps (Fig. 3 and Supp.Mat. 6).These present-day extension rates are maximum in southwestern Switzerland and in the northeastern French Alps, reaching ca.(1-2.5)x 10 -9 yr -1 .The forebulge regions (Rhone Valley, southern French Alps, eastern France) are associated with smaller shortening rates ca.(0.2-0.6) x 10 -9 yr -1 , with a shortening direction perpendicular to the Alpine arc.Transitional areas in between these extension and shortening domains show a variety of extension, strike-slip and shortening rates whose orientations and magnitudes vary strongly with the distance relative to the Würm icecap and with the model parameters (Fig. 3 and Supp.Mat. 6).Thus, expected GIA strain rates in this narrow intermediate region are poorly defined.
To a first order, the GIA model strain rates are consistent with the horizontal strain rates derived from GNSS data in the inner parts of the northern Western Alps (Fig. 3b vs. 3c).In the Swiss-French border region, GIA model and GNSS strain rates agree in both orientations (with a small rotation from N-S to NNW-SSE) and in amplitudes within their uncertainties (GIA model rates ca.(0.9-2.1) x 10 -9 yr -1 vs. GNSS rates ca.(2.5 ± 0.5) x 10 -9 yr -1 ).Similarly, in the French Alps GIA model rates are (0.4-0.8) x 10 -9 yr -1 , in agreement with the GNSS rate ca.(0.9 ± 0.4) x 10 -9 yr -1 .We estimate that GIA accounts for a minimum of 30% and up to 100% of the GNSS extension rates observed in the inner Western Alps.
In contrast, the GIA models do not reproduce the relatively high GNSS E-W extension and strike-slip rates ca.(1.5 ± 0.5) x 10 -9 yr -1 near the southern French-Italian border and in the western Po Plain (Fig. 3b vs. 3c).In this area, located directly south of the Adria / Eurasia rotation pole, the GNSS strain rates may be in part related to Adria micro-plate kinematics that would predict 0.3-0.5 mm.yr -1 of E-W extension in southern Western Alps.Similarly, the GIA models do not explain the GNSS strain rates observed in the French and Swiss Jura and the southern Upper Rhine Graben.There, GIA model strain rates are consistent in style (shortening rates) and orientations with GNSS data, but their amplitudes are about 3-4 time smaller (Fig. 3b vs. 3c).This difference is too large to be explained by variations in the GIA model parameters (cf.Supp.Mat.6), suggesting that the GNSS shortening rates cannot be solely attributed to Alpine GIA and that additional processes must contribute, such as possibly GIA effects from the Fennoscandian LGM ice sheets (e.g., Nocquet et al., 2005) or the deformation due to a mantle plume beneath the Eifel volcanic area (Kreemer et al., 2020).

GIA model vs. fault and seismicity depths
One of the main issues in using GIA models for studying stress perturbations on crustal faults is how the model setup relates to the actual Earth rheology and in particular its brittle / elastic domains.Although the upper elastic layer in GIA models is only a mechanical proxy for the crust and lithospheric mantle flexural rigidity, all studies make the simple assumptions that (1) this elastic layer corresponds to the crust or lithosphere, depending on its thickness, (2) its top corresponds to the Earth surface, and (3) seismogenic faults are embedded in the elastic layer starting at its top (e.g., Hetzel and Hampel, 2005;Steffen et al., 2014).More complex models integrating elastic, brittle and ductile behaviors show that, under long-term geological loads, the lithosphere flexural rigidity resides where its resistance to failure exceeds ca.10-20 MPa, which, in the upper crust, corresponds to depths below ca.1-4 km (Burov and Diament, 1995;Hyndman et al., 2009;Tesauro et al., 2009).To our knowledge, no similar study exists for shorter-term loads like glaciation cycles.Thus, matching "real Earth" depths, i.e., those of earthquakes or faults, with GIA model depths remains an issue.
In the Western Alps, most of the seismicity concentrates in the upper 10 km of the crust, except for the deeper cluster of extension seismicity ca.15-20 km depth east of the French-Italian border (cf.Section 2; Delacou et al., 2004;Mathey et al., 2021).In order to match our GIA model stress predictions to fault and seismicity depths, we consider the following points: (1) The range of GIA elastic plate thickness derived from GNSS data (he = 10-20 km, cf.Section 3) corresponds primarily to the flexural rigidity of the crust.
(2) The upper half of the elastic plate corresponds to a "real Earth" depth range starting ca.1-4 km depth (cf.discussion above) and extending to ca. 6-14 km depth (starting depth plus half the plate thickness).
(3) The lower half of the elastic plate extends from ca. 6-14 km to ca. 11-24 km "real Earth" depths.Thus, to a first order, the shallow (resp. deep) seismicity level can be matched with the upper (resp.lower) half of the GIA elastic plate.In the following, we use the maximum GIA stress perturbation derived at the top (resp.bottom) of the elastic plate to discuss the impact on fault rupture, keeping in mind that the closer to the plate center, the smaller the GIA-induced stress perturbation.
The impact of these stress perturbations on a given fault can be estimated using the variation of Coulomb Failure Stress (DCFS) derived from the fault geometry and the associated fault shear and normal stress (King et al., 1994): with µ' the fault effective friction coefficient, and Dt and Dsn the shear and normal stress perturbations on the fault.The faultslip style associated with the stress perturbation is given by the rake (r) (Bott, 1959): with td and ts the shear stresses in the fault dip and azimuth directions, respectively.
In a first step, we consider the Coulomb Failure Stress perturbation to discuss the tendency of present-day GIA effects to promote slip (DCFS > 0) or inhibit slip (DCFS < 0) on a subset of fault systems in the Western Alps (a full Coulomb Failure Stress analysis integrating region stress regimes is discussed in Section 4.3).We consider three fault structures directly below the former icecap (Fig. 4a): the Vuache Fault, the Belledonne Fault, and a series of unknown faults associated with normalfaulting earthquakes in the Briançonnais directly east of the Penninic Front (hereafter IMNF, Inner Massifs Normal Faults).
The fault geometries are from neotectonic databases (Grellet et al., 1993;Jomard et al., 2017) or assuming Andersonian geometries for the unknown structures.Due to the fault dip uncertainties, we test a large range of dip angles (25-90˚).
Seismicity associated with these faults concentrates in the upper crust ca.5-10 km depth (Mathey et al., 2021;Thouvenot et al., 2003a), associated with bending stresses in the upper half of the elastic plate (cf.Section 4.1).
Because the present-day GIA stress perturbations correspond to horizontal compressive stresses, faults with nearvertical dip angles (80-90˚) and standard effective friction coefficient (µ' = 0.6) tend to be clamped and are associated with negative DCFS values (i.e., slip inhibition).In contrast, faults with low effective friction (µ' = 0.1) or medium to low dip angles mostly show positive DCFS between 0.1 and 2 MPa, indicating a tendency of GIA to promote fault slip in these cases.However, these models are associated with fault-slip rakes corresponding to mostly reverse faulting (r = 50-130˚, Fig. 4b), opposite to that of earthquake normal-faulting mechanisms along the IMNF (r » -90˚, Fig. 4b, Mathey et al., 2021) and at odd with the strike-slip mechanisms on the Belledonne and Vuache Faults (r » 180˚ and r » 0˚, respectively, Fig. 4b, Rabin et al., 2018;Thouvenot et al., 2003).Thus, in all cases, the present-day GIA stress perturbations do not favor the observed seismicity on the main fault structures of the northern French Alps below the former Würm icecap.
The impacts of GIA along and outside the Würm icecap margin are more difficult to assess for two main reasons: present-day GIA stress perturbations are much smaller and, mostly, they are strongly sensitive to the model parametrization (elastic plate thickness, rheology assumptions).As an example, we estimate the impact on the Moyenne Durance Fault, one of the best characterized active faults in the southern French Alps (Cushing et al., 2008).There, present-day GIA stress perturbations correspond to horizontal NE-SW tension ranging from ca. 0.3-0.7 MPa at the northern end of the fault to ca.
0.1-0.3MPa at its southern end (Fig. 5a).DCFS are systematically positive (for a low effective friction) but, due to bends in the fault segments, the range of predicted rakes is very large: along the northern section, GIA favors right-lateral strike-slip, while favored rakes flip between right-and left-lateral along the central and southern sections (Fig. 5b).Focal mechanisms along the Moyenne Durance Fault correspond to an overall left-lateral kinematics (Cushing et al., 2008), indicating that the present-day GIA stress perturbations may favor or inhibit slip depending on the considered segment.

Present-day full Coulomb Failure Stress
The combination of GIA stress perturbations with the background stress field defines the full Coulomb Failure Stress (CFS) for any given fault: with t and sn the shear and normal stresses acting on the fault, with the associated fault-slip rake given by eq. ( 3).Depending on the geometric relationships between the stress perturbation, the background stress tensor, and the fault geometry, the CFS indicates that a fault should (CFS > 0) or should not (CFS < 0) slip in the local deformation regime.A major difficulty in computing CFS is the definition of the background stress tensor.Given the lack of direct measurements at seismicity depths (5-20 km), this requires several assumptions that can strongly influence the CFS estimation.A second difficulty in the case of the Western Alps is the strong variations in the deformation and stress regimes over distances of a few 10s km (Mathey et al., 2021;Delacou et al., 2004).
In order to test the conclusions drawn from the DCFS analysis (Section 4.2), we compute the full CFS on two major fault systems representative of the central Western Alps present-day seismicity: the northern Briançonnais normal-faulting system (IMNF) and the strike-slip Belledonne Fault (Fig. 2b).In each case, we construct a local background stress tensor using information from the focal mechanism stress inversions of (Mathey et al., 2021) and several additional assumptions for the IMNF and Belledonne systems: • the deformation regime is pure normal faulting (IMNF) or pure strike-slip faulting (Belledonne); • the stress tensor is Andersonian, with the vertical stress (sV, lithostatic pressure) corresponding to s1 (IMNF) or s2 (Belledonne); • the azimuth of the maximum horizontal stress sH = s2 (IMNF) or sH = s1 (Belledonne) is given by the focal mechanism stress inversion; • the differential stress (s1 -s3) is controlled by optimally oriented faults with an effective coefficient of friction µ' (eq. 4) and following the Mohr-Coulomb friction law for a normal (IMNF) or strike-slip (Belledonne) case (Jaeger and Cook, 1979); • the shape of the stress tensor (R = (s1 -s2) / (s1 -s3)) is set to R = 0.5 (i.e., s2 = 1/2 (s1 + s3)); • the same coefficient of friction (µ') is used for the local background stress and the tested faults.
Results of the full CFS computation (local background stress + GIA perturbation) at 5 km depth for the northern Briançonnais extension system are summarized in Figure 6.Assuming that the seismicity is associated with the reactivation of the Penninic Frontal Thrust (Bilau et al., 2021;Sue and Tricart, 1999), the variations of obliquity of the stress tensor to the fault geometry results in negative CFS from ca. -50 MPa to -10 MPa for a standard friction µ' = 0.6 (Fig. 6b).A small friction µ' = 0.1 results in smaller negative CFS (ca.-10 to 0 MPa).Relative to the modeled stress field, the northeastern section of the fault system is more favorably oriented than southern section and shows the largest (less negative) CFS, up to CFS » 0 MPa in the most favorable case (µ' = 0.1, fault dip of 60°).In contrast, if we assume that the seismicity is associated with faults corresponding to the average focal mechanism (azimuth N050, dip 50°; Mathey et al., 2021), this fault geometry is close to optimal orientation relative to the stress tensor and yields small, mostly negative CFS (ca.-5 to +1 MPa, Fig. 6c).In a few specific cases, the GIA stress perturbation pushes the fault to small positive CFS < 1 MPa.These cases correspond to models with a very thin elastic plate (he £ 10 km) and a small friction (µ' = 0.1) on the south-western fault section.
In nearly all cases, the GIA stress perturbations tend to diminish the CFS (render more negative) by a few MPa, which confirms the tendency of present-day GIA to inhibit the normal-faulting seismicity in the Briançonnais region of the Western Alps.This is easily explained by considering that the GIA stress perturbations correspond to horizontal compressive stresses that tend to inhibit normal faulting (Section 4.2).A few specific cases can result in a positive effect of the GIA perturbation on CFS, if the fault is nearly optimally oriented and the GIA stress perturbation correspond to a small subset of the tested parameter range.The case of the Belledonne Fault system shows more variability (Fig. 7).The local background stress is associated with a sH azimuth of N107°, strongly oblique to the strike of the Belledonne Fault.Thus, a standard friction µ' = 0.6 yields large negative CFS (ca.-80 to -60 MPa at 5 km depth, Fig. 7b), whereas a small friction µ' = 0.1 brings the fault -stress configuration closer to an optimal geometry with smaller negative CFS (ca.-10 to -8 MPa, Fig. 7b).The impact of the GIA stress perturbations depends on the assumed fault friction, dip and position within the local stress field: For a standard friction µ' = 0.6, the GIA perturbation systematically reduces the CFS, thus further inhibiting fault slip; for a small friction µ' = 0.1, the GIA perturbation can increase (render less negative) the CFS by 0.1-0.5 MPa if the fault dip is near vertical.
These two examples illustrate the complexity of full CFS estimations and their dependency on the fault and background stress parameters.Yet, their results generally confirm the main conclusion drawn from the DCFS analysis (Section 4.2): the present-day GIA stress perturbation tends to inhibit fault slip for the extension and strike-slip systems of the central Western Alps (below the former icecap).It is worth noting that because the GIA stress perturbations decrease exponentially with time, the crust will eventually return to its unperturbed (pre GIA) normal/strike-slip stress regime, which should translate in a slow increase of seismicity in the Western Alps.A few specific configurations of background stress, fault geometry, and fault friction can result in cases in which the GIA perturbation promotes fault slip.These particular configurations only apply to a very limited set of fault and earthquakes.
The present analysis was carried out for a depth of 5 km, assuming that the GIA perturbation is maximum (i.e., top of the elastic plate, Section 4.1).A smaller GIA perturbation (i.e., deeper in the upper half of the elastic plate) and a deeper analysis depth would simply result in a smaller impact of the GIA on the full stress field and full CFS computations.It is also worth noting the case of the deeper extension seismicity cluster (10-20 km depth) east of the Briançonnais cluster (cf.Section 2, Fig. 2b), for which the GIA stress perturbation would be opposite to that at shallower depths (horizontal tension, Section 4.1) and would likely favor the normal-faulting seismicity.

Discussion -Why GNSS deformation rates should not be directly compared to seismicity and fault slip rates in the Western Alps
As shown by several studies in the last decade, GNSS velocities and strain rates in the Western Alps are characterized, to a first order, by uplift rates ca.1-2 mm yr -1 and radial extension rates ca.(1-3) x 10 -9 yr -1 (Masson et al., 2019;Nguyen et al., 2016;Nocquet et al., 2016;Sánchez et al., 2018;Walpersdorf et al., 2018).These deformation rates are compatible with the GIA effects from the Last Glacial Maximum, provided the Western Alps region behave mechanically as a thin elastic plate over a low-viscosity upper mantle (cf.Section 3, Chéry et al., 2016;Mey et al., 2016).Although the exact contribution of GIA relative to other processes (e.g., slab breakoff, erosion, tectonics) in the total deformation rates is unclear, the good first-order agreement between GIA model predictions and the observed GNSS rates indicates that the former contributes to 30-80 % of the latter, at least in the region below the former Würm icecap (cf.Section 3, Sternai et al., 2019).
The increasing resolution of GNSS horizontal velocities and strain rates has led to direct comparisons of geodetic styles and rates with those of earthquake focal mechanisms.In most cases, these comparisons address the existence or absence of significant aseismic deformation rates captured by GNSS data, with direct implications on the geodynamics of the Western Alps but also eventually on the understanding of seismicity and the characterization of the associated seismic hazard.
Yet, our study shows that, for most faults in the Western Alps, the present-day GIA stress perturbations actually tend either to inhibit fault slip or to promote fault slip with the wrong mechanism compared to the seismicity deformation style.
Thus, at present-day in the Western Alps, GIA from the Last Glacial Maximum does not drive nor promote the observed seismicity, which must be driven by other processes.
This apparent paradox between "GNSS-GIA strain-rate agreement" and "seismicity-GIA stress disagreement" is easily resolved by considering that observed GNSS extension rates are merely a diminution over time of the finite shortening induced in the upper crust by its downward flexure under the Würm icecap.The gradual diminution of the finite shortening (i.e., extensional strain rate) corresponds to the transient return to the pre-ice situation controlled by the upper mantle relaxation time.This situation bears similarities to that observed in forearcs of subduction zones where finite stress and deformation styles can differ drastically from the transient strain rate patterns observed with GNSS data, the latter being due to the interseismic locking of the subduction fault (Mazzotti et al., 2002;Wang, 2000).A similar situation exists for the earlypostglacial reverse-faulting earthquakes in Fennoscandia, which occurred during a period of GIA extension strain rates, and hence must be the expression of compressive stress stored in the lithosphere over a long (tectonic) time (Craig et al., 2016;Muir-Wood, 2000).
This GIA strain rate / stress antinomy has a major implication for not using GNSS velocities and strain rates in direct comparisons with seismicity styles and rates.Because GIA does not promote the current seismicity, its surface expression, captured in GNSS velocities and strain rates, does not have a simple relationship with seismic moments rates or fault slip rates.
In other words, GNSS velocities and strain rates comprise a significant part of transient GIA deformation that does not directly contribute to the observed seismicity, but rather modulates the expression of the mechanisms (e.g., erosion, slab tear, tectonics) driving the current Western Alps geodynamics and seismicity.In this situation, GNSS data cannot inform models based on a steady-state tectonic process (e.g., far-field fault loading), but they can provide important constraints for models combining long-term forcing and GIA transient to study fault slip and seismicity variations during and after glaciations (e.g., (Steffen et al., 2014;Hetzel and Hampel, 2005)).This requires specific models integrating the complexities of the Alpine icecap history, of the regional crust and mantle rheology heterogeneities, and of the local fault characteristics.

Conclusion
Our analysis shows that GIA from the Würm icecap explains ca.30-80 % of the present-day uplift rates and horizontal extension rates observed by GNSS in the Western Alps, in agreement with previous studies (Chéry et al., 2016;Mey et al., 2016;Sternai et al., 2019).Yet, we also show that, for most of the major fault systems, present-day GIA stress perturbations either inhibit fault slip or promote it with the wrong mechanism (compared to the seismicity deformation style).Thus, at present-day in the Western Alps, GIA explains most of the geodetic rates but does not drive nor promote the observed seismicity.This apparent strain rate -stress paradox, resulting from the GIA relaxation and diminution over time of the upper crust flexure, has strong implications for seismicity and seismic hazard studies.First, because GNSS velocities and strain rates are mostly comprised of transient GIA deformation that tends to inhibit seismicity, they cannot be directly used to estimate seismic moment rates or fault slip rates (as would be done in standard "plate-boundary" interseismic studies).Second, because GIA stress perturbations are decreasing exponentially with a time constant ca.5000 yr, their inhibition effects on seismicity are slowly disappearing, which should result in an increase of seismicity over the next thousands of years.Thus, seismic hazard assessments in the Western Alps must integrate the complexities of Würm GIA, constrained through GNSS data and specific modeling, especially for long earthquake return periods and critical facilities.

Figure 1 :
Figure 1: Conceptual model of Glacial Isostatic Adjustment (GIA) induced stresses, strains and strain rates

Figure 3 .
Figure 3.Western Alps geodetic and GIA deformation rates.(a) Vertical velocity field from GNSS permanent networks (Masson et al., 2019).(b) Horizontal strain rate field extracted from Gaussian filter method (cf.text).(c) Best-fit GIA model strain rate field (he = 20 km, tr = 5500 yr).Strain rates are color-coded according to the deformation style.Grey line is Last Glacial Maximum ice extension (Mey et al., 2016).
elastic plate are computed from the second spatial derivative of the flexural response based on the Kirchhoff-Love theory (no vertical shear stress, only in-plane compression / tension stress).The time dependence of the response due to the upper mantle viscous relaxation is computed assuming an exponential decay controlled by a relaxation time (tr): (, , ) = (, )  +, -. / (1) For the ice model, we use the icecap reconstruction at LGM from Mey et al. (2016) based on an ice-flow model constrained by geomorphological ice extent and thickness data.The ice model covers the whole Alps, with a mean (resp.maximum) thickness at LGM ca.21 kyr BP of 415 m (resp.2445 m).In the Western Alps, the main ice load is located in the major valleys, up to 1000-1500 m ice thickness.Additional smaller ice lobes (few 100 m) were located in the Jura.Keeping in line with the simple modeling approach, we assume (1) that the Alpine icecap reached isostatic equilibrium at LGM, owing to the relatively short relaxation time of the viscous mantle (ca.5000 yr, cf.Section 3.2) and (2) an instantaneous deglaciation of the icecap at 15 kyr BP, considering that 80% has melted between 21 and 18 kyr BP and the remaining has disappeared by 12 kyr (Ivy-Ochs et al., 2008).