Basal stress controls ice-flow variability during a surge cycle of Hagen Bræ, Greenland

Abstract

Basal conditions play an essential role in the dynamics of outlet glaciers, but direct observations at the bed of glaciers are challenging to obtain. Instead, inverse methods can be used to infer basal parameters from surface observations. Here, we use a simple ice-flow model as a forward model in an inversion scheme to retrieve the spatio-temporally variable basal stress parameter for Hagen Bræ, North Greenland, from 1990 to 2020. Hagen Bræ is a surge-type glacier with up to an order of magnitude variability of winter velocities near the grounding line. We find that downstream changes in the basal stress parameter can explain most of the variation of flow velocity, and we further identify a region of high resistance ~20–40 km from the grounding line. We hypothesise that this region of high resistance plays an important role in controlling glacier discharge.

Introduction

Surge-type glaciers are characterised by fluctuations between long periods of slow ice flow (quiescent phase) and short periods of fast ice flow (surge phase) (Meier and Post, 1969). This behaviour has been observed in only 1% of the world's glaciers (Jiskoot and others, 1998; Sevestre and Benn, 2015). A surge causes a redistribution of mass and may as a consequence enhance ice loss as well as modify glacial run-off (e.g. Hewitt and Liu, 2013). In recent years, numerous surges have been documented through the use of remote sensing (e.g. Dunse and others, 2015; Steiner and others, 2018; Rashid and others, 2019; Bhambri and others, 2020; Paul, 2020; Solgaard and others, 2020).

A surge is commonly triggered by an imbalance in mass flux: during quiescence, mass gain exceeds mass loss and thus mass is accumulated, leading to a steepening of the surface slope. This causes the driving stress to exceed the resistive stress and consequently, the glacier speeds up – typically by at least an order of magnitude compared to the velocity during quiescence. This increase in velocities creates an ‘ice bulge’ that moves downstream as a kinematic wave (e.g. Kamb and others, 1985; Mayer and others, 2011; Adhikari et al., 2017). Two main hypotheses have been put forward to explain the underlying mechanisms of glacier surges: Clarke (1976) proposed that surges are regulated by thermal processes, i.e. a change in temperature of the basal ice of cold glaciers, whereas a later study by Kamb (1987) argued that surges are controlled by changes in basal water pressure leading to a decrease in basal friction. Sevestre and Benn (2015) presented a unified model based on coupled mass and enthalpy budgets. Benn and others (2019) applied this model quantitatively and predicted that enthalpy at the bed increases during the quiescent phase as heating exceeds enthalpy losses. Regardless of the specific surge mechanism, it is clear that the basal conditions of a surge glacier are an important factor in modulating the ice flow (Benn and others, 2019).

To get insight into which processes control the surge cycle of a glacier, it is important to gain information on its basal regime/characteristics. However, direct observations of basal conditions are difficult to obtain and only a few direct measurements exist (e.g. Lappegard and others, 2006; Ryser and others, 2014; Lefeuvre and others, 2015). Instead, information on changes in the basal system are typically inferred from surface observations – especially, surface velocity measurements or observations of surface topography and slope (e.g. Fatland and Lingle, 2002; Johnson and others, 2013; Haga and others, 2020). Another common approach for inferring basal conditions is the application of inverse methods, although this method has only seen limited applications to surge-type glaciers (e.g. Jay-Allemand and others, 2011; Beaud and others, 2021). An inverse method involves tuning the input parameters of ain ice-flow model to obtain the best fit to the observations. MacAyeal (1992, 1993) was the first to apply inverse techniques to ice-flow models in order to infer basal stress parameters. Subsequent studies have used variations of the inverse method, and recent studies have expanded on the use of inverse methods; for example, Ranganathan and others (2020) present a method to invert for the basal stress parameter and ice rheology simultaneously. Many of the studies employing inverse methods have focused on the dynamics of fast-flowing ice streams, aiming to retrieve information on basal shear stress (Joughin and others, 2001; Sergienko and others, 2008; Gillet-Chaulet and others, 2016). In such studies, the basal shear stress is typically assumed to be linearly related to the basal velocity, and results often show significant spatial variability in basal shear stress of several orders of magnitude between the front of the glacier and upstream areas (Joughin and others, 2001; Sergienko and others, 2008). Inverse methods have also revealed that the bed provides practically no frictional resistance to the fast flow for three of Greenland's largest marine-terminating outlet glaciers (Shapero and others, 2015). One of the few studies to invert for basal conditions of a surge-type glacier is the study by Jay-Allemand and others (2011) who retrieve the basal stress parameter during a surge cycle of Variegated Glacier. They interpret the changes in flow resistance to be caused by the variability of basal water volumes and hence changes in basal water pressure.

Here, we investigate Hagen Bræ, a surge-type glacier situated in North Greenland (81° 20′57" N, 28° 22′25″ W, depicted in Fig. 1a) where <5% of surge glaciers worldwide can be found (Sevestre and Benn, 2015). In contrast to the well-studied surge-type glaciers found in Alaska and the Himalayas (e.g. Jay-Allemand and others, 2011; Paul, 2020), Hagen Bræ is a marine-terminating glacier resting on a bed down to >400 m below sea level. From the grounding line and up to 70 km upstream, the bed is below sea level, and the first 20 km of the bed is retrograde (Fig. 1b). Hagen Bræ is 75 km long and 10 km wide and is one of North Greenland's major outlet glaciers. In the early 2000s, the glacier underwent a surge with variations in flow velocity of up to an order of magnitude. This was documented by Solgaard and others (2020) who proposed a periodicity of ~20–30 years.

Fig. 1. Overview of the data used. (a) Overview of Hagen Bræ (3rd August 2020. Contains modified Copernicus Sentinel-2 data 2020, processed by ESA). The extent of the image is shown in red on the inset of Greenland. The flight line is shown in red, while a flow line is shown in blue. The velocity transects show the lines along which the velocities are shown in Figure 7. The lightest green corresponds to 10 km from the downstream boundary, with the progressively darker each 10 km further upstream. The orange lines indicate the 1996 grounding line region (ESA Greenland Ice Sheet CCI, Grounding Lines from SAR Interferometry). The crosses denote the boundaries of the 55 km modelling domain. (b) A sample of the topography data used. The data are based on AeroDEM 1978, the 2002 ASTER DEM and three PROMICE flight campaigns (2007, 2011 and 2015) with linear interpolation in time to get a profile for each inversion year. All profiles are either measured along or nearest neighbour interpolated to the 2007 PROMICE flight line shown in red in (a). (c) Winter velocity datasets interpolated to the 2007 PROMICE flight line.

In this study, we investigate the variations in basal conditions of Hagen Bræ from 1991 to 2020. In contrast to the large ice streams investigated by Shapero and others (2015), we expect Hagen Bræ to be subject to frictional resistance at the bed, especially during its quiescent phase. We use the common relationship between the basal stress and basal velocity as suggested by Weertman (1957), where the basal shear stress relates to a power of the basal velocity with a factor denoted by the basal stress parameter, C. The basal stress parameter is a measure of resistance to sliding at the ice–bed interface, and our aim is to retrieve information on the basal shear stress and the basal drag, to document how basal resistance changed during and after the 2002 surge.

Data

We use observations of surface and bedrock topography and surface velocities as input to our forward model (see Fig. 1). As the focus in this paper is on annual and decadal-scale variation in ice flow and basal properties, we only use winter velocities. Including summer velocities would clutter the variations with seasonal change, and the noisier summer datasets would introduce instabilities in the forward model and inversion. We obtain surface topography from two digital elevation models (DEMs) and measurements along the glacier centre line. The oldest dataset, generated from oblique aerial photographs, is the AeroDEM from 1978 (Korsgaard and others, 2016a), with a reported vertical accuracy of 6 m. The second DEM is derived from two Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) scenes acquired on 30 June 2002. The DEM was extracted using MMASTER (Girod and others, 2017; Hugonnet and others, 2021), and bias corrected using the ArcticDEM mosaic (Porter and others, 2018).

Finally, surface topography was measured along the glacier centre line in August 2007, 2011 and 2015 during the Programme for Monitoring the Greenland Ice Sheet (PROMICE) flight campaigns. Reported uncertainties are 0.1 m for surface topography measured with a near-infrared laser scanner (Sørensen and others, 2018). Linear interpolation in time is applied to obtain a surface profile for all years. The bedrock topography was measured during the PROMICE 2007 flight campaign. Reported uncertainties are 35 m for the bedrock topography measured with ice-penetrating radar (Sørensen and others, 2018).

Ice velocities are based on synthetic aperture radar offset tracking products from various satellite campaigns (see Table 1 and Solgaard and others (2020) for details). The data sets span the period 1991/92 to 2019/20. The winter velocity datasets are combined from multiple image pairs spanning shorter ranges. Reported ice velocity uncertainties for the source datasets range from 10 to 35 m a⁻¹. The errors for the winter datasets are estimated from the ensemble of image pairs. We include the RMS of estimated uncertainty for each dataset in Table 1 to indicate typical values. The spatially varying errors are used in the inversion.

Table 1. All winter velocity datasets were combined from datasets spanning shorter ranges

The acquisition range of each winter dataset is listed in the table. The source datasets are referenced below. RMS(σ_obs) is the RMS of estimated observation error along the flight line for each winter dataset. Errors are estimated from the local standard deviation of the surface velocity maps generated by the offset tracking included in each velocity product.

References. a: Data from ESA Ice Sheets Greenland CCI project (www.esa-icesheets-greenland-cci.org) and Nagler and others (2015). b: Data from the PROMICE project (Solgaard and Kusk, 2021).

The 2-D velocity datasets are linearly interpolated onto the 2007 flight line (Fig. 1). In Figure 1b we see that the flow velocity is temporally and spatially variable. Approximately 35–55 km from the grounding line, the velocities are relatively similar, within a range of 200–300 m a⁻¹. Towards the grounding line, we see large variations in the different datasets, from nearly stagnant flow in the 90 s to velocities ~1000 m a⁻¹ at the peak of the surge in 2004/05.

Methods

An inverse problem is finding the set of model parameters that yields the best match between forward model output and observations. Here, the target model parameter is the basal stress parameter, the forward model is an ice-flow model, and the model output and observations are surface velocities. For a given input of basal stress parameters, a velocity field is returned, which can be compared to observations. The goal of the inversion is to find the set of basal stress parameters for each dataset that minimises the squared residual between the forward model and the observed velocities. The basal stress parameter is discretised spatially to retrieve the spatial variation, and the inversions are run independently for each winter velocity dataset to retrieve the temporal variation. Below, we describe the forward model and our inverse approach in more detail, and we motivate our use of a regularisation parameter.

Forward model

Weertman (1957) proposed that the relation between the basal shear stress, τ_b, and basal velocity, u _b, can be represented as a power law

(1)$$\tau_{\rm b} = - C u_{\rm b}^{1/m},\; $$

where m is the power-law exponent, and the basal stress parameter, C, is a positive scalar representative of local properties of the bed (roughness, water pressure, nature of the bed, etc.), assumed here independent of u and τ_b (cf. Gillet-Chaulet and others, 2016). Maier and others (2020) find this sliding relation to be applicable to most regions in Greenland. From Eqn (1) it is seen that a high C indicates high resistance to flow. In this study, C is assumed to vary spatially and temporally and is the target parameter for inversion.

This description of basal properties is implemented in an ice-flow model. Here, we use a 1-D formulation of the shallow shelf approximation (SSA) (MacAyeal, 1989). It is a vertically integrated flow model neglecting vertical shear and is most applicable to ice streams where ice motion is dominated by basal sliding (Joughin and others, 2001, Jouvet, 2014). This implies that we assume that the ice is flowing as plug flow, i.e. that observed surface velocities equal basal velocities. The SSA 1-D vertically integrated flow line stress balance (Tsai and others, 2015) with an added parameterisation of lateral drag (Van Der Veen and others, 2011) is

(2)$$\matrix{ & \underbrace{2 A^{-1/n} \left(H u_x\left\vert {u_x} \right\vert ^{1/n-1} \right)_x}_{\hbox{membrane}} - \underbrace{C u^{1/m}}_{\hbox{basal}} \cr & \quad- \underbrace{2Hw^{-1}\left(uA_\ast ^{-1}w^{-1}\right)^{1/n}}_{\hbox{lateral}} - \underbrace{\rho g H s_x}_{\hbox{driving}} = 0,\; }$$

where the derivative notation $u_x = {\partial {u( x) }\over \partial {x}}$ is used, A is the depth-averaged temperature-dependent rheological coefficient in Glen's flow law (also referred to as the creep parameter), n is the corresponding exponent, H is the ice thickness, u is the velocity, C is the basal stress parameter, w is the width of the glacier, $A_\ast$ is an increased A used in the lateral drag parameterisation to account for weakening of shear margins, ρ is the density of ice, g is the gravitational acceleration and s is the surface elevation. The four terms in Eqn (2) are membrane stresses held by viscous deformation, basal shear stress held at the base by till strength, lateral drag from the sides of the glacier and driving stress (Bueler and Brown, 2009). The following common values are assumed: A = 9.3 × 10⁻²⁵ s⁻¹ Pa⁻³, n = 3, m = 3, ρ = 900 kg m⁻³ and $A_\ast = 2A$ (Borstad and others, 2013; Bondzio and others, 2016). The value for the creep parameter, A, corresponds to a temperature of − 5°C (Cuffey and Paterson, 2010). The sensitivity of the results to the creep parameter, A, and power-law exponent m will be investigated. Although the value of $A_\ast$ has an effect on the absolute values of basal stress, it has almost no effect on the pattern of spatial variability, so the sensitivity of $A_\ast$ will not be further discussed.

The forward model (Eqn (2)) is solved numerically for velocity using a central difference scheme on staggered grids, with u and C staggered half a grid spacing from t and s _x. The non-linearity of the system is solved by applying a Picard iterative scheme (Langtangen, 2003, ch. 4), linearising the viscosity at each iteration. The velocity domain is discretised on N = 401 gridpoints yielding a spatial resolution of ≃ 140 m. The domain boundaries are slightly upstream of the grounding line and 55 km upstream (Fig. 1a). The focus in this study is not on grounding line dynamics, but rather the upstream regions and the changes in basal stress, and thus we apply Dirichlet boundary conditions on both boundaries.

Balise and Raymond (1985) find that spatial variation in basal friction over scales less than ice thickness is unlikely to have a significant effect on observable surface velocities. C is therefore discretised evenly on M = 15 gridpoints, corresponding to a spatial resolution of 3.9 km or ~ 3.9 times the maximum ice thickness of Hagen Bræ. Linear interpolation is used between these points to get a value for C at each velocity gridpoint, as required in the discretisation of Eqn (2).

The inverse problem

The inverse problem is to determine the basal stress parameter C that minimises the sum of squared differences between the observed velocities and forward model prediction

(3)$$J_0 = \sum_{i = 1}^{N}\left({{u}_i -{\,f_ i}( {\bf c}) \over \sigma_i}\right)^{2},\; $$

where N is the number of gridpoints, u _i are observed velocities, f _i(c) are predictions by the forward model, c is a vector with the discretised C and σ_i are reported errors of u _i. The minimisation of J ₀ alone is an ill-posed problem, i.e. there are many different solutions. To condition the problem and focus on larger scale variability, the spatial variation of C is restricted by introducing a regularisation term. Regularisation is a classic tool to stabilise the inversion, ensure the existence of a global minimum and reduce the risk of overfitting to noise and modelisation errors (Aster et al., 2013, ch. 4). First-order Tikhonov regularisation is used, minimising short wavelength variations of C, similar to Jay-Allemand and others (2011) and Gillet-Chaulet and others (2012). The regularisation term is

(4)$$J_{\hbox{reg}} = \sum_{i = 1}^{M-1}( c_{i + 1}-c_i) ^{2},\; $$

where M is the number of gridpoints for C, and c _i is the value of C at gridpoint i. The resulting cost function subject to minimisation is defined as the sum of data misfit and regularisation cost functions

(5)$$J = J_0 + \alpha J_{\hbox{reg}},\; $$

where the regularisation parameter α ≥ 0 controls the relative weight of J ₀ and J _reg under minimisation. A higher α preferentially weights the minimisation of J _reg, leading to a smoother inversion of C. The minimisation of this cost function is a well-posed problem, and is solved using Gauss–Newton iteration (Aster et al., 2013, ch. 10)

(6)$$\matrix{ {\bf c} & = {\bf c}_0 + \left[{\bf J}( {\bf c}_0) ^{\top} {\bf W} {\bf J}( {\bf c}_0) + \alpha {\bf L}^{\top} {\bf L}\right]^{-1} \cr & \quad\times \left[{\bf J}( {\bf c}_0) ^{\top} {\bf W}\left({\bf u} - {\bf f}( {\bf c}_0) \right)- \alpha {\bf L}^{\top} {\bf L} {\bf c}_0\right],\; }$$

where c and c₀ are the next and current vector of C, J is the Jacobian matrix whose elements $J_{ij} = {\partial {f_i( {\bf c}) }\over \partial {c_j}}$ are approximated as finite forward differences, W is a diagonal weighting matrix with elements $W_{ii} = \sigma _i^{-2}$, L is the first-order Tikhonov roughening matrix, u are observed velocities and f are the forward model outputs. The iteration is stopped when the absolute relative change in J from one iteration to the next ≤ 10⁻⁵, requiring 7–8 iterations for the majority of velocity datasets. This convergence threshold was determined by numerical experiments, and further iterations beyond yielded negligible changes in C.

Regularisation parameter

The regularisation parameter α is determined using L-curve analysis (Hansen, 1992). For each dataset, the inversion is run for a range of α-values, and log J ₀ vs log J _reg is plotted. We choose the α of maximum curvature. The rationale behind this choice is that at the point of maximum curvature, by further decreasing α, a small decrease in J ₀ can only be achieved through a large increase in J _reg. In other words, a small decrease in data misfit can only be achieved by increasing the spatial variation of C.

The L-curve analysis is shown in Figure 2, and is based on inversion with 29 values of α from 10^−6.5 to 10⁻³ Pa⁻² m^2/3 s^−2/3 (in the following, we omit units of α for brevity). Most datasets exhibit a broad region of maximum curvature, which defines the α that optimises the trade-off between the misfit norm J ₀ (Eqn (3)) and the regularisation norm J _reg (Eqn (4)). The success of the L-curve approach is dependent on a sharp corner feature, clearly marking the α of optimal trade-off. None of the L-curves in Figure 2 show a clearly determined corner. For some inversions, the maximum curvature yields inversions with a highly spatially variable C, and a higher α is chosen to yield spatial variations on the same scale as the more well-determined α's. The determined α's vary from 1.3 × 10⁻⁵ to 1.3 × 10⁻⁴, and are weakly negatively correlated with measurement errors, reflecting that higher uncertainties lead to lower values of J ₀ which are compensated by reducing α (Fig. 10). The α's sit at J ₀ between 2 × 10³ and 2 × 10⁵, and at J _reg between 6 × 10¹² and 2 × 10¹³.

Fig. 2. L-curve analysis for all datasets. Black circles indicate the used α. The associated α-values are denoted in the legend.

Results

The inversion results for the basal stress parameter C are shown in Figure 3a. There is clear spatial and temporal variation in C. We divide the results spatially into three regions: the upstream region from 40 to 55 km, the middle region from 20 to 40 km the downstream region from 0 to 20 km. Furthermore, we divide the results temporally into three different periods: the high resistance period during 1991–96 (blues), the low resistance period during 2002–16 and 2018–20 (yellow-purples) and the intermediate resistance period during 2016–18 (greens).

Fig. 3. Results for optimal α determined by the L-curve analysis in Figure 2. (a) Basal stress parameter. (b) Basal stress. (c) Data-model misfits.

In the upstream region, C increases towards the grounding line for all years, from 1–2 mPa m^−1/3 s^1/3 to 5–6 mPa m^−1/3 s^1/3. In the middle region, there is a region of high resistance, with the low resistance period peaking at ~25 km from the grounding line with C ≃ 6–7 mPa m^−1/3 s^1/3, and the intermediate and high resistance period peaking ~35 km from the grounding line with $C\simeq$4.5–6.5 mPa m^−1/3 s^1/3. In the downstream region, there is high variability in C. For the low resistance period, there is a decrease to $C\simeq$1.5–2.5 mPa m^−1/3 s^1/3. The intermediate resistance period decreases towards the grounding line towards the levels of the low resistance period but remains ~1 mPa m^−1/3 s^1/3 higher than the median of the intermediate resistance period. In the high resistance period, there are two distinct patterns; the years 1992/93 and 1993/94 remain around the peak value of ≃ 6 mPa m^−1/3 s^1/3, with a slight decrease from the peak and a slight increase up to the grounding line region. The years 1991/92 and 1995/96 increase to 8 mPa m^−1/3 s^1/3 ~25 km from the grounding line, decrease to 5.5 mPa m^−1/3 s^1/3 ~20 km from the grounding line, and increase to 8 mPa m^−1/3 s^1/3 towards the grounding line region.

The basal stress, τ_b, follows a similar pattern as C. In the upstream region, stress is increasing in all years from 20–50 to 110–130 kPa at ~40 km from the grounding line. In the middle region, the basal stresses span 70–130 kPa, with the majority of years having values of 100–110 kPa. In the downstream region, most years have decreasing stress, with values spanning 30–80 kPa near the grounding line. Up to the grounding line region, the highest stresses are seen for the high stress period, with values ~60–80 kPa. For the whole study period, there is a ‘sticky’ region of high resistance present ~20–40 km from the grounding line.

In Figure 1b, we see that these three resistance periods correspond to variations in the ice velocities. The high resistance period matches the decreased velocity from 200 to 250 m a⁻¹ at 40–55 km down to almost stagnant at the grounding line. The low resistance period shows the highest velocities, with a steady increase from 250 to 400 m a⁻¹ at 40 km, to 400–1000 m a⁻¹ at the grounding line. The intermediate resistance period shows velocities which are a combination of the high resistance and low resistance periods. In the upstream and middle regions, the intermediate resistance velocities are similar to the high resistance velocities. In the downstream region, the intermediate resistance velocities are similar to the low resistance velocities.

The data-model misfits are shown in Figure 3c. The majority of misfits are within 20 m a⁻¹. Most of the larger misfits occur in the same regions, and in particular, the model is consistently underestimating velocities ~40, 20 and 12 km, and overestimating velocities ~50, 30 and 18 km. The highest misfit values are seen for the years 2002/03, 2004/05, 2005/06 and 2006/07, where we also see the highest velocities. Our results are insensitive to further improvements in spatial resolution.

Discussion

In the following, we discuss the impact of the regularisation parameter on our results, the limitations of our forward model, and finally, how our results may be interpreted in the context of ice-flow dynamics.

Regularisation parameter

From a probabilistic view, the regularisation imposes a prior of a smooth solution to the inversion. The value of α determines the weighting of this prior and thus has a significant influence on the resulting C distribution. To investigate the robustness of C for different α, we further investigate results from 1991/92 and 2006/07. In Figure 2 we see that these two L-curves represent two outlying cases; the 1991/92 curve shows no corner features and is almost linear, whereas the 2006/07 curve shows a more well-defined corner region. These two cases represent a poorly determined and better determined α, respectively. The points along the L-curve are spaced further apart around the chosen α in J ₀–J _reg-space for 1991/92, compared to the results for 2006/07. This indicates that the 2006/07 inversions are relatively similar for all the α's in the corner region, compared to 1991/92 where the inversions around the chosen α vary more. In Figure 4 the α determined from the L-curve has been scaled by 3 and 1/3 (equivalent to four steps along the L-curve in either direction) and the effect on the resulting C, basal stress and residuals are shown.

Fig. 4. Robustness of basal stress parameter results for a range of regularisation parameters, α. The main lines represent the result for α based on the L-curve analysis. The dotted and dashed lines represent less fierce and more fierce regularisation, respectively. The area between the dotted and dashed lines has been shaded. (a) Basal stress parameter. (b) Basal stress. (c) Data-model residual.

The rationale behind the L-curve analysis is that by choosing the point of maximum curvature, a decrease in α will lead to a small decrease in misfit cost for a large increase in regularisation cost. Equivalently, an increase in α will lead to a large increase in misfit cost for a small decrease in regularisation cost. In Figure 4 we see that this is somewhat true for 2006/07. Decreasing α (dotted lines) leads to a less smooth resolution in C for a small decrease in misfit. Increasing α (dashed lines), gives a slightly smoother solution but with a notable increase in data misfit. This reflects the corner feature of the 2013/14 L-curve analysis and supports the choice of α for this dataset. For 1991/92 a decrease in α (dotted) allows for a less smooth solution with a notable decrease in misfit. A smoother solution is returned for an increase in α (dashed), with a notable decrease in misfit. Thus for 1991/92, the results do not align well with the rationale of applying the L-curve, making the choice of α difficult. To ensure that the conclusions drawn from C are actual variations at the bed, results for more fiercely regularised inversions than the L-curve determines are shown in Figure 9. The analysis regarding the distribution is still valid in this figure, and thus we find that the decrease in C during a surge is a real temporal variation.

We expect to find a negative correlation between α and observation error estimates for the datasets. A decreased observation error leads to a larger data misfit in the cost function and must be balanced by an increased weight on regularisation to yield similar inversions. This is found to be the case, shown in Figure 10 where we compare α to the RMS of velocity observation error estimates.

Forward model

By definition, the SSA neglects vertical shearing. This assumption is valid for ice flow where the horizontal motion is dominated by basal ice motion, but this is unlikely to be the case for upstream regions of the Greenland glaciers. Furthermore, we emphasise that this study's results attribute variation in the ice flow from 1991 to 2020 to changes in basal conditions along with topography changes. Thus, it is important to highlight that the results in this study are consistent with the dynamics of the SSA and the applied surface observations, but do not necessarily reflect actual physical till properties (Habermann and others, 2013). Even so, we argue that basal changes dominate the temporal variation in flow velocity, and thus the model is adequate for understanding temporal changes in flow resistance at the base. It is noted that while the 1-D SSA is valid for flow lines, we use velocities from the flight line that differs from a flow line (see Fig. 1a). From the velocity transects (see Fig. 7), we find that up to 55 km from the grounding lines, the velocity variations across the centre of the flow are negligible, supporting the use of flight line velocities.

The forward model is sensitive to changes in the physical input parameters. Gillet-Chaulet and others (2016) performed a parameter study with a similar friction law for Pine Island Glacier (Antarctica), and found that variation in m (power-law exponent) yields very different forward model fits. For C values constant in time, they found that linear and slightly non-linear friction laws (m = 1,  3) yield a poorer fit to data than a highly non-linear friction law, with the best fit to data for m ≥ 5. In this study, we have used m = 3. Our findings’ robustness has been tested for values of m = 1 and m = 5, with the results presented in Figure 5. We find that our analysis is robust to different degrees of weak non-linearity, with negligible difference in data-misfit. The presence of the sticky region, however, is less pronounced when applying the linear friction law.

Fig. 5. Sensitivity analysis for m, using the datasets 1992/93 and 2013/14 for m = 1,  3 and 5. (a) Basal stress parameter. Results for each m are normalised to the grounding line value of C for the corresponding 1992/93. (b) Basal stress. (c) Data-model residual.

The creep parameter A is assumed constant in this study, and the results are obtained using A = 9.3 × 10⁻²⁵ s⁻¹ Pa⁻³, corresponding to a temperature of − 5°C. The sensitivity of our results to changes in this parameter has been tested with values of A corresponding to temperatures of − 20 and 0°C, corresponding to the lower limit (approximate average winter surface air temperature in the region) and upper limit. The results of this are presented in Figure 6. It is seen that for both the datasets, a lower resistance (both in terms of basal stress parameter and basal stress) is found for colder ice. In particular, the location of the highest resistance for 2013/14 is ~40 km for the colder ice and 30 km for the warmer ice. For 2013/14, the warmer model is less capable of fitting observations in the downstream region.

Fig. 6. Sensitivity analysis for A, using the datasets 1992/93 and 2013/14. The model has been run for A = 1.2 × 10⁻²⁵ s⁻¹ Pa⁻³ (corresponding to − 20°C), A = 9.3 × 10⁻²⁵  s⁻¹ Pa⁻³ (corresponding to − 5°C) and A = 2.4 × 10⁻²⁴ s⁻¹ Pa⁻³ (corresponding to 0°C). − 20°C is around the winter average surface temperature for glaciers in the region. (a) Basal stress parameter. (b) Basal stress. (c) Data-model residual.

Fig. 7. Transverse velocities, corresponding to the transverses shown in Figure 1a, where darker lines are closer to the grounding line. Velocities are normalised to the flight line velocity. Positive distance from flight line is towards southeast.

Basal stress parameter

Our results show that the surge of Hagen Bræ coincides with marked changes in the basal stress parameter. The changes in the basal stress parameter also take place in our results with a constant topography (see Fig. 8). Thus, we hypothesise that the surge is not driven by topography alone but is rather accompanied or triggered by changes in basal stress. The most likely mechanism for changing basal conditions is variations in basal water pressure at the bed of the glacier. This is in line with Solgaard and others (2020) who find that the surge of Hagen Bræ was triggered by changes in englacially stored meltwater. We investigate this by comparing the results from winters 2016/17 and 2017/18 with surface meltwater production (see Fig. 11 (Mankoff and others, 2020)). Winters 2016/17 and 2017/18 exhibit lower velocities than the preceding and following winters but do not have a notably lower surface meltwater production. Thus, although the surface meltwater likely modulates basal processes there is not a direct one-to-one relationship between the amount of meltwater input during summer and the magnitude of the velocities the following winter. Thus, we hypothesise that in addition to surface meltwater volumes, the basal sliding is not only governed by local processes, but is also influenced by factors that vary over larger spatio-temporal scales. This has been suggested for other surge-type glaciers, for example, Vallot and others (2017) found that calving affects buttressing and thereby the glacier flow field, although for Hagen Bræ specifically, we suggest that a delayed release of stored englacial water is more likely (cf. Lingle and Fatland, 2003). Based on surface velocity data, the Hagen Bræ surge peaks in 2004/05. Although the velocities are higher for almost the entire domain, the main differences in C are seen in the downstream region. This could indicate that while the basal stress of the whole base of the glacier is an important dynamic factor for its flow velocity, the basal stress decrease of the downstream region enhances the flow in the entirety of the domain. Our results indicate that the release of meltwater starts affecting the basal conditions ~35 km from the grounding line, where we see the onset of temporal variations in C. Habermann and others (2013) find similar structures for Sermeq Kujalleq (Jakobshavn Isbræ), where the basal stress parameter builds up to a similar sticky region, after which the different velocity regimes show very different basal stresses. We hypothesise that variations in the available amount of water at the sticky region could have a disproportionately large impact on the glacier flow of Hagen Bræ.

Fig. 8. Results for constant topography (2007 PROMICE). (a) Spatial distribution of C for each dataset. (b) Basal stress for each dataset. (c) Data-model misfits for each dataset.

Fig. 9. Results for optimal α multiplied by 6. (a) Spatial distribution of C for each dataset. (b) Basal stress for each dataset. (c) Data-model misfits for each dataset.

Fig. 10. α plotted against the RMS of reported observational velocity uncertainties on log–log scale.

Fig. 11. Freshwater run-off for Hagen Bræ from Mankoff and others (2020) for the years 2015–19. Output from two regional climate models are displayed: MAR and RACMO.

Conclusions

We use an inverse approach to investigate basal conditions of Hagen Bræ, a surge-type glacier, over a surge cycle. We find that the onset of the early 2000s surge is not driven by driving stress alone, but is related to changes in basal conditions ~0–30 km from the grounding line. Around the grounding line, the basal stress parameter C ranges from 6 to 8 mPa m^−1/3 s^1/3 for the quiescent period and 1.5 to 2.5 mPa m^−1/3 s^1/3 for the following surge period. The basal stress around the grounding line is 55–65 kPa for the quiescent period (1991–96), ~75 kPa at the peak of the surge (2002/03) and 30–50 kPa for the remainder of the studied period (2004–20). Our results are insensitive to further increases in spatial resolution, and our conclusions are robust for various degrees of non-linearity in the basal sliding law. We hypothesise that the flow of the glacier is sensitive to changes in basal conditions, particularly from the downstream end of the sticky region ~20–40 km from the grounding line, where the basal stress parameter is highly variable between the quiescent and surge periods.