Fast uplift in the southern Patagonian Andes due to long- and short-term deglaciation and the asthenospheric window underneath

An asthenospheric window underneath much of the South American continent increases the heat flow in the southern Patagonian Andes where glacial–interglacial cycles drive the building and melting of the Patagonian Icefields since the latest Miocene. The Last Glacial Maximum (LGM) was reached ∼26 000 yr BP (years before present). Significant deglaciation onsets between 21 000 and 17 000 yr BP were subject to an acceleration since the Little Ice Age (LIA), which was ∼400 yr BP. Fast uplift rates of up to 41±3 mm yr⁻¹ are measured by global navigation satellite system (GNSS) around the Southern Patagonian Icefield and are currently ascribed to post-LIA lithospheric rebound, but the possible longer-term post-LGM rebound is poorly constrained. These uplift rates, in addition, are 1 order of magnitude higher than those measured on other glaciated orogens (e.g. the European Alps), which raises questions about the role of the asthenospheric window in affecting the vertical surface displacement rates. Here, we perform geodynamic thermo-mechanical numerical modelling to estimate the surface uplift rates induced by post-LIA and post-LGM deglaciation, accounting for temperature-dependent rheologies and different thermal regimes in the asthenosphere. Our modelled maximum post-glacial rebound matches the observed uplift rate budget only when both post-LIA and post-LGM deglaciation are accounted for and only if a standard continental asthenospheric mantle potential temperature is increased by 150–200 °C. The asthenospheric window thus plays a key role in controlling the magnitude of presently observed uplift rates in the southern Patagonian Andes.

1 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, 1996, 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, 1996, 2004; Mitrovica and Forte, 1997; Kaufman and Lambeck, 2002; Stuhne and Peltier, 2015; van der 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, 1995, 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).

The southern Patagonian Andes in the South American continent are located above a transition zone between the subducting Antarctic and Nazca plates and a wide asthenospheric window where hot buoyant asthenospheric mantle upwells (Fig. 1a; Cande and Leslie, 1986; Breitsprecher and Thorkelson, 2009; Russo et al., 2010, 2022; Dávila et al., 2018; Ávila and Dávila, 2018, 2020; Mark et al., 2022; Ben-Mansour et al., 2022). The Chile triple junction (CTJ) at ∼46° S delimits the surface tip of the asthenospheric window, which opened during the last ∼16 Myr from south to north (Ramos and Kay, 1992; Breitsprecher and Thorkelson, 2009). First-order effects of the asthenospheric flow on the surface continental geology are the inhibition of arc volcanism in favour of retroarc magmatism, reduction in shortening to null or very minor, and rock uplift and exhumation (Ramos and Kay, 1992; Ramos, 2005; Lagabrielle et al., 2004, 2010; Breitsprecher and Thorkelson, 2009; Guillaume et al., 2009, 2013; Scalabrino et al., 2010; Lange et al., 2014; Georgieva et al., 2016, 2019; Ávila and Dávila, 2020; Mark et al., 2022; Ávila et al., 2023; Muller et al., 2023). Rock uplift due to asthenospheric upwelling occurs through dynamic and thermal effects (Guillaume et al., 2009, 2013; Conrad and Husson, 2009; Flament et al., 2013; Sternai et al., 2016b; Dávila et al., 2018; Ávila and Dávila, 2020; Faccenna and Becker, 2020; Mark et al., 2022). Dynamic uplift occurs above zones of viscous convection of the asthenospheric mantle, generating long wavelength (>300 km) deformation of the lithosphere through vertical stresses (Hager and O'Connell, 1981; Flament et al., 2013, 2015; Sternai et al., 2016b; Ávila and Dávila, 2020; Faccenna and Becker, 2020). This effect is difficult to measure because vertical stresses in the lithosphere also occur due to lithospheric tectonics and to the surface mass redistribution of glaciers, lakes, and sediments (Lachenbruch and Morgan, 1990; Molnar and England, 1990; Watts, 2001). Dynamic uplift was estimated to be between 0.02 and 0.15 mm yr⁻¹ in the last 3 Myr over an area of about 100 000 km² around the CTJ latitude (Guillaume et al., 2009, 2013; Flament et al., 2015; Ávila and Dávila, 2020; Ávila et al., 2023). The thermal component is expressed by an increase in temperatures and by a shallowing of the lithosphere–asthenosphere boundary where asthenospheric mantle upwells (Ávila and Dávila, 2018, 2020; Russo et al., 2010, 2022; Mark et al., 2022; Ben-Mansour et al., 2022), decreasing the integrated elastic lithospheric thickness and generating uplift and a higher surface heat flow than in normal subduction zones (Ranalli, 1997; Flament et al., 2015; Ávila and Dávila, 2018, 2020; Ávila et al., 2023). The heat flow was calculated as >100 mW m⁻² near the CTJ, ∼70–90 mW m⁻² in the centre of the asthenospheric window (∼50° S), and 50–60 mW m⁻² near its northern boundary (∼46° S) (Ávila and Dávila, 2018). Uplift due to lithospheric thinning was estimated to be ∼0.3 mm yr⁻¹ since the middle Miocene in the southern Patagonian Andes (Pedoja et al., 2011; Ávila and Dávila, 2020; Ávila et al., 2023; Ding et al., 2023).

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 from DeMets 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 in Lange et al. (2014).

The Patagonian Ice Sheet covered the southern Patagonian Andes between ∼47 000 and ∼17 000 yr BP, extending from latitudes 38 to 55° S with an estimated area of ∼490 000 km² (Fig. 1a); volume of ∼550 000 km³; and average and maximum thickness of 1100 and 2500 m, respectively, based on preserved glacial geomorphology, stratigraphy, palaeoecology, and geochronological data (Moreno et al., 1999, 2015; McCulloch et al., 2000, 2005; Hulton et al., 2002; Rabassa, 2008; Glasser et al., 2004, 2005, 2008, 2016; Glasser and Jansson, 2008; Hein et al., 2010; Boex et al., 2013; Strelin et al., 2014; Bourgois et al., 2016; Martinod et al., 2016; Kaplan et al., 2016; Bendle et al., 2017; Thorndycraft et al., 2019; Reynhout et al., 2019; Davies et al., 2020; Yan et al., 2022). The LGM in southern Patagonia is estimated around 26 000 yr BP, but the beginning of significant glacial retreat occurred between 21 000 and 17 000 yr BP (Hulton et al., 2002; Hein et al., 2010; Glasser et al., 2011; Davies and Glasser, 2012; Moreno et al., 2015; Bendle et al., 2017; Reynhout et al., 2019; Davies et al., 2020). The long-term ice loss rate is uncertain, but more than 75 % of ice was certainly lost since the LGM, and some models predicted more than 95 % of ice loss with separation between the Southern Patagonian Icefields (SPI) and the Northern Patagonian Icefields (NPI) in the first 5000 to 10 000 years of post-LGM deglaciation (McCulloch et al., 2000; Hulton et al., 2002; Boex et al., 2013; Bourgois et al., 2016; Thorndycraft et al., 2019; Davies et al., 2020). A glacial minimum must have been attained around 13 000 yr BP, but several glacial advances were recorded since that time, and the last recorded one was the Little Ice Age (LIA) with an apex around 1630 CE, well dated by terminal moraines around the present-day NPI and SPI (Ivins and James, 1999, 2004; McCulloch et al., 2000; Glasser et al., 2004, 2008, 2011; Davies and Glasser, 2012; Strelin et al., 2014; Kaplan et al., 2016; Reynhout et al., 2019; Davies et al., 2020). Recent mass balance measurements in the Patagonian Icefields – e.g. the Shuttle-Radar Topography Mission (SRTM) or Gravity Recovery and Climate Experiment (GRACE) – often present discrepancies but consistently show an increasing ice loss from ∼15 Gt yr⁻¹ between ∼1940–2000 to ∼25 Gt yr⁻¹ between ∼2000–2012 (Aniya, 1996; Aniya et al., 1997; Rignot et al., 2003; Chen et al., 2007; Ivins et al., 2011; Jacob et al., 2012; Willis et al., 2012; Gómez et al., 2022). Currently, the SPI covers an area of ∼13 219 km² with a volume of 3632±675 km³, whereas the NPI covers an area of ∼3976 km² with a volume of 1124±260 km³ (Fig. 1). The present-day ice thickness may reach up to ∼2000 m in deep glacial valleys (Millan et al., 2019).

GNSS-measured data show ongoing vertical rock uplift rates between 18±3 and 41±3 mm yr⁻¹ in the northern part (18–50.5° S) of the SPI (Fig. 1b), decreasing to values between 2±6 and 17±5 mm yr⁻¹ 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³ 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¹⁸ Pa s) and a thin elastic lithosphere (∼35 km thick) (Ivins and James, 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 and James, 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 post-glacial 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.

2 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).

2.1 Numerical model

We use a fully coupled thermo-mechanical, visco-elasto-plastic 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), and Muller et al. (2022). The continuity equation allows for the conservation of mass during the displacement of a geological continuum:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\left(\mathit{\rho }\mathit{v}\right)=\mathrm{0},\end{array}$$

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:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial x_i}}+\mathit{\rho }{g}_i=\mathit{\rho }\left({\displaystyle \frac{\partial {\mathit{v}}_i}{\partial t}}+v_j{\displaystyle \frac{\partial {\mathit{v}}_i}{\partial x_j}}\right),\end{array}$$

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:

$$\begin{array}{}\text{(3)}& \mathit{\rho }{C}_P{\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}-\mathrm{div}\left(c\mathrm{\nabla }T\right)+\mathit{v}\mathrm{\nabla }T=H_{\mathrm{r}}+H_{\mathrm{s}}+H_{\mathrm{a}}+H_{\mathrm{l}},\end{array}$$

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. $H_{\mathrm{a}}\propto \frac{\mathrm{D}P}{\mathrm{D}t}$, $H_{\mathrm{s}}={\mathit{\sigma }}_{ij}^{\prime }{\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{viscous}\right)}^{\prime }$, where ${\mathit{\sigma }}_{ij}^{\prime }$ is the deviatoric stress tensor and ${\dot{{\mathit{\epsilon }}^{\prime }}}_{ij\left(\mathrm{viscous}\right)}$ is the viscous deviatoric strain rate tensor. Ductile deformation is thermally activated generating viscous flow, which is calculated according to the material shear viscosity, η:

$$\begin{array}{}\text{(4)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{A}_d\mathit{\sigma }{{}^{\prime }}_{\mathrm{II}}^{n-\mathrm{1}}}}\exp \left({\displaystyle \frac{E_{\mathrm{a}}+P{V}_{\mathrm{a}}}{RT}}\right),\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{II}}^{\prime }$ 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. ${\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{viscous}\right)}^{\prime }$ is computed as

$$\begin{array}{}\text{(5)}& {\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{viscous}\right)}^{\prime }={\displaystyle \frac{{\mathit{\sigma }}_{ij}^{\prime }}{\mathrm{2}\mathit{\eta }}}-{\mathit{\delta }}_{ij}{\mathit{\eta }}_{\mathrm{bulk}}{\dot{\mathit{\epsilon }}}_{kk},\end{array}$$

where δ_ij is the Kronecker delta, ${\dot{\mathit{\epsilon }}}_{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, ${\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{elastic}\right)}^{\prime }$, as

$$\begin{array}{}\text{(6)}& {\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{elastic}\right)}^{\prime }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\mu }}}{\displaystyle \frac{\mathrm{D}{\mathit{\sigma }}_{ij}^{\prime }}{\mathrm{D}t}},\end{array}$$

where μ is the shear modulus and $\frac{\mathrm{D}{\mathit{\sigma }}_{ij}^{\prime }}{\mathrm{D}t}$ is the objective co-rotational 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:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{yield}}=C+\sin \left(\mathit{\phi }\right)P,\end{array}$$

where C is cohesion and φ is the effective internal friction angle. The plastic strain rate tensor, ${\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{plastic}\right)}^{\prime }$, is defined as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{plastic}\right)}^{\prime }=\mathrm{0}\mathrm{for}{\mathit{\sigma }}_{\mathrm{II}}<{\mathit{\sigma }}_{\mathrm{yield}},\\ & {\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{plastic}\right)}^{\prime }=X{\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}^{\prime }}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{II}}}}\mathrm{for}{\mathit{\sigma }}_{\mathrm{II}}\ge {\mathit{\sigma }}_{\mathrm{yield}},\end{array}\end{array}$$

where X is the plastic multiplier which satisfies the plastic yielding condition σ_II=σ_yield. The bulk strain rate tensor, ${\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{bulk}\right)}^{\prime }$, integrates the viscous, elastic, and plastic deformation:

$$\begin{array}{}\text{(9)}& {\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{bulk}\right)}^{\prime }={\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{viscous}\right)}^{\prime }+{\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{elastic}\right)}^{\prime }+{\dot{\mathit{\epsilon }}}_{ij\left(\mathrm{plastic}\right)}^{\prime }.\end{array}$$

2.2 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 and Dá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 and Dávila, 2018, 2020). The initial geotherm is piecewise linear, resulting from an adiabatic temperature gradient of 0.5 °C km⁻¹ 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., 2016, 2019; Stevens Goddard and Fosdick, 2019; Ávila et al., 2023; Muller et al., 2023). Rocks' rheological properties are listed in Table 1.

Table 1 Material properties used in the numerical experiments.

${\mathit{\rho }}_{\mathrm{0}}^{\mathrm{s}}$ is the standard densities of solid rocks; E_a is the activation energy; V_a 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; C_p is the specific heat capacity; H_r 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 (Ranalli, 1995; Hirschmann, 2000; Johannes, 1985; Turcotte and Schubert, 2002).

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).

Figure 2 Reference numerical model setup. (a) Thermo-mechanical numerical model domain with rheological layers (Table 1), isotherms (white lines), and yield strength ($\mathrm{\Delta }\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{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).

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⁻³ (Harvey et al., 2017) (Table 1), and we compute the surface load through time, L, as

$$\begin{array}{}\text{(10)}& L={\mathit{\rho }}_{\mathrm{ice}}g{h}_{\mathrm{ice}},\end{array}$$

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). We run two sets of experiments for post-LGM deglaciation. In Model set 1, 75 % of ice loss occurs in 20 000 years (i.e. 1500 m drop of ice thickness, Fig. 2b), thus assuming a conservative estimate of ice loss since the beginning of the LGM until the present day, simplifying the several glacial retreats and re-advances since the LGM (e.g. Glasser et al., 2004, 2008, 2011; Davies and Glasser, 2012; Strelin et al., 2014; Kaplan et al., 2016; Reynhout et al., 2019). In Model set 2, 95 % of ice loss occurs in 10 000 years (i.e. 1900 m drop of ice thickness, Fig. 2b), thus assuming faster deglaciation rates of the Patagonian Ice Sheet in the first half of post-LGM deglaciation (McCulloch et al., 2000; Hulton et al., 2002; Boex et al., 2013; Bendle et al., 2017; Thorndycraft et al., 2019; Davies et al., 2020). For the post-LIA deglaciation, we simulate 10 % of ice loss in 400 years (i.e.  200 m drop of ice thickness, Fig. 2c) using estimates of ice loss rates occurring since the nineteenth century (Aniya, 1996; Aniya et al., 1997; Rignot et al., 2003; Ivins and James, 1999, 2004; Chen et al., 2007; Dietrich et al., 2010; Ivins et al., 2011; Jacob et al., 2012; Willis et al., 2012). The pseudo-ice cap is 200 km wide for the post-LGM model sets 1 and 2, based on estimates of the LGM extent of the Patagonian Ice Sheet (e.g. McCulloch et al., 2000; Hein et al., 2010; Thorndycraft et al., 2018; Davies et al., 2020), and it is 70 km wide for the post-LIA model set, based on the estimates of the LIA maximum extent of the SPI (e.g. Glasser et al., 2011; Strelin et al., 2014; Kaplan et al., 2016; Reynhout et al., 2019) (Fig. 2a).

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_{\mathrm{curr}}-z_{\mathrm{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 and Dávila, 2018, 2020; Sternai, 2020; Mark et al., 2022) to mimic the presence of a slab window at depth.

3 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 (Figs. 4–6). 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²²–10¹⁹ Pa s in simulations with TA equal to 0 (reference model), 50, and 100 °C and between 10¹⁹–10¹⁶ 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²² to 10²⁰ Pa s), thereby decreasing the integrated lithospheric strength.

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 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 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 6 Surface uplift rates vs. distance for the post-LIA deglaciation model set. Panel (a) shows t=100 years of deglaciation, panel (b) shows t=200 years of deglaciation, panel (c) shows t=300 years of deglaciation, and panel (d) shows 400 years of deglaciation. Different line colours correspond to different TA.

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. Figure 7 is a plot of the maximum uplift rate vs. time calculated for each time step in all numerical models.

In Model set 1 for post-LGM deglaciation, when TA equals 0 (reference model) the maximum uplift rate is <1 mm yr⁻¹ during the first 5000 years of the deglaciation, increasing gradually up to 9.5 mm yr⁻¹ 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⁻¹, 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⁻¹ 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⁻¹ 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⁻¹ when TA ≤ 100 °C, whereas it drops to 0 mm yr⁻¹ almost immediately when TA is 150 or 200 °C (Fig. 7a).

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.

In Model set 2 for post-LGM deglaciation, the maximum uplift rate is less than 2 mm yr⁻¹ 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⁻¹ 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⁻¹, respectively, when TA is 0, 50, and 100 °C, whereas it reaches up to between 36 and 41 mm yr⁻¹ 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⁻¹ (Fig. 7b and Table 2b), whereas it quickly drops to <2 mm yr⁻¹ 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⁻¹ 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⁻¹ 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⁻¹ when TA is equal to 0, 50, and 100 °C and ∼14 and ∼25.5 mm yr⁻¹ 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⁻¹ 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 higher-magnitude post-glacial rebound which, however, takes much longer to re-equilibrate to 0 mm yr⁻¹ after the end of the deglaciation than a cooler and more viscous asthenosphere.

4 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 and James, 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 and James, 1999, 2004; Hulton et al., 2002; Klemann et al., 2007; Ivins et al., 2011; Boex et al., 2013). Measurements of regional erosion rates since the LGM are between 0.02 and 0.83 mm yr⁻¹ (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 low-temperature 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), and Muller et al. (2023) suggest using values between 0.1 and 1 mm yr⁻¹ from 7 to 4 Ma, followed by a period of erosional quiescence (<0.1 mm yr⁻¹) and a possible increase to 1 mm yr⁻¹ 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⁻³. 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., 2013, 2018; Herman and Brandon, 2015; Fernandez 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⁻¹; 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 thermo-mechanical geodynamic experiments provide realistic uplift rates (Figs. 4–7) that one can compare to current geodetic observations. Following the example of previous studies (Ivins and James, 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 and Dá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⁻¹ 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⁻¹ and 20 mm yr⁻¹, 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¹⁸ Pa s (Ivins and James, 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¹⁹ Pa s, with the lowest viscosity value of 10¹⁶ 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), occurs. Under these conditions, however, our experiments provide max uplift rates between 14 and 26 mm yr⁻¹ 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⁻¹. 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⁻¹ 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-than-normal 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 and Dá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., 2010, 2022; Mark et al., 2022; Ávila and Dá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 and James, 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)

$$\begin{array}{}\text{(11)}& {\mathit{\tau }}_{\mathrm{r}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }v}{g\mathit{\lambda }}},\end{array}$$

where v is the asthenosphere viscosity, λ is the width of the ice sheet, and g is the gravity acceleration. Using 10${}^{\mathrm{16}}<v<{\mathrm{10}}^{\mathrm{18}}$ Pa s and λ=200 km leads to $\sim \mathrm{2000}<{\mathit{\tau }}_{\mathrm{r}}<\sim \mathrm{200}\mathrm{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⁻¹ in the next century if the currently observed ice loss rate of at least −20 Gt yr⁻¹ in the SPI (Willis et al., 2012) will continue until the total meltdown of the ice sheet in ∼200 years.

5 Conclusions

We propose that rock uplift rates of up to 40 mm yr⁻¹ 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¹⁶ Pa s). Our thermo-mechanical visco-elasto-plastic 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, 2022; Ávila and Dávila, 2018, 2020; Mark et al., 2022; Ben-Mansour et al., 2022).