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, Reference Meier and Post1969). This behaviour has been observed in only 1% of the world's glaciers (Jiskoot and others, Reference Jiskoot, Boyle and Murray1998; Sevestre and Benn, Reference Sevestre and Benn2015). 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, Reference Hewitt and Liu2013). In recent years, numerous surges have been documented through the use of remote sensing (e.g. Dunse and others, Reference Dunse2015; Steiner and others, Reference Steiner, Kraaijenbrink, Jiduc and Immerzeel2018; Rashid and others, Reference Rashid, Majeed, Jan and Glasser2019; Bhambri and others, Reference Bhambri2020; Paul, Reference Paul2020; Solgaard and others, Reference Solgaard2020).
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, Reference Kamb1985; Mayer and others, Reference Mayer, Fowler, Lambrecht and Scharrer2011; Adhikari et al., Reference Adhikari, Ivins and Larour2017). Two main hypotheses have been put forward to explain the underlying mechanisms of glacier surges: Clarke (Reference Clarke1976) 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 (Reference Kamb1987) argued that surges are controlled by changes in basal water pressure leading to a decrease in basal friction. Sevestre and Benn (Reference Sevestre and Benn2015) presented a unified model based on coupled mass and enthalpy budgets. Benn and others (Reference Benn, Fowler, Hewitt and Sevestre2019) 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, Reference Benn, Fowler, Hewitt and Sevestre2019).
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, Reference Lappegard, Kohler, Jackson and Hagen2006; Ryser and others, Reference Ryser2014; Lefeuvre and others, Reference Lefeuvre, Jackson, Lappegard and Hagen2015). 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, Reference Fatland and Lingle2002; Johnson and others, Reference Johnson, Larsen, Murphy, Arendt and Zirnheld2013; Haga and others, Reference Haga2020). 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, Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011; Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2021). An inverse method involves tuning the input parameters of ain ice-flow model to obtain the best fit to the observations. MacAyeal (Reference MacAyeal1992, Reference MacAyeal1993) 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 (Reference Ranganathan, Minchew, Meyer and Gudmundsson2020) 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, Reference Joughin, Fahnestock, MacAyeal, Bamber and Gogineni2001; Sergienko and others, Reference Sergienko, Bindschadler, Vornberger and MacAyeal2008; Gillet-Chaulet and others, Reference Gillet-Chaulet2016). 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, Reference Joughin, Fahnestock, MacAyeal, Bamber and Gogineni2001; Sergienko and others, Reference Sergienko, Bindschadler, Vornberger and MacAyeal2008). 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, Reference Shapero, Joughin, Poinar, Morlighem and Gillet-Chaulet2015). One of the few studies to invert for basal conditions of a surge-type glacier is the study by Jay-Allemand and others (Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011) 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, Reference Sevestre and Benn2015). In contrast to the well-studied surge-type glaciers found in Alaska and the Himalayas (e.g. Jay-Allemand and others, Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011; Paul, Reference Paul2020), 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 (Reference Solgaard2020) 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 (Reference Shapero, Joughin, Poinar, Morlighem and Gillet-Chaulet2015), 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 (Reference Weertman1957), 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, Reference Korsgaard2016a), 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, Reference Girod, Nuth, Kääb, McNabb and Galland2017; Hugonnet and others, Reference Hugonnet2021), and bias corrected using the ArcticDEM mosaic (Porter and others, Reference Porter2018).
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, Reference Sørensen2018). 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, Reference Sørensen2018).
Ice velocities are based on synthetic aperture radar offset tracking products from various satellite campaigns (see Table 1 and Solgaard and others (Reference Solgaard2020) 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−1. 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 (Reference Nagler, Rott, Hetzenecker, Wuite and Potin2015). b: Data from the PROMICE project (Solgaard and Kusk, Reference Solgaard and Kusk2021).
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−1. 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−1 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 (Reference Weertman1957) proposed that the relation between the basal shear stress, τb, and basal velocity, u b, can be represented as a power law
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, Reference Gillet-Chaulet2016). Maier and others (Reference Maier, Gimbert, Gillet-Chaulet and Gilbert2020) 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, Reference MacAyeal1989). 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, Reference Joughin, Fahnestock, MacAyeal, Bamber and Gogineni2001, Jouvet, Reference Jouvet2014). 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, Reference Tsai, Stewart and Thompson2015) with an added parameterisation of lateral drag (Van Der Veen and others, Reference Van Der Veen, Plummer and Stearns2011) is
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, Reference Bueler and Brown2009). The following common values are assumed: A = 9.3 × 10−25 s−1 Pa−3, n = 3, m = 3, ρ = 900 kg m−3 and $A_\ast = 2A$
(Borstad and others, Reference Borstad, Rignot, Mouginot and Schodlok2013; Bondzio and others, Reference Bondzio2016). The value for the creep parameter, A, corresponds to a temperature of − 5°C (Cuffey and Paterson, Reference Cuffey and Paterson2010). 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, Reference Langtangen2003, 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 (Reference Balise and Raymond1985) 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
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 0 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., Reference Aster, Borchers and Thurber2013, ch. 4). First-order Tikhonov regularisation is used, minimising short wavelength variations of C, similar to Jay-Allemand and others (Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011) and Gillet-Chaulet and others (Reference Gillet-Chaulet2012). The regularisation term is
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
where the regularisation parameter α ≥ 0 controls the relative weight of J 0 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., Reference Aster, Borchers and Thurber2013, ch. 10)
where c and c0 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−5, 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, Reference Hansen1992). For each dataset, the inversion is run for a range of α-values, and log J 0 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 0 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−3 Pa−2 m2/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 0 (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−5 to 1.3 × 10−4, and are weakly negatively correlated with measurement errors, reflecting that higher uncertainties lead to lower values of J 0 which are compensated by reducing α (Fig. 10). The α's sit at J 0 between 2 × 103 and 2 × 105, and at J reg between 6 × 1012 and 2 × 1013.
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 s1/3 to 5–6 mPa m−1/3 s1/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 s1/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 s1/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 s1/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 s1/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 s1/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 s1/3 ~25 km from the grounding line, decrease to 5.5 mPa m−1/3 s1/3 ~20 km from the grounding line, and increase to 8 mPa m−1/3 s1/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−1 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−1 at 40 km, to 400–1000 m a−1 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−1. 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 0–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, Reference Habermann, Truffer and Maxwell2013). 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 (Reference Gillet-Chaulet2016) 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−25 s−1 Pa−3, 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−25 s−1 Pa−3 (corresponding to − 20°C), A = 9.3 × 10−25 s−1 Pa−3 (corresponding to − 5°C) and A = 2.4 × 10−24 s−1 Pa−3 (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 (Reference Solgaard2020) 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, Reference Mankoff2020)). 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 (Reference Vallot2017) 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, Reference Lingle and Fatland2003). 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 (Reference Habermann, Truffer and Maxwell2013) 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 (Reference Mankoff2020) 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 s1/3 for the quiescent period and 1.5 to 2.5 mPa m−1/3 s1/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.
Acknowledgements
Ice velocity maps were produced as part of the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) using Copernicus Sentinel-1 SAR images distributed by ESA, and were provided by the Geological Survey of Denmark and Greenland (GEUS) at http://www.promice.dk. PROMICE is funded by the Geological Survey of Denmark and Greenland (GEUS) and the Danish Ministry of Climate, Energy and Utilities under the Danish Cooperation for Environment in the Arctic (DANCEA), and is conducted in collaboration with DTU Space (Technical University of Denmark) and Asiaq, Greenland. Ice velocity maps before 2015 are provided by the ESA Greenland Ice sheet CCI project (www.esa-icesheets-greenland-cci.org). Model output and information on how to obtain the data shown in Figures 1b and 1c are available at GEUS Dataverse (doi: 10.22008/FK2/SVMNIK). The authors thank James Lea, University of Liverpool for helpful discussions on the availability of DEMs.
Appendix
Transect velocities
A sample of transect velocities is shown in Figure 7. The transects are shown in Figure 1a as the green lines. From this, we argue that it is appropriate to use flight line velocities in the SSA flow line model.
Constant topography
In Figure 8 we show the main results for a constant topography, using the 2007 PROMICE topography. The results are slightly different from results using variable topography, but our conclusion that driving stress is not the controlling factor of the surge remains.
Increased regularisation
In Figure 9 we show the main results for an increased α for all datasets. As the main features of our conclusions are still present in this inversion, we argue that while the optimal α is determined with varying degrees of success using the L-curve, our conclusions are robust to increases in α.
Regularisation parameter
In Figure 10 we show the chosen values of α plotted against the corresponding RMS of the reported observation uncertainties. We find a weakly negative correlation as expected.
Freshwater discharge
In Figure 11 we show results from Mankoff and others (Reference Mankoff2020) of freshwater discharge for Hagen Bræ. We find no immediate relation between the decreased velocities for 2016/17 and 2017/18 and the variation in freshwater discharge.
Open access
