Fast uplift in the Southern Patagonian Andes due to long and 1 short term deglaciation and the asthenospheric window 2 underneath

. An asthenospheric window underneath much of the South American continent 14 increases the heat flow in the Southern Patagonian Andes, where glacial-interglacial cycles 15 drive the building and melting of the Patagonian Icefields since the latest Miocene. The Last 16 Glacial Maximum (LGM) was reached ~26000 years Before Present (BP). Significant 17 deglaciation onsets between 21000 and 17000 years BP, subject to an acceleration since the 18 Little Ice Age (LIA), ~400 years BP. Fast uplift rates of up to 41 ± 3 mm/yr are measured by 19 GNSS around the Southern Patagonian Icefield and currently ascribed to post-LIA 20 lithospheric rebound, but the possible longer-term post-LGM rebound is poorly constrained. 21 These uplift rates, in addition, are one order of magnitude higher than those measured on 22 other glaciated orogens (e.g., the European Alps), which raises questions about the role of the


Introduction
Vertical displacements of the Earth's surface with respect to the geoid occur in response to the motion of crustal and mantle rock masses due to mantle dynamics, plate tectonics, and the redistribution of sediments, water, and ice by surface processes (e.g.Molnar and England, 1990;Watts, 2001;Champagnac et al., 2012;Sternai, 2023;Cloetingh et al., 2023).For instance, an excess of topography in orogenic regions due to convergence, crustal shortening, and magmatism deflects the lithosphere downward, whereas surface unloading by erosion and ice melting causes upward deflection of the lithosphere, known as "isostatic" adjustment (e.g.Peltier and Andrews, 1976;Peltier, 1996Peltier, , 2004;;Mitrovica and Forte, 1997;Butler and Peltier, 2000;Kaufman and Lambeck, 2002;Watts, 2001;Turcotte and Schubert, 2002;Sternai et al., 2016a).Glacial isostatic adjustment (GIA) models study the viscoelastic response of the solid Earth to the building and melting of regional ice sheets, and they commonly use GNSS and/or satellite-measured rock uplift rates in regions subject to deglaciation to estimate, for instance, the regional mantle rheology and sea-level changes (e.g.Peltier and Andrews, 1976;Peltier, 1996Peltier, , 2004;;Mitrovica and Forte, 1997;Kaufman and Lambeck, 2002;Stuhne and Peltier, Published by Copernicus Publications on behalf of the European Geosciences Union.Wal et al., 2015;Peltier et al., 2018;Whitehouse, 2018).Most of the GIA studies address the post-LGM (Last Glacial Maximum) around 21 000 yr BP (years before present) deglaciation as a trigger for increasing uplift rates in glaciated regions (e.g.Peltier, 2004;Whitehouse, 2018).The magnitude of uplift rates is set primarily by the lithosphere and asthenosphere viscosities, which depend, amongst other factors, on the thermal field at depth (McKenzie and Richter, 1981;McKenzie and Bickle, 1988;Gurnis, 1989;Ranalli, 1995Ranalli, , 1997;;Kaufman et al., 1997;Watts, 2001;Turcotte and Schubert, 2002).While the GIA theory is well developed, few studies use thermo-mechanical visco-elasto-plastic (non-Newtonian Earth layers) geodynamic models to estimate uplift rates in response to surface load changes to be compared with GNSS data.Here, we use this approach to constrain the role of the solid Earth rheology in setting the rates of surface vertical displacement in southern Patagonia, which hosts the biggest continental ice sheets outside of Antarctica and presents ongoing very high rock uplift rates (Ivins and James, 2004;Dietrich et al., 2010;Lange et al., 2014;Richter et al., 2016;Lenzano et al., 2023).
GNSS-measured data show ongoing vertical rock uplift rates between 18 ± 3 and 41 ± 3 mm yr −1 in the northern part (18-50.5°S) of the SPI (Fig. 1b), decreasing to val-ues between 2 ± 6 and 17 ± 5 mm yr −1 in its southern part (50.5-51.5°S)(Ivins and James, 2004;Dietrich et al., 2010;Lange et al., 2014;Richter et al., 2016;Lenzano et al., 2023).Such outstandingly high uplift rates, especially in the northern part of the SPI, are currently ascribed to the lithospheric viscoelastic GIA following the LIA, which was responsible for an ice loss of 503 ± 101.1 km 3 in the SPI (Glasser et al., 2011).To match the very high observed uplift rate budget, previous GIA studies infer a low asthenosphere viscosity (in the order of 10 18 Pa s) and a thin elastic lithosphere (∼ 35 km thick) (Ivins andJames, 1999, 2004;Klemann et al., 2007;Dietrich et al., 2010;Lange et al., 2014;Richter et al., 2016;Ávila and Dávila, 2020;Mark et al., 2022;Lenzano et al., 2023).Although this is consistent with abnormally high asthenospheric mantle temperatures, viscosity estimates from these previous studies are untied to the regional thermal regime, which prevents a more thorough characterisation of the role of the asthenospheric window underneath the SPI in affecting the observed uplift rates.In addition, the contribution of post-LGM deglaciation to present-day rock uplift rate was marginally addressed (Ivins andJames, 1999, 2004;Klemann et al., 2007).Here, we perform fully coupled thermo-mechanical numerical geodynamic experiments forced by surface unloading scaled on post-LIA and post-LGM ice melting to evaluate their relative contribution to the observed regional uplift rates.Numerical experiments account for a range of positive thermal anomalies in the asthenosphere to further assess the role of the asthenospheric window in setting the mantle viscosity and associated postglacial rebound.Focusing on the magnitude rather than on the pattern of the inferred surface uplift rates due to limited information on the spatial-temporal variations in the ice net mass balance and thickness since the LGM (e.g.Davies et al., 2020), we use the observed budget of rock uplift rate to constrain plausible thermal and viscosity structures at depth and the timing of post-glacial rebound.

Methodology
As reference, we used the GNSS-derived data from 31 GPS stations installed at 380 km in north-south directions and 130 km in east-west directions around the SPI since 1996, published in Lange et al. (2014).The observed and estimated regional aseismic viscoelastic uplift rates presented in that study are shown in Fig. 1b.Details on the GPS data acquisition and analysis are given in the reference study (Lange et al., 2014).

Numerical model
We use a fully coupled thermo-mechanical, visco-elastoplastic numerical geodynamic model to quantify the effect of thermal anomalies in the asthenospheric mantle on the magnitude of surface uplift rates due to deglaciation.We provide a short overview of the governing equations hereafter, while a detailed description of numerical technique can be found, for instance, in Gerya and Yuen (2007), Gerya (2019), Sternai (2020), Sternai et al. (2021), andMuller et al. (2022).The continuity equation allows for the conservation of mass during the displacement of a geological continuum: where ρ is the local density, t is time, v is the velocity vector, and ∇ is the divergence operator.The momentum equation describes the changes in velocity of an object in the gravity field due to internal and external forces: where σ ij is the stress tensor, x i and x j are spatial coordinates, and g i is the ith component of the gravity vector.The energy equation allows for the conservation of energy during advective and conductive heat transfer in the continuum: where P is pressure; T is temperature; C P is specific heat capacity at a constant P ; c is the thermal conductivity; and H r , H s , H a , and H l are the volumetric heat productions by radiogenic, shear, adiabatic, and latent heat, respectively.
, where σ ij is the deviatoric stress tensor and ε ij (viscous) is the viscous deviatoric strain rate tensor.Ductile deformation is thermally activated generating viscous flow, which is calculated according to the material shear viscosity, η: where σ II is the second invariant of the deviatoric stress tensor, n is the stress exponent, A d is the pre-exponential factor, E a is the activation energy, V a is the activation volume, and R is the gas constant.ε ij (viscous) is computed as where δ ij is the Kronecker delta, εkk is the bulk strain rate in response to irreversible volume changes, and η bulk is the bulk viscosity.Recoverable deformation is defined by the elastic deviatoric strain rate tensor, ε ij (elastic) , as where µ is the shear modulus and Dt is the objective corotational time derivative of the deviatoric stress tensor.The plastic deformation, brittle and localised, occurs at a low temperature when the absolute shear stress limit, σ yield , is reached: where C is cohesion and ϕ is the effective internal friction angle.The plastic strain rate tensor, ε ij (plastic) , is defined as where X is the plastic multiplier which satisfies the plastic yielding condition σ II = σ yield .The bulk strain rate tensor, ε ij (bulk) , integrates the viscous, elastic, and plastic deformation:

Reference model setup and modelling approach
The model domain is 700 km wide and 120 km thick to account for realistic lithosphere and asthenospheric mantle thicknesses of South America at the latitudes of the SPI (Van der Meijde et al., 2013;Ávila andDávila, 2018, 2020).From top to bottom, the model accounts for 10 km of "sticky" air, 30 km of continental crust (with rheology of quartzite, Ranalli, 1995), 30 km of lithospheric mantle, and 50 km of asthenospheric mantle (with rheology of dry dunite, Ranalli, 1995); this is in agreement with literature data (e.g.Van der Meijde et al., 2013;Ávila andDávila, 2018, 2020).The initial geotherm is piecewise linear, resulting from an adiabatic temperature gradient of 0.5 °C km −1 in the asthenosphere (Turcotte and Schubert, 2002) and from thermal boundary conditions equal to 0 °C at the surface and 1327 °C at the bottom of the lithosphere, with a nil horizontal heat flux across the vertical boundaries.The rheologic and thermal structure of the reference model gives a lithospheric elastic thickness, T e (sensu Burov and Diament, 1995), of ∼ 30 km, comparable to previous estimates underneath the SPI based on GIA models (Ivins and James, 1999;Dietrich et al., 2010;Lange et al., 2014), heat flow data (Ávila and Dávila, 2018), waveform inversion (Robertson Maurice et al., 2003), and low-temperature thermochronology data (Thomson et al., 2010;Guillaume et al., 2013;Georgieva et al., 2016Georgieva et al., , 2019;;Stevens Goddard and Fosdick, 2019;Ávila et al., 2023;Muller et al., 2023).Rocks' rheological properties are listed in Table 1.
The numerical model uses the finite differences with the marker-in-cell technique, resolved by 51 × 61 nodes in the horizontal, x, and vertical, y, directions, respectively, distributed on an Eulerian grid that accounts for a maximum resolution of 1 km along the y direction in the upper part of the model domain and ∼ 13 km in the x direction.Along the x and y dimensions, 400 × 400 Lagrangian markers are randomly distributed and used for advecting the material properties (Gerya and Yuen, 2007;Gerya, 2019).The material properties carried by Lagrangian markers are then interpolated onto the Eulerian grid via a fourth-order Runge-Kutta interpolation scheme.An internal free surface is simulated through the 10 km thick layer of sticky air.The velocity boundary conditions are free slip at all boundaries (x = 0 and x = 700 km; y = 0 and y = 120 km).
On the top of the crust and in the middle of the model domain we impose a 2 km thick pseudo-ice cap to simulate lithospheric unloading during deglaciation (Fig. 2a).The pseudo-ice cap has an initial density, ρ ice , of 920 kg m −3 (Harvey et al., 2017) (Table 1), and we compute the surface load through time, L, as where g is the gravity acceleration and h ice is the ice cap thickness.The load change due to the deglaciation occurs by gradually and uniformly reducing h ice in time (Fig. 2b and c).
In the models, the lateral extent of the pseudo-ice cap does not change throughout the deglaciation.Although this simplification may affect the inferred pattern of post-glacial rebound, it greatly facilitates the simulation of deglacial lithospheric unloading without significantly affecting the magnitude of post-glacial rebound, which is the main focus here.All simulations account for some spin-up time before the deglaciation begins so that the lithosphere-asthenosphere system adjusts to the pseudo-ice cap's initial load.The uplift rate during the deglaciation is calculated through time as the surface elevation change which results from the modelled strain field divided by the viscoelastic time step (i.e.U = (z curr − z prev )/t, where z curr is the modelled topography at the considered time step, z prev is the modelled topography at the previous time step, and t is the viscoelastic time step duration).Given the geologically short time window investigated here, we neglect deformation related to longer-term tectonic forces (Breitsprecher and Thorkelson, 2009;Guillaume et al., 2013;Eagles and Scott, 2014;Muller et al., 2021).The parametric study focuses on the asthenospheric mantle potential temperature (sensu McKenzie and Bickle, 1988) which accounts for positive thermal anomalies, TA, of up to 200 °C in steps of 50 °C; this is then added to the reference asthenospheric mantle potential temperature of 1265 °C (McKenzie and Bickle, 1988;Currie and Hyndman, 2006;Ávila andDávila, 2018, 2020;Sternai, 2020;Mark et al., 2022) to mimic the presence of a slab window at depth.

Results
Results are shown in Table 2 and Figs.4-7.In agreement with the theory of lithospheric flexure (e.g.Turcotte and Schubert, 2002), the deglaciation triggers uplift in the region covered by the melting pseudo-ice cap and subsidence in the neighbouring regions .Overall, increasing the asthenospheric mantle potential temperature decreases the asthenospheric viscosity, with significant effects on the magnitude of the modelled surface velocity field.The asthenosphere viscosity ranges between 10 22 -10 19 Pa s in simulations with TA equal to 0 (reference model), 50, and 100 °C and between 10 19 -10 16 Pa s in simulations with TA equal to 150 and 200 °C (Fig. 3a-d).Lithospheric warming due to increasing asthenospheric mantle potential temperature also leads to a reduction in the lower lithosphere viscosity (from 10 22 to 10 20 Pa s), thereby decreasing the integrated lithospheric strength.
In Model set 1 for post-LGM deglaciation, when TA equals 0 (reference model) the maximum uplift rate is < 1 mm yr −1 during the first 5000 years of the deglaciation, increasing gradually up to 9.5 mm yr −1 in the later stages of the deglaciation (i.e.20 000 years, Fig. 4).When TA equals 50, 100, 150, and 200 °C, the maximum uplift rates can reach up to ∼ 2, ∼ 5, ∼ 12, and ∼ 15 mm yr −1 , respectively, already in the first 1000 years of the deglaciation (Fig. 4a).When TA is 50 and 100 °C, the maximum uplift rate is subject to a protracted increase in time, reaching up to ∼ 12 and ∼ 14 mm yr −1 after 20 000 years of deglaciation (Figs.4b-d  and 7a).For TA equal to 150 and 200 °C, the maximum uplift rate reaches a plateau between 11 and 17 mm yr −1 during the 20 000 years of deglaciation (Figs. 4 and 7a, Table 2a).After the end of the deglaciation, the maximum uplift rate takes longer than about 5000 years to re-equilibrate to 0 mm yr −1 when TA ≤ 100 °C, whereas it drops to 0 mm yr −1 almost immediately when TA is 150 or 200 °C (Fig. 7a).
In Model set 2 for post-LGM deglaciation, the maximum uplift rate is less than 2 mm yr −1 during the first 1000 years of deglaciation when TA is 0, 50, and 100 °C, whereas it reaches up to ∼ 22 and ∼ 30 mm yr −1 during the first 1000 years of deglaciation when TA is 150 and 200 °C (Figs.5a and 7b, Table 2).Between 5000 and 10 000 years of deglaciation, the maximum uplift rate increases to ∼ 19, ∼ 25, and ∼ 36 mm yr −1 , respectively, when TA is 0, 50, and 100 °C, whereas it reaches up to between 36 and 41 mm yr −1 between 50 000 and 1000 years of deglaciation when TA equal to 150 and 200 °C.The maximum uplift rate decreases slower if TA is 0, 50, and 100 °C, taking longer than 5000 year after the deglaciation to drop to values < 5 mm yr −1 (Fig. 7b and Table 2b), whereas it quickly drops to < 2 mm yr −1 when the deglaciation is over and TA is 150 and 200 °C (Figs.5b-d and 7b).Overall, a warmer and less viscous asthenosphere generates a higher magnitude and a fast-changing post-glacial rebound than a cooler and more viscous asthenosphere.
In the post-LIA model set, the maximum uplift rate is ∼ 1.4, ∼ 2.3, and ∼ 2.2 mm yr −1 during the first 100 years of deglaciation when TA is respectively 0, 50, and 100 °C, whereas it reaches ∼ 8.3 and ∼ 23 mm yr −1 during the same interval when TA is respectively 150 and 200 °C (Figs.6a  and 7c, Table 2c).Between 200 and 400 years of deglaciation, the maximum uplift rate reaches ∼ 1.9, ∼ 2.5, and ∼ 3 mm yr −1 when TA is equal to 0, 50, and 100 °C and ∼ 14 and ∼ 25.5 mm yr −1 when TA is 150 and 200 °C, respectively (Figs. 6c-d and 7c, Table 2c).When the deglaciation ends, the maximum uplift rate drops to ∼ 0 mm yr −1 in ∼ 100 years when TA ≤ 100 °C, whereas it takes longer than 1000 years when TA equals 150 or 200 °C (Fig. 7c).Overall, a warmer and less viscous asthenosphere generates a highermagnitude post-glacial rebound which, however, takes much longer to re-equilibrate to 0 mm yr −1 after the end of the deglaciation than a cooler and more viscous asthenosphere.

Discussion
Our modelling is simplistic in that we impose a linear and uniform ice loss instead of a more realistic ice sheet melting pattern in space and time (Fig. 2b and c).Although the stratigraphic and geochronologic record is fairly precise for the post-LGM ice extent (e.g.Lagabrielle et al., 2004;Rabassa, 2008;Glasser et al., 2011;Davis and Glasser, 2012;Strelin et al., 2014;Kaplan et al., 2016;Martinod et al., 2016;Bendle et al., 2017;Reynhout et al., 2019;Davies et al., 2020), information about melting velocities and associated ice thickness and redistribution of the surface masses is limited for the time windows investigated here.GNSS, SRTM, and GRACE data, constraining the net ice mass balance only during the last few decades, still showing some discrepancies (e.g.Aniya, 1996;Aniya et al., 1997;Rignot et al., 2003;Ivins andJames, 1999, 2004;Chen et al., 2007;Dietrich et al., 2010;Ivins et al., 2011;Jacob et al., 2012;Lange et al., 2014;Willis et al., 2012;Richter et al., 2016;Gómez et al., 2022;Lenzano et al., 2023).Tracing back the post-LGM or Holocene ice loss rate from current measurements is difficult considering that climate was at least 6 °C colder during the LGM (Hulton et al., 2002;Sugden et al., 2002;Seltzer et al., 2021;Yan et al., 2022).As a result, previous models have assumed simple deglaciation histories as well (e.g.Ivins andJames, 1999, 2004;Hulton et al., 2002;Klemann et al., 2007;Ivins et     2011; Boex et al., 2013).Measurements of regional erosion rates since the LGM are between 0.02 and 0.83 mm yr −1 (Fernandez et al., 2016).However, given the short time intervals investigated here, it seems reasonable to assume that the eroded material is still in the transport zone and therefore does not significantly contribute to unloading the surface of the orogen.If one refers to erosion rates from lowtemperature thermochronology, even though these measures quantify erosion rates over millions of years and not millennia, Fosdick et al. (2013), Herman and Brandon (2015), Fernandez et al. (2016), andMuller et al. (2023) suggest using values between 0.1 and 1 mm yr −1 from 7 to 4 Ma, followed by a period of erosional quiescence (< 0.1 mm yr −1 ) and a possible increase to 1 mm yr −1 in the last ∼ 2 Myr in the SPI region (Muller et al., 2023).Supposing that these erosion rates still apply in the last ∼ 20 000 years, this would translate into 2-20 m of rocks eroded on average since the LGM, leading to local unloading of approximately 60-600 kPa if one assumes a crustal density of 3000 kg m −3 .Such stress change is approximately equivalent to the melting of about 6-60 m of ice, whereas we simulate the melting of 200-1500 m of ice in our simulations.The forcing of global cooling in increasing erosion rates during the Quaternary is, however, debated, and not widely quantified in Patagonia nor worldwide (Valla et al., 2012;Champagnac et al., 2014;Herman et al., 2013Herman et al., , 2018;;Herman and Brandon, 2015;Fernan-dez et al., 2016;Georgieva et al., 2019;Yan et al., 2022).
Even if long-term erosion rates contribute to present-day uplift rates (Herman et al., 2018), since they are comparable to those of the European Alps, for instance, we assume a similar contribution to regional uplift rates (i.e.generally a fraction of a mm yr −1 ; Sternai et al., 2019); this is a negligible contribution in the context of the southern Patagonian Andes.We also assume a homogeneous lithosphere and neglect lateral viscosity variations in the asthenosphere, despite the long-term southern Andean orogenic history (Cande and Leslie, 1986;Ramos, 2005;Breitsprecher and Thorkelson, 2009;Muller et al., 2021) and the suggested contribution from lateral rheological heterogeneities (Klemann et al., 2007;Richter et al., 2016).Overall, notwithstanding these limitations in the model, our fully coupled numerical thermomechanical geodynamic experiments provide realistic uplift rates (Figs.4-7) that one can compare to current geodetic observations.Following the example of previous studies (Ivins andJames, 1999, 2004;Klemann et al., 2007;Dietrich et al., 2010;Lange et al., 2014;Richter et al., 2016;Lenzano et al., 2023), we discuss our results assuming that GNSS-measured rock uplift rates are mostly related to the deglaciation history and only marginally controlled by the longer-term geodynamics (e.g.Ramos, 2005;Breitsprecher and Thorkelson, 2009;Eagles and Scott, 2014;Muller et al., 2021).The elastic thickness of the lithosphere (T e ) varies between the simulations according to the imposed asthenospheric thermal anomaly, but it is generally lower than 30 km, resulting in a decoupled lithospheric rheology (sensu e.g.Burov and Diament, 1995), as shown by the yield stress envelope in Fig. 2a.This results in predominant elastic deformation in the upper crust (below the ∼ 300 °C isotherm) and upper mantle lithosphere (below the ∼ 700 °C isotherm) and viscous deformation in the lower crust and lower lithospheric mantle and asthenosphere (Fig. 3).We remark that, when we impose higher temperatures in the asthenospheric mantle, shallower 300 and 700 °C isotherms decreases T e and increases the isostatic surface uplift rates.Lithospheric thinning due to the asthenospheric window underneath southern Patagonia thus affects the regional uplift rates as previously suggested (Ávila andDávila, 2018, 2020;Mark et al., 2022;Ben-Mansour et al., 2022;and Avila et al., 2023).
The inferred maximum post-LIA uplift rate of up to a few mm yr −1 from experiments with or without a low asthenospheric thermal anomaly (TA ≤ 100 °C, Fig. 7c) is within the same order of magnitude of maximum uplift rates measured in collisional orogens such as the European Alps (Sue et al., 2007;Serpelloni et al., 2013;Walpersdorf et al., 2015;Sternai et al., 2019) and the Himalayas (Larson et al., 1999).Since these collisional orogens are characterised by a thicker lithosphere (Geissler et al., 2010;Ravikumar et al., 2020), they are likely less sensitive to mantle dynamics than the southern Patagonian Andes.When we consider lithospheric unloading due to post-LGM deglaciation of a wider ice sheet, however, the inferred maximum uplift rate via Model set 1 and Model set 2 reaches up to 10 mm yr −1 and 20 mm yr −1 , respectively, even without an asthenospheric thermal anomaly (Fig. 7a and b).This suggests a likely contribution from a long-term post-glacial rebound to the present-day uplift rates measured in the SPI.
In the southern Patagonian Andes, GIA models estimated the regional asthenosphere viscosity to be between 1.6 and 8 × 10 18 Pa s (Ivins andJames, 1999, 2004;Dietrich et al., 2010;Willis et al., 2012;Lange et al., 2014;Richter et al., 2016;Lenzano et al., 2023).Similarly, the asthenosphere viscosity from our models when TA > 100 °C is < 10 19 Pa s, with the lowest viscosity value of 10 16 Pa s imposed where partial melting, supported by the regional Holocene volcanism (Stern and Kilian, 1996) and by geophysical data (e.g.shear wave velocity data by Mark et al., 2022)  vide max uplift rates between 14 and 26 mm yr −1 toward the end of the LIA deglaciation (Fig. 7c).Even with a very low viscosity asthenosphere, the rebound due to short-term post-LIA deglaciation does not reach the presently observed maximum uplift rates of 41 ± 3 mm yr −1 .Experiments that account for a low viscosity asthenosphere and long-term post- LGM deglaciation lasting for 20 000 and 10 000 years reach up to ∼ 25 and ∼ 42 mm yr −1 of an uplift rate during the final stages of the deglaciation (Fig. 7a and b), respectively, which is comparable to present-day values.Results, therefore, indicate that the outstanding observational budget of rock uplift in the SPI is matched only when accounting for higher-thannormal asthenospheric mantle temperatures, thereby highlighting the relevance of the regional asthenospheric window.Consistently, even though the higher heat flow is currently further north from our study region near the CTJ (46-48°S) (Ramos, 2005;Breitsprecher and Thorkelson, 2009;Ávila andDávila, 2018, 2020;Ben-Mansour et al., 2022), increased asthenospheric temperatures beneath southern Patagonia are highly supported by the geophysical data (e.g.Russo et al., 2010Russo et al., , 2022;;Mark et al., 2022;Ávila andDávila, 2018, 2020;Ben-Mansour et al., 2022).
Because of the limited knowledge regarding the timing and amount of ice loss since the LGM (e.g.Ivins andJames, 1999, 2004;Hulton et al., 2002;Klemann et al., 2007;Boex et al., 2013;Davies et al., 2020), it is difficult to position in time present-day uplift rate measurements within the investigated deglaciation scenarios to assess the contribution of post-LGM, post-LIA, and present-day deglaciation to the maximum uplift rate budget.In the faster post-LGM deglaciation scenario (Model set 2) the observed maximum uplift rate budget is attained in about 10 000 years of deglaciation, but only minor residual rebound could be observed today regardless of the amount of ice loss (Fig. 7b).If post-LGM deglaciation occurred more slowly (Model set 1), this event may contribute up to 40 % to the present-day uplift rate budget (Fig. 7a).Although it is difficult to reconcile this scenario with the geomorphological and geochronological evidences (Hulton et al., 2002;Boex et al., 2013;Davis and Glasser, 2012;Martinod et al., 2016;Bendle et al., 2017;Thorndycraft et al., 2019;Davies et al., 2020), it appears that post-LIA rebound alone cannot cover the entire budget of the observed uplift rates even with the highest tested TA, which points to a non-negligible contribution from post-LGM deglaciation.This latter conclusion is reinforced by estimates of the mantle relaxation time, τ r , as (Turcotte and Schubert, 2002) where v is the asthenosphere viscosity, λ is the width of the ice sheet, and g is the gravity acceleration.Using 10 16 < v < 10 18 Pa s and λ = 200 km leads to ∼ 2000 < τ r <∼ 200 000 years, a time range considerably longer than the post-LIA deglaciation and including full Pleistocene glacial-interglacial cycles (Ruddiman et al., 1986).Although an increasingly negative ice mass balance in the last ∼ 50 years has contributed to the elastic lithospheric uplift rates (Dietrich et al., 2010;Lange et al., 2014), a longer-term contribution from the viscous lithosphere is necessary to explain the GNSS-measured uplift rates (Ivins and James, 2004;Dietrich et al., 2010;Lange et al., 2014;Richter et al., 2016;Lenzano et al., 2021).As a final consideration, our models suggest that we shall measure regional uplift rates in the order of the tens of cm yr −1 in the next century if the currently observed ice loss rate of at least −20 Gt yr −1 in the SPI (Willis et al., 2012) will continue until the total meltdown of the ice sheet in ∼ 200 years.

Conclusions
We propose that rock uplift rates of up to 40 mm yr −1 in the southern Andes are due to both post-LIA and long-term post-LGM lithospheric rebound, as postulated for other glaciated orogens (e.g. the European Alps, Fennoscandia, and North America; Peltier et al., 2018).We also propose that currently observed uplift rates in the southern Andes are enhanced by a mantle thermal anomaly of at least 150 °C due to the regional asthenospheric window.Asthenospheric thermal anomalies higher than 200 °C are unlikely and would decrease the asthenospheric viscosities to unrealistic values (less than 10 16 Pa s).Our thermo-mechanical visco-elastoplastic forward modelling approach thus helps to constrain the increase in temperature in geodynamic asthenospheric upwelling contexts, such as in southern Patagonia (Russo et al., 2010(Russo et al., , 2022;;Ávila andDávila, 2018, 2020;Mark et al., 2022;Ben-Mansour et al., 2022).

Figure 1 .
Figure 1.Regional context and uplift rate data.(a) Map of southern Patagonia with the Southern Patagonian Icefield (SPI), Northern Patagonian Icefield (NPI), and the Cordillera Darwin Icefield (CDI) in light blue.The map also shows the approximate extension of the Patagonian Ice Sheet at the Last Glacial Maximum (LGM) (adapted from Thorndycraft et al., 2019) and the approximate extension of the present-day asthenospheric window (dashed region) beneath the South American (SAM) Continent (adapted from Breitsprecher and Thorkelson, 2009).In the Pacific Ocean, the spreading ridges (s.r.; thick-black lines) and transform faults (t.f.; thin-black lines) separate the Nazca (NZ) and the Antarctic (AT) plates.The subduction trench is also highlighted in black.The arrows show the approximate rate and direction of subduction of the oceanic plates (adapted fromDeMets et al., 2010).(b) A zoomed-in view of the SPI with GNSS-measured rock uplift rates (colour-coded disks) used to estimate the viscoelastic uplift rates inLange et al. (2014).
standard densities of solid rocks; Ea is the activation energy; Va is the activation volume; n is the stress exponent; C is cohesion; ϕ eff is the effective internal friction angle; c is thermal conductivity; µ is the shear modulus; Cp is the specific heat capacity; Hr and H l are the radiogenic and latent heat productions, respectively; α is thermal expansion; and β is compressibility.Qz and Ol are quartzite and olivine, respectively.All rheological and partial melting laws/parameters are based on experimental rock mechanics and petrology

Figure 2 .
Figure 2. Reference numerical model setup.(a) Thermo-mechanical numerical model domain with rheological layers (Table1), isotherms (white lines), and yield strength ( σ = σ 1 − σ 3 ) profile (yellow line).The yield strength ( σ ) profile is not scaled and aims to show the proportionality of the yield strength amongst the layers, dependent on the temperature and composition (Eq.4).(b, c) Ice thickness vs. time used in the numerical models to simulate post-LGM deglaciation in two model sets (b) and in post-LIA deglaciation (c).

Figure 3 .
Figure 3. Distribution of viscosity and velocity vectors in the numerical models.(a, c) Reference models without an asthenospheric thermal anomaly, TA = 0 °C, in the last time step of post-LIA deglaciation (a) and of Model set 1 of post-LGM deglaciation (c).(b, d) Model with the higher simulated asthenospheric thermal anomaly, TA = 200 °C, in the last time step of post-LIA deglaciation (b) and of Model set 1 of post-LGM deglaciation.Model set 2 has a very similar viscosity and velocity vectors distribution to Model set 1 in the last deglaciation time step.Velocity vectors do not have the same scaling and are only meant for visualisation purposes.

Figure 4 .
Figure 4. Surface uplift rates vs. distance for Model set 1 of post-LGM deglaciation.Panel (a) shows t = 1000 years of deglaciation, panel (b) shows t = 5000 years of deglaciation, panel (c) shows t = 10 000 years of deglaciation, and panel (d) shows 20 000 years of deglaciation.Different line colours correspond to different TA.

Figure 5 .Figure 6 .
Figure 5. Surface uplift rates vs. distance for Model set 2 of post-LGM deglaciation.Panel (a) shows t = 1000 years of deglaciation, panel (b) shows t = 5000 years of deglaciation, panel (c) shows t = 10 000 years of deglaciation, and panel (d) shows 20 000 years of deglaciation.Different line colours correspond to different TA.

Figure 7 .
Figure 7. Maximum uplift rates vs. time for model sets of deglaciation with different TA.(a) Model set 1 of post-LGM deglaciation accounting for 75 % of ice loss in 20 000 years; deglaciation starts at 400 000 years.(b) Model set 2 of post-LGM deglaciation accounting for 95 % of ice loss in 10 000 years; deglaciation starts at 100 000 years.(c) The post-LIA deglaciation model set accounting for 10 % of ice loss in 400 years (shaded-blue region); deglaciation starts at 100 000 years.Shaded-blue regions highlight the modelled deglaciation intervals.Please note that the time axis in (a), (b), and (c) are different, and the post-LGM models account for longer timescales.

Table 2 .
Maximum uplift rates derived from the numerical models with a thermal anomaly (TA) of 0, 50, 100, 150, and 200 °C for Model set 1 (a) and Model set 2 (b) of post-LGM deglaciation and the post-LIA deglaciation model set (c).The t = 0 is the time step immediately before the beginning of deglaciation, and the other selected time steps show how the uplift rates change during the deglaciation until it is over for the post-LGM (a, b) and post-LIA (c) deglaciation intervals.Figure7is a plot of the maximum uplift rate vs. time calculated for each time step in all numerical models.= 1000 yr t = 5000 yr t = 10 000 yr t = 15 000 yr t = 20 000 yr t = 25 000 yr = 1000 yr t = 5000 yr t = 10 000 yr t = 15 000 yr t = 20 000 yr t = 25 000 yr , occurs.Under these conditions, however, our experiments pro- https://doi.org/10.5194/se-15-387-2024Solid Earth, 15, 387-404, 2024