2.1 Digital elevation models (DEMs)

Table 1Different DEMs used in this study produced from satellite data.

Download Print Version | Download XLSX

We use seven different DEMs derived from different satellite missions between 2000 and 2016 (see Table 1). The SRTM C-band radar DEM from mid-February 2000 (Farr et al., 2007) is used as the steady-state reference of the two glaciers for reconstructing bedrock topography. A Pléiades optical satellite stereo DEM from 1 October 2016, after the collapse, allows us to evaluate the modelled bedrock reconstruction over the detachment zone. We compute ice emergence velocities by differencing pre-collapse high-resolution DEMs from TanDEM-X, Spot7, and WorldView data and correct these for mass balance following the approach described in Gilbert et al. (2016) (Fig. 2). The effect of uncertainty linked to radar penetration in the TanDEM-X data should be minimized when comparing same wavelength data (X-band) at similar times of the year. Changes in penetration depth between the TanDEM-X data of 2011 (early June) and 2013 (mid April) due to different snow wetness should also be limited because surface melting in the accumulation area of the Aru glaciers only occurs from around mid-June (Kääb et al., 2018). X-band penetration into glacier ice (i.e. the Aru glaciers ablation areas) is very limited anyway (Dehecq et al., 2016). Comparing Spot7 (2015) and TanDEM-X (2014) elevations likely introduces uncertainty from TanDEM-X penetration in the accumulation area leading to higher apparent emergence velocities in this part (visible in Fig. 2). This problem only influences our friction reconstruction in the upper parts of the glaciers but not in the detachment areas. Details on DEM accuracies and acquisition methods can be found in Kääb et al. (2018).

Figure 2Mean emergence velocities obtained by differencing elevation changes from repeat DEMs and modelled mass balances during different periods prior to the collapses. Steady-state velocities in the first panel are modelled. Orange dashed lines indicate the detachment outline.

Download

2.2 Field observations

We investigated glacier till properties by analysing samples collected from the Aru-1 avalanche deposit in the gorge close to the former glacier tongue. We collected these samples 1 year after the collapse on 18 July 2017. Rainy conditions on that day highlighted the behaviour of the surrounding lithology that quickly turned into a soft and unstable slurry in the presence of water. Additional information about our samples can be found in Figs. S1–S3 in the Supplement.

3.1 Mass balance

Our mass balance model for the two Aru glaciers is based on a degree-day model described in Gilbert et al. (2016). It has been calibrated for the Aru glaciers by using satellite-derived glacier mass balances and is fed by ERA-Interim climate reanalysis (Kääb et al., 2018). The model output, taken from Kääb et al. (2018), provides the spatio-temporal distribution of surface mass balance, firn thickness, and available surface meltwater for percolation/refreezing in the firn to constrain the thermo-mechanical model below.

3.2 Thermo-mechanical model

Our thermo-mechanical ice-flow model is based on the Stokes equation coupled with an energy equation using the enthalpy formulation (Aschwanden et al., 2012; Gilbert et al., 2014). Changes in the glacier geometry are computed using a free surface equation (Gilbert et al., 2014). We adopt a pure viscous isotropic ice rheology following Glen's flow law (Cuffey and Paterson, 2010). The model is solved using the finite-element software Elmer/Ice (Gagliardini et al., 2013). Parameters and variables of our model set-up are summarized in Table 2.

Table 2Variables and parameters of the thermo-mechanical model.

^a Cuffey and Paterson (2010). ^b Calonne et al. (2011). ^c Gilbert et al. (2014).

Download Print Version | Download XLSX

We adopt a linear friction law as a basal boundary condition for the Stokes equation:

where τ_b is the basal shear stress (MPa), u_s is the sliding velocity (m yr⁻¹), and β is the friction coefficient (MPa yr m⁻¹). This coefficient is inverted using a control inverse method to minimize a cost function defined from the misfit with measured surface data and a regularization term (Gillet-Chaulet et al., 2012; Gagliardini et al., 2013). Following Gilbert et al. (2016) we used the emergence velocity $U_{N_z}^{\mathrm{obs}}$ (m yr⁻¹) to compute this cost function:

where $U_{N_z}=(\mathit{u}\cdot \mathit{N})N_z$ is the modelled emergence velocity (m yr⁻¹), u is the surface velocity vector (m yr⁻¹), $\mathit{N}=(N_x,N_y,N_z)$ is the unit vector normal to the surface, Γ_s is the surface boundary, J_reg is the regularization term, Γ_b is the bedrock surface boundary, and λ is a positive number. The emergence velocity is obtained by removing the mean modelled mass balance from the elevation change rate measured from our repeat satellite-derived DEMs over the same periods (Fig. 2):

where $\partial h/\partial t$ is the measured elevation change rate (m yr⁻¹) and M is the mean surface mass balance during the corresponding period (m yr⁻¹).

The surface boundary condition is set as a stress-free boundary for the Stokes problem and a Dirichlet condition for the enthalpy equation. In order to take into account water percolation and refreezing within the firn, we follow the approach by Gilbert et al. (2015), using a 6-month time step. Latent heat due to refreezing is released every year during the summer time step. The firn-thickness distribution is estimated from the mass balance model following Gilbert et al. (2016) and the firn density is computed using a linear density profile set to the following:

where ρ is the density (kg m⁻³) at depth d (m), ρ₀ is the surface density, ρ_ice is the ice density, and H_firn is the firn thickness (m).

The lateral boundary condition is set to a no-flux condition for both the Stokes and enthalpy equations. We assume a basal heat flux of $\mathrm{8.0}\times {\mathrm{10}}^{-\mathrm{2}}$ W m² for the enthalpy equation according to heat flux measured in boreholes at the Guliya ice cap (6200 m a.s.l., 200 km north of the Aru Range) (Thompson et al., 1995) and modelled geothermal heat fluxes for the region (Tao and Shen, 2008).

3.3 Modelling strategy for the steady-state glaciers

The first step of modelling the dynamics and thermal regime of the Aru glaciers is to establish a steady-state glacier as the initial condition for 1970 (start of the climatic reanalysis used). Landsat satellite images of the glacier area and the mass balance model indicate that the two glaciers were close to equilibrium from 1970 to 1995 (Kääb et al., 2018). We therefore assume that the surface topography measured in February 2000 by the SRTM mission (oldest available DEM) is representative of the glaciers, being in equilibrium with the mean climate over this period, although the positive mass balance between 1995 and 2000 probably thickened the glacier by a few metres in the accumulation area. We use the mean mass balance between 1980 and 1995 as an equilibrium mass balance considering that modelled mass balance is close to steady state during this period, before becoming positive from 1995 to 2008 (Kääb et al., 2018).

We first run the model on a 2-D flow line until a steady state is reached, deriving bedrock topography in the detached parts from a post-collapse Pléiades DEM and by reconstructing the bed at the upper glacier parts, assuming a constant basal shear stress (plastic approximation; Cuffey and Paterson, 2010). This initial step allows for the first approximation of the steady-state thermal regime which we presented in Kääb et al. (2018). We then use the 10 m depth temperature modelled with the flow line model to define the steady-state surface enthalpy as a function of elevation, which includes meltwater percolation and refreezing effects. This relationship is used to define a Dirichlet surface boundary condition for enthalpy in order to solve the steady-state thermal regime of the glaciers in 3-D in the bedrock inversion procedure (Sect. 3.3.1). Because the effects of meltwater percolation and refreezing are already included in the surface enthalpy value, there is no need to solve for these effects in diagnostic runs. The final 3-D steady-state glacier solution is obtained by running a transient simulation using the inverted bedrock topography and solving water percolation and refreezing until surface topography and the enthalpy field reach equilibrium with the imposed climatic condition.

3.3.1 Reconstructing bedrock topography

Using constant climatic conditions associated with the balanced glacier conditions corresponding to the SRTM DEM, we determined the bedrock topography, allowing the best match between modelled and observed (i.e. SRTM DEM) surface topography (van Pelt et al., 2013). For this purpose, we ran a 3-D transient simulation assuming no sliding, fixed surface topography (SRTM DEM), and constant surface forcing (mass balance and enthalpy). The no-sliding assumption is likely a good assumption in 2000 since the glacier was not surging at this time (Kääb et al., 2018). Mesh horizontal resolution is set to about 50 m with 15 vertical layers. The evolution of the free surface is taken into account by varying the basal mesh elevation instead of the surface elevation. The mesh surface topography thus remains constant, while the bed topography is updated by solving the equation

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \left(z_{\mathrm{bed}}\right)}{\partial t}}+{\mathit{v}}_{\mathrm{s}}\cdot \mathbf{\nabla }\left(z_{\mathrm{bed}}-z_{\mathrm{bed}\mathrm{0}}\right)={\mathit{v}}_{\mathrm{s}}\cdot \mathbf{\nabla }{z}_{\mathrm{s}\mathrm{0}}-M-w_{\mathrm{s}},\end{array}$$

where z_bed is the bedrock elevation (m), v_s is the surface horizontal velocity (m yr⁻¹), z_bed0 is the initial bedrock topography (m), M is the surface mass balance (m yr⁻¹), w_s is the vertical component of the surface velocity (m yr⁻¹) and z_s0 is the measured surface elevation (m). The right side of Eq. (6) vanishes once bedrock topography satisfies the required topography to keep z_s0 constant for a given mass balance M. The advective term on the left side of Eq. (6) allows z_bed to be smoothed in the flow direction. The enthalpy field is solved at each time step by solving the steady-state equation for the current velocity field and mesh. We start with a uniform 200 m ice thickness (rough maximum expected thickness on the glacier) and run the model until a steady bedrock topography is reached. This generates a new z_bed0 value to re-run the model until a new steady state is reached. After only two iterations, we validate the modelled bedrock topography by running the model with the new fixed bedrock topography and free surface evolution. The resulting surface topography is in excellent agreement with the measured one (Fig. 3), indicating that our method to infer the bedrock topography works well for these two glaciers.

Figure 3Pléiades image of Aru-1 (left) and Aru-2 (right) glaciers after collapse, with topographic profiles 1 to 6 plotted for both glaciers (Copyright CNES 2016, Distribution Airbus D&S). The topographic profiles 1 to 6 show the measured surface topography in 2000 (SRTM, in red) and 2016 after the collapse (Pléiades, in yellow). These profiles are compared with the modelled bedrock (in purple) and surface (in green) topographies. The coloured dots on the Aru-2 glacier show the location of specific points of the profiles in the Pléiades image. Those points correspond to locations where our reconstruction matches the Pléiades DEM and where bedrock should thus be visible on the Pléiades image (no ice debris). Grey shading indicates the detached parts according to the Pléiades DEM compared to SRTM.

Download

We use the opportunity provided by the exposed detachments to compare the reconstructed bedrock topography with the measured Pléiades DEM from after the collapses (Fig. 3). On the Aru-2 glacier, the points at which bedrock is clearly visible in the Pléiades images match well with the locations where our reconstructed bedrock topography matches the Pléiades DEM (dots in Fig. 3). Elsewhere in the Aru-2 glacier detachment zone, the modelled bedrock is deeper than the observed surface elevation; this is likely due to the remaining ice debris overlying the actual bedrock, so the Pléiades DEM elevations are expectedly higher. This is confirmed by the good continuity between the ground topography measured outside of the former glacier tongue and the one inferred from our bedrock reconstruction (see Fig. 3, profile 6). On the Aru-1 glacier, the reconstructed bedrock on profiles 2, 3, and 4 is systematically deeper than the Pléiades DEM, even on the steep side, close to the margin of the glacier where no ice remained after the collapse. This means that ice flow is not accurately modelled in this part, likely due to the premise of no sliding, which is probably not accurate considering that the glacier may have been temperate at its base here (see Sect. 4.1). The error in the modelled bed topography of the Aru-1 glacier is, however, <30 m and will only slightly affect the absolute value of the friction coefficient inferred during the instability development (see Sect. 4.2), but not its relative changes, which are the focus of this study. The assumption of no sliding should also affect the result on the Aru-2 glacier, which has a similar thermal regime, but where the no-sliding condition seems to work. This indicates the existence of different sliding conditions for the two glaciers prior to collapse, which is also supported by the friction inversion analysis presented in Sect. 4.2. In the upper parts of both glaciers, the no-sliding assumption is, however, supported by the friction inversion analysis.

4.1 Steady state configuration of the two glaciers

The Aru glaciers are representative of a cold and semi-arid climate regime, and thus would normally show little dynamic behaviour under mostly cold-ice conditions (below the pressure melting point). The steady-state equilibrium line is located around 5750 m a.s.l. (minimum glacier elevations around 5200 m a.s.l., and maximum elevations around 6100 m a.s.l.) with a maximum accumulation of 0.6 m.w.e. yr⁻¹ at 6100 m a.s.l. and a maximum ablation of −2.5 m.w.e. yr⁻¹ on the tongue (Fig. 4b). Both glaciers are composed of two similar catchments characterized by a smaller western branch that joins the main stream in the ablation area. The western branch of each glacier is thinner and less dynamic than the main branch that collapsed in summer 2016 (Fig. 4a). Maximum surface horizontal velocity reached 20 m yr⁻¹ in the Aru-2 glacier, which has to accommodate higher ice flux than the Aru-1 glacier due to a wider accumulation area (1.7 km² vs 1.2 km²) converging in a similarly narrow gorge.

Figure 4Modelled steady-state horizontal surface velocities (a) and surface mass balance (b) for Aru-1 and Aru-2 glaciers. The black contours in (a) are modelled steady-state glacier outlines. The white contour in (b) is the glacier outline as mapped from 2015 images. Orange dashed lines indicate the detachment outline.

Download

As previously concluded by Kääb et al. (2018), our results show that the main branches of the two glaciers are characterized by a polythermal structure with a cold accumulation zone above 5900 m a.s.l. and a temperate-based ablation area surrounded by cold ice (Fig. 5). However, through the more accurate bedrock topography derived in this study and the 3-D approach, we show here that the temperate zones likely extended into significantly larger areas beneath the detachments than previously thought. Temperate ice forms in the lower part of accumulation zones due to a significant amount of percolation and refreezing of meltwater, which increases the temperature of the near-surface firn. This warmer ice is then advected into the ablation zone contributing, together with basal heat flux, to temperate basal conditions in the lower parts of the two glaciers. Cold surface conditions due to absence of water percolation in the ablation zone (cold impermeable ice) lead to a significant cold surface layer that eventually reaches the glacier base in the shallowest zones of the glacier tongues (Fig. 5). The western branches of the two glaciers have a significantly smaller temperate area with an ablation zone that is almost entirely cold-based (Fig. 5). This thermal structure may explain why the western branches remained stable after the collapses even though each branch lost its downstream supporting buttress formed by the detached glacier tongues. The modelled spatial extent of the temperate basal ice under steady-state conditions coincides with the detached areas and indicates that friction changes leading to the collapse occurred in temperate ice rather than being produced by a change from cold to temperate thermal conditions at the glacier beds. However, the large amount of cold ice, especially along the side of the gorge, could have provided significant lateral drag that built up high driving stress, which was able to balance gravitational force under the frictional change at the temperate parts of the beds.

Figure 5Modelled steady-state temperature on the Aru-1 and Aru-2 glaciers. Left panel shows basal temperature with black hatched lines showing temperate areas. The inset highlights temperate-based (red) and cold-based (blue) areas. Orange dashed lines indicate the detachment outline. Panel (b) shows 2-D temperature profiles 1 to 6 as indicated in panel (a) (red lines). Profiles include Pléiades 2016 elevation profiles (orange lines). The dashed black lines indicate the cold-temperate transition surface. Note that vertical scale is exaggerated in profiles 1 and 4.

Download

4.2 Basal friction change since 2011

The surge-like behaviour of the two glaciers identified from DEM comparison in Kääb et al. (2018) documents a change in the glacier dynamics during the 5 years prior to the twin collapses. By removing the elevation change due to surface mass balance we quantify the emergence velocity for constraining the basal friction parameter (Gilbert et al., 2016) for different periods: 2011–2013, 2013–2014, 2014–2015, and September–November 2015 (Table 1, Fig. 2). Our results highlight contrasting behaviour between the Aru-1 and Aru-2 glaciers, where friction decreased progressively in magnitude through time in both glaciers but over significantly different areas (Fig. 6). Frictional changes over the 5 years prior to collapse are also more significant on the Aru-1 glacier, resulting in a higher increase in surface velocity than on the Aru-2 glacier (Fig. 7). Similarly inferred friction for the Aru-2 glacier for annual means (2011–2013 and 2013–2014) and a 2-month mean (September–November 2015) indicates low seasonal variability of the basal condition. Similarly, modelled surface velocities on the Aru-1 glacier in September–November 2015 correlate well with those measured for January–April 2016 using satellite image correlation (Kääb et al., 2018) (Fig. 7f), also indicating low seasonal variability.

Figure 6Friction coefficient β inferred from emergence velocity during different periods prior to the collapse. Orange dashed lines indicate the detachment outline.

Download

Figure 7(a–e) Modelled mean horizontal velocities for the two glaciers at steady state and for the periods 2011–2013, 2013–2014, 2014–2015, and September to November 2015. (f) Measured horizontal velocities from Sentinel-2 image correlation between January and April 2016 (adapted from Kääb et al., 2018). Orange dashed lines indicate the detachment outline.

Download

4.3 Force balance analysis

To evaluate how resisting forces acted and evolved to balance the driving forces, we compute the mean basal shear stress during different periods from the inverse method. We therefore assume that basal shear stress is mainly constrained by the global stress balance and should not be influenced by the sliding law that we used (Eq. 1) (Joughin et al., 2004; Minchew et al., 2016). The steady-state condition shows a basal shear stress between 100 and 200 kPa in both glaciers with mean shear stresses of 137 and 150 kPa for the Aru-1 and Aru-2 glaciers (Fig. 8a). In comparison, mean driving stresses are 152 kPa (Aru-1) and 213 kPa (Aru-2), indicating that 10 % (Aru-1) and 30 % (Aru-2) of the driving force are accommodated by normal force along the sidewalls. These levels of driving stress are at the higher end of the observed range of driving stresses on mountain glaciers (Cuffey and Paterson, 2010) and reflect the presence of strong resisting forces due mainly to cold-ice conditions combined with the resistance of the valley walls.

Figure 8Modelled basal shear stress at steady state (a) and in November 2015 (b). Panel (c) is the difference between (b) and (a).

Download

The inversion of mean elevation changes between September and November 2015 (Fig. 8b) reveals that basal shear stresses on the Aru-1 glacier decreased to only 20 to 10 kPa in large areas, and basal resistance was mainly achieved by a few sticky spots (Stokes et al., 2007) in the detachment zone where shear stresses exceeded 250 kPa. Along the left bank of the glacier, close to the terminus of the Aru-1 glacier, shear stress was about 6–7 kPa and was not more than 15 kPa at the terminus. In comparison, the Aru-2 glacier behaved differently with more localized friction changes that produced a smaller change in the distribution of the basal shear stress during the same period (Fig. 8c).

The analyses of the dynamics and force-balance evolution on an area restricted to the detachment zone (dashed red lines in Fig. 8) reveal both similarities and differences between the two events (Fig. 9). Further references below to “lateral stress” apply to the detachment zone and not to the whole glacier and refer to the stress provided by the shearing interface between the stable and the unstable part of the glacier (visible in Fig. 9d). On the one hand, as already highlighted in Fig. 7, the mean detachment velocity prior to collapse behaved differently for the two glaciers (Fig. 9a). While the Aru-1 glacier detachment significantly accelerated, following behaviour typical for slope failure (Voight, 1990) over several years (blue dashed line in Fig. 9a), the Aru-2 glacier showed very little acceleration. On the other hand, force balances evolved similarly in the two detachments with a large increase in lateral stresses along the detachment margin due to both an increase in the driving stress and reduction in basal friction (Fig. 9b). Interestingly, lateral resistance overcomes basal resistance in both detachments with a delay time (81 days) close to the actual delay between the two final collapses (66 days) (Fig. 9c). On the Aru-2 glacier, it seems that smaller changes in friction were compensated by a higher change in driving stresses resulting in a similar increase in stress at the detachment margin compared to the Aru-1 glacier (Fig. 9b). The difference in surface velocity response to these similar stress transfers was a consequence of different basal drag repartitions in the two glaciers. Basal drag decreased uniformly across the entire detachment of the Aru-1 glacier with the appearance of localized sticky spots, whereas basal drag decreased only in the higher part of the detachment of the Aru-2 glacier. This led to more intense bulging and a lower velocity increase (Kääb et al., 2018) due to the high-friction patch remaining in the tongue (Fig. 6).

Figure 9(a) Mean sliding speed of the detachment zone for Aru-1 (black) and Aru-2 (red) glaciers. The dashed blue lines show predicted speed following an empirical law of slope failure (Voight, 1990). (b) Force balance acting on the detachment of Aru-1 (black) and Aru-2 (red). (c) Ratio of resisting force over driving force for Aru-1 (black) and Aru-2 (red). Vertical lines show collapse dates in the three panels. (d) The Aru-2 glacier detachment zone on 4 October 2016 (Picture from T. Yao). The side resistance computed in the force balance analysis arises from lateral shearing on the right side of the detachment (black arrow).

Download

To evaluate the impact of the friction change on the mechanical property of the ice, we compute the maximal principal Cauchy stress and compare it with a threshold value set to 0.1 MPa (Krug et al., 2014) to identify the damage production (crevasse opening) (Krug et al., 2014; Pralong and Funk, 2005) (Fig. 10). The modelled stress fields clearly highlight zones where a progressive intensification of cracks opened around the detachment zone of the Aru-1 glacier (Fig. 10c) as observed on satellite images (Kääb et al., 2018); these fractures led to its final collapse. In comparison, the Aru-2 glacier again behaves differently, with less damage (cracks) that only affect the upper part of the detachment (Fig. 10c). This means that damage at the shear margin would have occurred suddenly in the Aru-2 glacier in 2016, which is confirmed by the observed sharp crack surrounding the detachment only a few days before the collapse (Kääb et al., 2018). In sum, the Aru-1 and Aru-2 glaciers underwent similar stress transfers, transitioning from basal drag to lateral shearing in their respective detachments, but showed different responses in terms of damage (i.e. crack production) and sliding speed due to different basal drag repartition. The Aru-1 glacier progressively evolved towards collapse, whereas the Aru-2 glacier accumulated stresses until a sudden release led to collapse. This indicates that critical stress transfers, precursory to such collapses, may occur without observable phenomena (i.e. surface velocity increase, crevassing) in the preceding years.

Figure 10Maximum principal Cauchy stress excess above damage initiation threshold at steady state (a) and prior to collapse (b) at the glacier surface. Panel (c) shows the difference between (a) and (b). Background image in (a) is a WorldView image from 2 December 2011 when the instability just started. Background image in (b) and (c) is a Spot 7 image from 21 September 2015, 1 year before collapse (copyright Airbus D&S). These results show a good match between predicted and observed crevasse formation in response to frictional changes.

Download

5.1 Result uncertainties

The modelled thermal regime is sensitive to basal heat flux, which is poorly constrained. However, sensitivity tests (see Supplement, Fig. S4) show that the temperate area remains stable for a basal heat flux between $\mathrm{6.0}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{1}}$ W m⁻² and disappears only at $\le \mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{2}}$ W m⁻². Measurements in the Guliya Ice Cap (Thompson et al., 1995) and reconstructions from Tao and Shen (2008) both give a value of $\mathrm{8.0}\times {\mathrm{10}}^{-\mathrm{2}}$ W m⁻², making a low value of $\le \mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{2}}$ W m⁻² very unlikely. Kääb et al. (2018) have also shown that firn thickness has a great influence on the modelled thermal regime around 5900 m a.s.l. where melting occurs in the accumulation zone. Firn thickness is, however, hard to estimate without field investigation; following Kääb et al. (2018), we applied an intermediate scenario where firn thickness linearly increases from the ELA to the glacier top, where it reaches a 15 m maximum. The sensitivity test showed that only very little firn thickness (<5 m at 6000 m a.s.l.) would lead to an almost cold glacier (Supplement Fig. S5). Nevertheless, the modelled thermal regime and the friction reconstruction, which are both almost independent from each other, are in good agreement with the localization of sliding and modelled temperate areas, lending confidence in our results despite uncertainties in basal heat flux and surface boundary conditions (see Sect. 5.2).

The uncertainty in the reconstruction of basal friction mainly depends on the accuracy of measured elevation changes, which is generally higher over longer time periods (increased signal-to-noise ratio), making the 2011–2013 reconstruction the most reliable one. The measured September–November 2015 elevation change is subject to a lower signal-to-noise ratio and is thus poorly resolved in the accumulation area (see Fig. 2). However, it leaves the reconstruction reliable on the detachment area where elevation changes are much more significant. A similar conclusion applies to the 2014–2015 reconstruction, where the upper part of the glacier is affected by penetration of the X-band signal, leading to an overestimation of the emergence velocities (see Sect. 2.1). However, the influence of this uncertainty on the modelled mass balance used to compute emergence velocity is also low in the detachment zone, since elevation changes due to surface mass balance are relatively small compared to the dynamical height changes linked to the surge-like instability (<20 %) (Fig. S6). Our results are, therefore, least affected by uncertainties and most reliable in the detachment area, which is the focus of this study. In addition, bedrock topography is well constrained in the detachment areas from the post-collapse Pléiades DEM, giving additional confidence in the friction reconstruction there.

Using emergence velocities to constrain basal friction is not a commonly used method and has been successfully tested only on a slow-flowing ice cap by Gilbert et al. (2016). We therefore provide additional validation of this method in the Supplement (Figs. S7–S9) by inverting the friction on both glaciers using horizontal velocities inferred from offset tracking obtained from repeat TerraSAR-X data in December 2013. This test reveals good agreement between our emergence-based approach and the more common method based on horizontal velocities. In particular, sliding zones are similarly localized in both methods. Using the inversion based on horizontal velocities as a reference, we estimate a sliding speed magnitude accuracy of 0.036 m day⁻¹ in the emergence-based inversion. Our additional validation test also indicates that using emergence velocities may provide for an improved constraint of the friction coefficient in accumulation areas. The reason for this is that the underlying data used in generating the emergence velocities (i.e. DEMs, modelled mass balance) are often more spatially resolved and cover larger areas on small mountain glaciers as opposed to measurements of horizontal displacements, which have problems over accumulation areas.

5.2 Frictional changes

Our results indicate that low friction below the Aru glaciers was not linked to seasonal variability of water pressure, which is often observed on glaciers elsewhere (Bartholomaus et al., 2008; Vincent and Moreau, 2016). Rather, it is likely associated with sustained change in the basal conditions caused by an accumulation of liquid water over several years prior to the collapse. Over a hard bed (Cuffey and Paterson, 2010), this would likely result in the existence of a subglacial lake, which is very unlikely here because the low friction zone on the Aru-1 glacier extended to the tongue and the lake should have drained in such a case. Furthermore, in temperate ice, high water pressure conditions are unstable over long time periods because they lead to channel formation that can efficiently drain water and decrease the pressure (Schoof, 2010). High water pressure in a cavity network would be also difficult to maintain in the Aru cases, since increasing sliding speeds tend to increase cavity size and decrease water pressure. These arguments suggest that basal friction under the Aru glaciers was probably controlled by processes associated with soft bed properties (Cuffey and Paterson, 2010).

Figure 11(a, b, c) Modelled temperate area (hatched zones) and −2 ^∘C isotherm at steady state compared with sliding speeds over different periods (background colour map). (d) Comparison of sliding location for the different periods.

Download

Comparison between sliding speed evolution and modelled basal steady-state temperature reveals a good correlation between the zones of sliding and temperate ice conditions and shows that the size of sliding areas remains similar over time (Fig. 11). This confirms that friction reduction since 2011 mainly occurred within zones that were already temperate areas, and that friction reduction is therefore not linked to a simple change from cold to temperate basal conditions. However, contrary to the Aru-1 glacier, Aru-2 appears to have been affected by a high-friction zone under its lower tongue, which the modelled basal temperatures are not able to explain as they indicate temperate, not cold conditions (Fig. 11). This zone of high friction explains the different behaviour observed from the two glaciers in terms of surface velocities and glacier advance. Indeed, a few months before the collapse, the Aru-2 glacier velocities were still low compared to the Aru-1 glacier (Fig. 7) and, while the Aru-1 glacier has advanced almost 200 m since July 2015, the front of the Aru-2 glacier remained stationary until the collapse (Kääb et al., 2018). Although the high-friction zone may have delayed the collapse of the Aru-2 glacier, it did not prevent it.

5.3 Role of the bedrock lithology and glacier till

Field observations after collapse and inspection of the detachment zones showed no presence of a hard-bed lithology beneath the glaciers, and no large boulders were observed in the forefields and avalanche deposits. Rather, extensive deposits of soft, unconsolidated and fine-grained lithologies were identified (Figs. S1–S3). We collected till samples from the Aru-1 glacier avalanche deposit and measured their grain-size distribution (Fig. 12). Mean values over the four samples in the avalanche path (Fig. 12) indicate till consisting of 14 % clay, 24 % silt, 44 % sand, and 18 % gravel. These samples are representative of the material found in the deposit and are likely also representative of the glacier till on which the glacier rested. We also observed a rather smooth-surfaced failure interface (i.e. detachment plane) suggesting a low bedrock roughness at the macro-scale (>1 m).

Figure 12Grain-size distribution measured in four glacier till samples collected in the Aru-1 deposit area (numbers 1 to 4 in the right panel). Background image from Google Earth.

Download

These findings confirm that the Aru glaciers rested on a soft bed, which likely played an important role in controlling the behaviour of the glaciers from the surge initiation to the collapse. For such bed types, basal motion is not controlled by ice flow around bedrock bumps (Weertman, 1964; Lliboutry, 1968) but rather by deformation in till (Truffer et al., 2001). The sustained very low basal drag under the Aru glaciers (<20 kPa) may be similar to ice stream mechanisms whereby water-saturated till enables fast flow at low driving stresses (≈20 kPa) (Cuffey and Paterson, 2010). It has been shown that glacier till behaves with a plastic rheology with a shear strength that is strongly dependent on the effective normal stress (Iverson et al., 1998; Clarke, 2005; Iverson, 2010). The use of a linear friction law in our inversion can be viewed as a parameterization, where β includes these physics and is only valid at the time of the inversion. Therefore, the change in friction coefficient β can be interpreted in terms of a plastic till. This behaviour was found to be well described by a Coulomb-type friction law (Boulton and Jones, 1979; Clarke, 2005) as follows:

where τ_u is the ultimate shear strength, c is the cohesion parameter, N is the normal effective stress and φ is the friction angle. This kind of law, where shear stress is independent of the sliding velocity, allows for unstable behaviour leading to failure. In a general case, glaciers remain stable because till and water pressure are not equally distributed at the basal interface, leading to sticky spots where stress concentrates to balance gravitational forces with lateral drag at glacier margins (Cuffey and Paterson, 2010). The Aru collapses would thus be an example where basal shear stress becomes limited to till strength in such large zones that resisting forces can no longer balance gravity, eventually leading to catastrophic failure. The latter would happen for a bedrock with low roughness, which provides less vestigial resistance (surface-normal stresses) to constrain ice velocity in case of failure in the till. We propose that the change in effective normal stress due to increasing pore pressure in the till beneath the temperate zones of the Aru glaciers quickly reduced the ultimate strength of glacier till and limited the basal shear stress to a maximum value lower than the driving stress. The glacier shape could not adjust fast enough to reduce the driving stress due to strain-rate limitation in the cold-based zones, leading then to the accumulation of large stresses in the remaining sticky spots (Fig. 8) until their sudden rupture. The sticky spots were likely remnants of frozen-stiff till rather than solid rock irregularities, rendering them susceptible to failure under high stress and vulnerable to thaw from water-saturated temperate surroundings and increasing deformational heat. The high clay and silt content measured in the till is indicative of lithologies having unique low friction-angle properties (Iverson et al., 1998) and higher sensitivity of the shear strength to changes in water pressure.

One likely scenario for the development of the now collapsed Aru glaciers is that they grew in the past (pre-industrial climate) in colder conditions with low melting rates in summer, allowing for rigidity in the basal till to support high driving stress (see Sect. 4.3). The low water pressure meant that there was likely very little sliding and therefore very little production of till at that time. Upon commencement of some sliding, which may have occurred gradually over an increasing area of the bed during the past century, till production increased and the local glacier deformation regime tended to adapt to the distribution of till and liquid water reaching the bed. At this stage, percolation into the glaciers and accumulation of meltwater beneath the glaciers increased so rapidly in recent years (Kääb et al., 2018) that the glaciers could not keep balance with the changing conditions at the bed. Sliding may also have contributed to an increase in the water pressure in a positive feedback by destroying any efficient drainage system (Clarke et al., 1984). The contributory role of soft-bed lithology in the collapses is therefore likely threefold by (i) deforming with a plastic rheology when shear strength is reached, (ii) providing for low roughness at the ice-bed contact, and (iii) maintaining high water pressure, while sliding speed increases; a known process that accounts for surge behaviour (Clarke et al., 1984; Raymond, 1987; Fowler et al., 2001). High content in clay and silt probably also leads to low hydraulic conductivity favourable to higher water pressure in the till (Fowler et al., 2001).

We suggest that the existence of a high friction area under the Aru-2 glacier tongue prior to collapse is due to both higher basal normal stresses (Fig. S10), which increased the till strength, and higher lateral drag along the west side of the detachment which decreased the basal shear stresses compared to the Aru-1 glacier (Fig. 9). In this way, and contrary to the Aru-1 glacier, basal shear stress under the tongue of the Aru-2 glacier only approached the ultimate shear strength of the till just before the final collapse in response to both decreasing resistance by the lateral margin (due to crevassing) and increasing driving stresses (due to bulging).

5.4 Till-strength controlled glacier collapses

The Aru collapses, and in retrospect the Kolka Glacier collapse, define a newly recognized type of avalanching glaciers characterized by an underlying failing substrate. These ice slides could occur on glaciers with fairly low angles, therefore involving potentially large volumes of material and presenting serious consequences in terms of hazard potentials. The high sensitivity of the ultimate shear strength of the substrate to pore water pressure, combined with low bed roughness, allows for a dramatic and sustained change in basal friction conditions capable of driving this kind of instability.

The Kolka event in 2002 in the Caucasus Mountains is probably another example of this type of instability in which the maximum shear strength of the till is exceeded by a sudden increase in basal shear stress at constant effective normal stress. Indeed, during the few weeks before this collapse, significant mass was added on top of the glacier by rock and ice-fall activity, increasing basal shear and normal stresses (Huggel et al., 2005; Evans et al., 2009). This could change the surface slope, which reached the till friction angle and triggered the failure. If the till was saturated with water and had low hydraulic conductivity, increasing pore water pressure could also have compensated for the rising normal stress, keeping the normal effective stress constant. Once maximum shear strength is exceeded, failure is triggered (Evans et al., 2009). This hypothesis is plausible since Kolka Glacier is known to be a surging glacier, able to store large amounts of liquid water, and high water content and pressure were observed before its 2002 collapse (e.g. unusual ponds observed on the glacier prior to collapse) (Kotlyakov et al., 2004).

However, changes in till strength in response to changing in water pressure are likely also involved in the surge mechanisms of temperate glaciers without the large majority of surges turning into gigantic collapses. This renders sudden changes in till strength as a necessary but not fully sufficient condition, for collapses controlled by till strength. The necessary secondary condition for catastrophic failure is a sustained high driving stress with low bed roughness, coincident with weakening till. This means that the glacier has to grow over time, atop of a more stable substrate capable of supporting higher driving stresses. In particular, freezing conditions allow for the development of relatively thick glaciers on slopes that would otherwise be unable to support such high shear stress under the presence of liquid water. This makes the spatio-temporal interplay of soft-bed characteristics and the polythermal glacier regime a prerequisite of the Aru collapses, whereas for the Kolka Glacier the additional loading over a short time should have caused a fast increase in shear stress, significantly exceeding the glacier's normal conditions.

In many of the world's glacierized regions, ongoing atmospheric warming increases surface melt and the amount of water reaching glacier beds, thereby modifying the till shear strength. This development is therefore in theory capable of driving more till-strength-controlled instabilities and collapses. The most impacted glaciers would be those flowing on soft and highly erodible bed lithologies at high driving stress, particularly those with heterogeneous thermal structure (polythermal glaciers). Such glaciers are mostly localized in cold and dry climates, where a small increase in temperature results in a relatively large change in melting conditions such that the amount of water reaching the glacier base can significantly increase instability. In reality, however, an array of factors and their specific (and to this point rare) interplay in time and space are necessary to catalyse glacier collapsing as observed for the Aru and Kolka glaciers.