Validating ensemble historical simulations of Upernavik Isstrøm (1985–2019) using observations of surface velocity and elevation

Abstract

The future of tidewater glaciers in response to climate warming is one of the largest sources of uncertainty in the contribution of the Greenland ice sheet to global sea-level rise. In this study, we investigate the ability of an ice-sheet model to reproduce the past evolution of the velocity and surface elevation of a tidewater glacier, Upernavik Isstrøm, by prescribing front positions. To achieve this, we run two ensembles of simulations with a Weertman and a regularised-Coulomb friction law. We show that the ice-flow model has to include a reduction in friction in the first 15 km upstream of the ice front in fast-flowing regions to capture the trends observed during the 1985–2019 period. Without this process, the ensemble model overestimates the ice flow before the retreat of the front in 2005 and does not fully reproduce its acceleration during the retreat. This results in an overestimation of the total mass loss between 1985 and 2019 of 50% (300 vs 200 Gt). Using a variance-based sensitivity analysis, we show that uncertainties in the friction law and the ice-flow law have a greater impact on the model results than surface mass balance and initial surface elevation.

1. Introduction

Due to recent climate change, sea level is rising at an accelerating rate from 1.9 mm a⁻¹ between 1971 and 2006 to 3.7 mm a⁻¹ between 2006 and 2018 (IPCC, 2021). Compared to other sources, the dynamic contributions of the Antarctic ice sheet, and to a lesser extent, Greenland, remain significant sources of uncertainty in climate projections (IPCC, 2019). This is mainly not only due to our limited understanding of the interaction between ice dynamics and climate or ocean forcings but also due to the uncertainties inherent in ice-sheet models (Goelzer and others, 2018, 2020). First, external forcings, such as increased ocean temperatures (Wood and others, 2021) or surface melting (Hofer and others, 2020), can cause changes in boundary conditions, resulting in significant alterations in glacier dynamics. In addition, the future of these forcings and their interactions with the dynamics are complex and highly non-linear (Flowers, 2018; Davison and others, 2019). Second, the initial state of the model, parameterisations and input parameters remain poorly constrained, resulting in uncertainties that propagate into future projections of Greenland mass loss.

In geophysical models, parameterisations are employed to establish a connection between small-scale and large-scale physics. This allows for the reduction of computing time while still accurately representing small-scale processes. In ice-sheet models, two important parameterisations are the flow law, which describes ice rheology and friction law, which describes the interaction between the ice and its bed. Uncertainties attached to these parameterisations depend on the input parameters (such as the ice rate factor or friction coefficient) as well as their form and how they relate to other prognostic variables (e.g. isotropic vs anisotropic flow laws, Weertman vs effective pressure-dependent friction laws). Many ice-sheet models now use time-independent inverse methods to constrain these input parameters from surface observations and initialise the model to present day (Gillet-Chaulet and others, 2012; Larour and others, 2012; Cornford and others, 2013; Pattyn, 2017; Quiquet and others, 2018). By construction, it is generally not possible to constrain the form of these parameterisations with a single inversion (Joughin and others, 2004; Gillet-Chaulet and others, 2016), and assessing the performance of the method in recovering the best parameters is challenging. This is due to the fact that the inversion is often restricted to a single parameter, the others being fixed (Babaniyi and others, 2021). Expanding the number of controlled parameters could exacerbate the ill-posed nature of the problem and potentially result in mixtures between these parameters (Ranganathan and others, 2021). Additionally, to increase spatial coverage, observations are often compiled from several dates, which can lead to inconsistencies between the datasets (Seroussi and others, 2011) and decreases the number of independent datasets available to validate the models. All of these factors raise questions about the robustness of the forecasts produced by ice-sheet models for climate projections.

The recent increase in satellite observations, mainly of surface elevation and velocities, provides new opportunities to investigate the performance of ice-sheet models on decadal timescales. Several studies have undertaken comparisons between historical model outputs and observational data, examining specific cases such as the effects of alterations on the ice front position (Bondzio and others, 2017; Haubner and others, 2018), grounding line retreat (Joughin and others, 2019; Åkesson and others, 2021) and variations in surface mass balance (SMB) (Peyaud and others, 2020). Nevertheless, establishing suitable quantitative measurements for categorising and distinguishing model results poses a challenge. Aschwanden and others (2013) demonstrated that validating a model is more effective when utilising spatially dense data, such as velocity, altitude and altitude change, as opposed to global variables. Otherwise, there is a risk of mischaracterising the thermal and dynamical state of the glacier.

Properly quantifying the reliability of a model or discriminating between different model formulations or parameterisations requires a proper assessment of the modelling uncertainty. Uncertainty can arise from various sources, including random noise in the observational data, sparsity in space or time and the use of a simplified or incomplete mathematical model of system dynamics. To quantify and propagate uncertainties, ensemble modelling is increasingly used in geophysical models (Brown and others, 2010). This can be done with a single model by perturbing the model initial state, parameters, boundary conditions or forcings (Schlegel and others, 2018; Hill and others, 2021). The advantage resides in the ease of conducting numerous simulations and assembling a substantial ensemble, a contrast to multi-model ensembles. However, the primary constraint is the computational time required. Nevertheless, these ensembles, in general, are unable to account the uncertainties associated with the model structure (Weisheimer and others, 2011). For this reason, multi-model ensembles allow for a wider range of uncertainties to be explored, such as model formulation and numerical implementation, and are widely used in climate projections, including the most recent forecasts of ice-sheet evolution over the 21st century (Goelzer and others, 2020; Seroussi and others, 2020). The results for the Greenland ice sheet (GrIS) showed that it would contribute 90 ± 50 mm of sea-level rise by 2100 (Goelzer and others, 2020). Of this 100 mm of sea-level rise uncertainty, the model uncertainty explains a similar portion of the ensemble gap (40 mm of sea-level rise by 2100) as the atmospheric forcing uncertainty and twice as much as the ensemble gap due to the oceanic forcing uncertainty (19 mm). More significant differences between ice-flow models appear near tidewater glaciers. At these locations, the amount of ice available for calving varies between initial model states, as does the transmission of dynamic perturbations on the margins to inland, which is strongly dependent on the non-linearity of the friction law (Price and others, 2008).

The objective of this study is to test the ability of the Elmer/Ice model, in a configuration identical to the one used at large spatial scale, to reproduce past trends in available satellite observations. For this purpose, we study the evolution of the Upernavik Isstrøm (UI) catchment during the period 1985–2019. UI is a tidewater glacier in the northwestern sector of Greenland that is now divided into five different branches which are named here as in Mouginot and others (2019), from north to south: UI-NN, UI-N, UI-C, UI-S and UI-SS (Fig. 1). Due to the very different dynamics of the three main branches (UI-N, UI-C and UI-S), UI allows us to conduct several tidewater glacier studies in one. Furthermore, this glacier has experienced significant mass loss since 1985, contributing to a sea-level rise of 0.61 mm, which accounts for ~4% of Greenland's total contribution during this period, highlighting notable temporal changes. The selection of this region, marked by significant spatio-temporal variations in mass and heterogeneous, dynamic behaviours, serves to mitigate overconfidence in the model results. Finally, the large amount of satellite observations collected during the period 1985–2019 makes UI a good case study to evaluate the performance of a large-scale ice-sheet model in reproducing the available observations of a local glacier. After describing the model setup and the observational datasets, we investigate in particular the possibility of using observations to discriminate the performance of two friction laws in producing realistic velocity and surface elevation evolutions. For each friction law, the uncertainties related to the initial state and other model input parameters are taken into account using an ensemble approach. Subsequently, we explore the sensitivity of the simulated ice velocity, surface elevation, ice volume and ice discharge to the input parameters and the possibilities to better constrain them through a statistical analysis.

Figure 1. Left: GrIS drainage basins with the catchment of UI in red. The blue box is the validation area shown in the right with the different catchments (UI-NN, UI-N, UI-C, UI-S and UI-SS), the front positions between 1985 and 2018 (Wood and others, 2021) and the surface ice velocities (Mouginot and others, 2019) overlaid on a Landsat image (2017-08-13).

2. Data and method

We study the evolution of UI for the period 1985–2019. Historically, all UI branches formed one single outlet glacier from at latest 1849, which corresponds to the first reported observation and the maximum extent of the Little Ice Age (Weidick, 1958), until 1931 where the branches separated. In the 1930s, the three northern glaciers (UI-NN, UI-N and UI-C) separated from the southern ones (UI-S and UI-SS). As of 1976, UI-C had separated from UI-N and UI-NN, both of which underwent a split in 2010 (Andresen and others, 2014). For this last retreat, an increase in ice flow (King and others, 2018; Mankoff and others, 2019; Mouginot and others, 2019) and a rapid decrease in surface elevation was documented. The position of the UI-S front has remained stable over the period 1985–2019. In total, UI has retreated more than 35 km in the last 170 years. The choice to start the study in 1985 is based on the fact that the fronts are available from that date and that a DEM is available at that time (Korsgaard and others, 2016). The notations and constants employed in this study are detailed in appendix Tables 1 and 2, respectively.

2.1 Model description

We use the parallel finite-element code Elmer/Ice (Gagliardini and others, 2013) to model the evolution of UI from 1985 to 2019. The model domain corresponds to the UI catchment shown in Figure 1.

We use an anisotropic mesh adaptation scheme that equidistributes the interpolation error of the observed surface velocities and thickness (Frey and Alauzet, 2005). This scheme leads to decreasing the element size in the directions of the highest curvatures, and the resulting mean element size ranges from 150 to 600 m for the first 50 km from the margin and is equal to ~5 km further upstream. A time step of 1 d is used. The position of the front is forced at each time step by interpolating the annual observations of Wood and others (2021). The mesh is fixed and the effective ice–ocean boundary is determined by the edges between glaciated and deglaciated elements, and thus changes discretely over time. The deglaciated elements are then deactivated and not numerically resolved.

To address the force balance, we employ the Shelfy-stream approximation (MacAyeal, 1989), incorporating Glen's constitutive flow law (Glen and Perutz, 1955). This non-linear relationship depends on the Glen exponent n, fixed at a value of 3, and the rate factor A, which is a function of ice temperature. Due to the complexity to initialise a temperature field consistent with other model variables (e.g. Schäfer and others, 2012; Zhao and others, 2018) and the fact it depends on processes that are very poorly represented in models (e.g. cryo-hydrologic warming; Phillips and others, 2010), the thermo-mechanical coupling is neglected and we assume that A is constant in time. Furthermore, Bondzio and others (2017) found that the gradual warming of the shear margin contributed only 5–10% to the acceleration of Jakobshavn Isbræ between 1985 and 2016. As for GrIS simulations with the model Elmer/Ice (Goelzer and others, 2018), A is initialised using a present-day 3-D ice temperature field computed with SICOPOLIS (Greve, 1997) after a palaeo-climatic spin-up and using the values given by Cuffey and Paterson (2010) for the prefactors and activation energies. Uncertainties associated with this flow law are usually accounted for by the enhancement factor E, a scale factor to A. For further details, see the Supplementary materials.

2.1.1 Friction laws

In this study, two different friction laws, relating the velocity u and the basal shear stress τ_b, are used for grounded areas:

• • a Weertman friction law (Weertman, 1957):

  (1)$${\boldsymbol \tau_{\rm b}} = -\beta_{\rm W} \Vert {\bf u}\Vert ^{{1}/{m}} {{\bf u}\over \Vert {\bf u}\Vert }$$

• • a regularised-Coulomb friction law (Joughin and others, 2019):

  (2)$${\boldsymbol \tau_{\rm b}} = - \beta_{{\rm RC}} \left({\Vert {\bf u}\Vert \over \Vert {\bf u}\Vert + u_0}\right)^{{1}/{m}} {{\bf u}\over \Vert {\bf u}\Vert }$$

Both Eqns (1) and (2) depend on a friction coefficient, β _W or β _RC, m is a positive exponent and u ₀ is a threshold velocity (Joughin and others, 2019). As the datasets are too incomplete in 1985, the friction coefficients β _W and β _RC are initialised with the adjoint inverse method using BedMachine v3 (Morlighem and others, 2017) for the topography and surface velocity datasets that cover the period 2010–2018.

Equation (2) exhibits two asymptotic behaviours depending on the threshold velocity u ₀. When ‖u‖ ≪ u ₀, the relation tends towards a Weertman regime with $\beta _{\rm W} = \beta _{{\rm RC}} u_0^{-1/m}$. Conversely, when ‖u‖ ≫ u ₀, it approaches a Coulomb regime with τ _b = −β _RC (u/‖u‖). In contrast to other similar friction laws that explicitly depend on the effective pressure N (Schoof, 2005; Gagliardini and others, 2007) and therefore require a model for the basal water pressure, here, its effect is accounted for in the parameters β _RC and u ₀.

However, in Joughin and others (2019), β _RC is made proportional to the height above flotation, h _af, only when h _af drops below a specified threshold. This limitation confines the effect to the region immediately upstream of the grounding line. This parameterisation is equivalent to assuming a perfect hydrological connection between the subglacial drainage system and the ocean in these areas. It results in a smooth decrease of the basal shear stress τ _b to 0 as the ice column approaches flotation.

2.1.2 Friction parameterisation

Here, we found that this parameterisation was not suitable, as the β _RC coefficients were initialised based on recent observations. This initialisation would lead to a rise in the coefficients associated with the initial state back in 1985, a period characterised by greater ice thicknesses and high uncertainties. As a consequence, the outcomes would become notably susceptible to inversions in regions where the current ice column closely approaches flotation. However, to maintain a reduction of the basal shear stress as ice approaches flotation, β _RC is made a function of the distance to the ice front, d, as

(3)$$\beta_{{\rm RC}} = \beta_{{\rm ref}} + \beta_{{\rm lim}} {d\over d + d_{{\rm lim}}}$$

where β _ref is a time-independent reference field obtained from the inversion, and β _lim and d _lim are two parameters that are constrained from the results of the inversion under the central flowlines of the three main outlet glaciers. This parameterisation constrains the changes to fast-flowing areas, with magnitudes consistent with the findings of Habermann and others (2013) regarding the evolving basal conditions during the speed-up of Jakobshavn Isbræ. Further details are provided in the Supplementary materials.

2.1.3 External forcing

For the evolution of the bottom and top free surfaces, we solve the continuity equation for the ice thickness (Supplementary materials) using the flotation condition. As we do not resolve the thermo-mechanical coupling, we neglect the basal melt rate in grounded areas. We also set it to 0 in floating areas, as the floating areas remain small during our simulations. For the SMB $\dot {a}_{\rm s}$, we prescribe the annual means from two different regional climate models, RACMO (Noël and others, 2018) and MAR (Fettweis and others, 2017). The influence of ocean forcing is implicitly taken into account by the prescribed position of the front. This position is indeed determined by processes such as melting and calving that occur in contact with the ocean, in addition to the velocity of the glacier at the front.

2.2 Validation data

To validate our model, we compile surface velocities and elevations obtained from various airborne and spaceborne sensors between 1985 and 2019 in our validation area (see Fig. 1). Both datasets are averaged annually to produce time series of the annual mean posted on a common regular grid with a horizontal resolution of 150 m. This arithmetic mean may differ from the mean annual velocity or elevation if the distribution is not uniform, i.e. if there are more observations in winter, the arithmetic mean will be more representative of the winter mean. The velocities are produced from speckle or feature tracking using 16 different satellite sensors (Joughin and others, 2010; Mouginot and others, 2019; Derkacheva and others, 2020; Howat, 2020). For the surface elevation, we compiled existing airborne and spaceborne laser altimetry data, interferometric single-pass DEMs and photogrammetric DEMs based on aerial photographs (Studinger, 2014; Zwally and others, 2014; Fenty and others, 2016; Korsgaard and others, 2016; Howat and others, 2017; Blair and Hofton, 2019; OMG, 2020). We also produced a time series of photogrammetric DEMs based on ASTER observations. For each dataset, time-series plots, details of spurious value filtering and additional information on the ASTER observations are provided in Supplementary materials.

Finally, we compiled published ice discharges from King and others (2018), Mankoff and others (2019) and Mouginot and others (2019). The ice discharges are computed as the ice flux crossing flux gates near the outlets, usually assuming uniform velocity profiles along the ice thickness. The results depend on the velocity and thickness datasets, accounting for differences between them. Notably, Mouginot and others (2019) employ thicknesses directly obtained from radar measurements, whereas the other two studies rely on BedMachine data. The annual ice mass loss is then derived from these ice discharges using the input–output method at the catchment scale, as described by Mouginot and others (2019). This is achieved by subtracting the previously calculated discharge from the SMB modelled by the regional atmospheric model RACMO. The total observed ice mass loss is the integration of this annual ice mass loss observed since 1985.

3. Ensemble modelling

In this section, we outline the method used to assess the model's ability to replicate observations of surface velocities and elevations. To achieve this, we conduct two ensembles of numerical simulations for the period between 1985 and 2019, employing the two different friction laws described above (Eqns (1) and (2)). These ensembles are referred to as the Weertman ensemble (WE) and the regularised-Coulomb ensemble (RCE), respectively. In the Supplementary materials, we provide further details on some of the issues involved in this method.

To create both ensembles, we performed small perturbations of model parameters, initial conditions and external forcings. We examine the response of the ensembles to the retreat of the glacier front in terms of ice-flow velocity, surface elevation and ice flux to the ocean, and compare them to the validation datasets. To that purpose, we defined metrics and methods to evaluate the two ensembles and their members, which are described at the end of Section 3.

3.1 Model parameters

Uncertainties in the flow law may originate from its form (isotropic vs anisotropic; Lliboutry and Duval, 1985, from the Glen exponent; Glen and Perutz, 1955; Gillet-Chaulet and others, 2011; Millstein and others, 2022), or the initial temperature field and the Arrhenius parameters used to compute A. Because of the assumption that A depends exclusively on temperature, its spatial variability does not go beyond its dependence on temperature. Additionally, we refrained from varying the Glen exponent n, as modifying n would also impact A, given that A's units are intricately tied to the value of n. Consequently, distinguishing between the effects of altering n and those of altering A would become challenging. However, the flow law uncertainty is explored through the enhancement factor E using a continuous uniform distribution (see Fig. 2) between 0.5 and 5 (based on Huybrechts, 1990, Ma and others, 2010, Le clec'h and others, 2019). We therefore choose not to consider the temporal and spatial variability of the flow law parameters.

Figure 2. Probability distributions for uncertain parameters included in our analysis. Red box: parameters related to the ice dynamics; orange box: parameters related to initial surface field; green box: parameters related to historical forcing.

For the two friction laws, we investigate the dependence on the initial friction coefficient fields, and on the parameters m (with a discrete uniform distribution between 1 and 5) and u ₀ ($\log _{10}( u_0) \sim {{\cal N}}( \log _{10}( 300) ,\; \, 0.25)$). This parameter range for m and u ₀ is similar to that used by Joughin and others (2019) and u ₀ is assumed to be uniform in space.

3.1.1 Ensemble of friction fields

For the friction coefficient field β _W, we use an ensemble generated at the GrIS scale using classical inverse methods (Gillet-Chaulet and others, 2012). The goal of these inverse methods is to minimise the difference between observed and modelled surface velocities. These inversions are conducted using the Weertman friction law (Eqn (1)) with different pairs of values for E and m. The results are sensitive to the observed data and their uncertainty, as well as the regularisation weights λ _reg and λ _div, obtained through L-curve methods (see Gillet-Chaulet and others, 2012, for more details). As λ _reg increases, the solution becomes smoother, while increasing λ _div minimises the difference between the observed and modelled ice flux divergence. From an L-curve analysis, we choose to use the following continuous distributions for the weights λ _reg and λ _div: $\log _{10}( \lambda _{{\rm reg}}) \sim {\cal N}( 5,\; \, 1),\ \log _{10}( \lambda _{{\rm div}}) $ $\sim {\cal N}( -5,\; \, 0.5)$ (Fig. 2). To account for the uncertainty in the input surface velocity observations u _obs, we use five different datasets: a mosaic compiled from observations between 1995 and 2015 (Joughin and others, 2016), and four yearly datasets from 2015 to 2018 (Mouginot and others, 2017). We then randomly select from these five fields (Fig. 2).

β _RC is derived from the same inversions, assuming that the basal shear stress τ _b and basal velocity u from the inversion are identical in both friction laws, leading to:

(4)$$\beta_{{\rm RC}} = \beta_{\rm W} ( \Vert {\bf u}\Vert + u_0) ^{{1}/{m}}$$

3.1.2 Initialisation of the friction coefficients in front retreat areas

These inversions require the geometry to be prescribed. As the surface altitude observations were incomplete before the 2000s (see Supplementary material, Fig. S4), and in order to be consistent with the velocity datasets, we used the topography given by BedMachine v3. The position of the fronts in this dataset roughly corresponds to our observed front positions in the 2010s. Therefore, it is necessary to initialise the friction coefficients in the areas that were not covered by ice during the inversion, particularly where the front has retreated since 2005. This is a common challenge for studies using inversions to initialise the model with data more recent than the initial state (e.g. Haubner and others, 2018). Due to the non-linearities of the friction laws, there is no single solution that would lead to the same initial state in 1985 for the Weertman and regularised Coulomb laws.

For WE, β _W is extrapolated by using a linear fit between the results of the inversion and the bedrock elevation. Indeed, the results of the inversions show a correlation with the bedrock elevation, friction being generally lower at the bottom of the valleys in the deepest areas. In 1985, the front of UI-N and UI-C was located in a place deeper than it actually is, so the extrapolation leads to a friction that is almost zero; however, for UI-S the front was located in a bedrock high, leading to a small but non-negligible friction in the area where the front has retreated.

For RCE, we assume β _RC,ref = 0 in the front retreat areas as no additional information is available. In this area, the friction coefficient depends solely on the distance to the front according to Eqn (3). In 1985, this leads to slightly lower friction for UI-S and slightly larger friction for UI-N and UI-C compared to the extrapolation for WE. However, importantly the friction coefficient for RCE in this area evolves with time according to Eqn (3).

3.2 Initial topography

For the bed elevation, we use BedMachine v3. In the validation area, the bed is reconstructed using the mass conservation method and reported uncertainties are up to 200 m. However, incorporating the uncertainty of the bed would introduce an additional layer of complexity. Generally, generating representative ensembles of beds necessitates advanced geostatistical techniques (e.g. MacKie and others, 2021). For simplicity, we neglect this uncertainty.

As the initial dataset for 1985 elevations is incomplete, missing values have been filled with more recent observations. To address these initial inconsistencies, the model is run for a short time period under constant SMB forcing, termed ‘relaxation’. We vary both the SMB $\dot {a}_{\rm s}$ and the relaxation time t _relax. For $\dot {a}_{\rm s}$, we use ten different fields from two regional climate models (SMB _model: RACMO and MAR) averaged over different time periods (SMB _date: 10 year averages: 1959–68, 1969–78, 1979–88, 1989–98, 1999–2008); for the relaxation time, it varies from 5 to 10 years (see Fig. 2), which is an empirical compromise between the dissipation of anomalies and the strong influence of the SMB.

After random selection of t _relax, SMB _model and SMB _date, a relaxation is performed for WE utilising the Weertman friction law. Subsequently, for RCE, a new relaxation is initiated by restarting from the initial topography state of WE and extended for 2 years with the regularised-Coulomb friction law. This methodology is applied to harmonise and minimise disparities in the initial states between the two ensembles.

3.3 Historical forcings

For $\dot {a}_{\rm s}$ between 1985 and 2019, we use the annual values given by the two regional climate models MAR and RACMO. Both products have already been downscaled to 1 km resolution. Both models have been forced by global reanalyses and differences between the values integrated over the model domain and period are in the order of 10%. In addition, the SMB from RACMO showing greater spatial variability than MAR as a function of altitude. We employ the same SMB model for both the historical simulation from 1985 to 2019 and the relaxation. There is no surface-elevation feedback.

3.4 Ensemble evaluation

For each ensemble, the parameters are sampled from the probability distributions discussed in the previous section (Fig. 2), using a Latin hypercube sampling (McKay and others, 1979).

3.4.1 Performance metrics

For each ensemble using a different friction law, we establish 120 members. Each member is initialised in 1985 using the methodology described above and then propagated forward in time to 2019. At the end of each year, we employ a bi-linear interpolation method to map our current snapshot model states onto the observation grid, utilising a netCDF data structure. Then, the results from each member and the ensemble mean are evaluated against satellite observations and scored using the following metrics, with Q a given physical quantity: the bias (bias_Q), the RMSE_Q and the mean absolute error (MAE_Q); bias_Q is the averaged difference between a predicted quantity Q _m and its observed equivalent Q _o; RMSE_Q and MAE_Q express the accuracy, i.e. the ability of the modelled quantity to match an observation, but the MAE_Q is less sensitive to large mismatch:

(5)$$\rm {bias}_Q = {1\over k}\sum_{i = 1}^{k}{\left(Q_{\rm m}^i - Q_{\rm o}^i \right)}$$

(6)$$\rm {RMSE}_Q = \sqrt{{1\over k}\sum_{i = 1}^{n}{\left(Q_{\rm m}^i - Q_{\rm o}^i \right)^2}}$$

(7)$$\rm {MAE}_Q = {1\over k}\sum_{i = 1}^{n}{\left\vert Q_{\rm m}^i - Q_{\rm o}^i \right\vert }$$

where k can be either a number of observations in time or space, or a combination of both.

To evaluate the performance of the whole ensembles, we use the continuous rank probability score (CRPS_Q), which highlights the accuracy and the sharpness (opposite of uncertainty/spread) of an ensemble, meaning that lower values are obtained when the ensemble mean is close to the observations and the ensemble spread similar to the observation uncertainty (Brown, 1974; Matheson and Winkler, 1976; Unger, 1985; Bouttier, 1994; Hersbach, 2000):

(8)$$\rm {CRPS}_Q = {1\over k}\sum_{i = 1}^{k}{\int_{{\opf R}} \left(F_{\rm m}^i( Q) - F_{\rm o}^i( Q) \right)^2\, {\rm d}Q}$$

where for a quantity Q at point i, $F_{\rm m}^i( Q)$ is the cumulative distribution function of the ensemble model and $F_{\rm o}^i( Q)$ is the cumulative distribution function of the observation. At a given time and location, we have a single observation and it is common to neglect the uncertainty associated with the observation to compute the CRPS (Brown, 1974; Matheson and Winkler, 1976; Unger, 1985; Bouttier, 1994; Hersbach, 2000) therefore $F_{\rm o}^i( Q)$ is a Heaviside step function: $F_{\rm o}^i( Q) = 1$ if Q > Q _obs and $F_{\rm o}^i( Q) = 0$ if Q ≤ Q _obs.

Our metrics do not take into account the observation uncertainties. We made this choice to avoid giving even more weight to recent data, which are already more numerous. This may give too much weight to the old data in relation to their accuracy, but it allows us to examine the model's ability to reproduce the whole trend, rather than looking at its ability to reproduce current data very accurately.

Ensemble modelling is computationally resource-intensive, making it crucial to assess the added value of an ensemble compared to a single member. This evaluation can be conducted by comparing MAE_Q and CRPS_Q. For a one-member ensemble, such as a deterministic simulation, CRPS_Q equals MAE_Q. Consequently, if the CRPS_Q closely resembles the MAE_Q of the ensemble mean, it indicates that the dispersion of the ensemble does not offer additional information. This similarity could result from biased observations or an overly confident ensemble. In the latter case, the confidence interval is too narrow, suggesting the underestimation of uncertainty in one or more model parameters, or that our model assumptions are too stringent.

3.4.2 Sensitivity analysis

Finally we also conduct a sensitivity analysis to identify the parameters that exert the most influence on the modelled quantities. In this analysis, we use the first-order Sobol indices at a given time (S _i), expressing the sensitivity of an output Y to the input parameter X _i:

(9)$$S_i = {{{\rm Var}}\, {( {{\rm E}}\, [ Y\vert X_i] ) }\over {{\rm Var}}\, {Y}}$$

where Var Y is the variance of Y computed across the entire ensemble. To compute Var (E [Y|X _i]), we first create subgroups of members with a similar input parameter X _i, then average each subgroup and finally calculate the variance of these different subgroups means. The Sobol indices are normalised so that their sum is equal to 1, neglecting the influence of Sobol indices of higher order. A high Sobol index indicates that the input parameter explain a large fraction of the ensemble spread.

4. Results

4.1 Impact of the friction laws

4.1.1 Global variables

Figure 3 shows the total ice discharge and mass change from 1985 to 2019. WE is therefore unable to reproduce the initial state or the increase in ice discharge between 2005 and 2010 as reconstructed by King and others (2018), Mankoff and others (2019) or Mouginot and others (2019). The initial WE flux is ~3 Gt a⁻¹ (or ~30%) too high compared to the observations. Moreover, the discharge increases by only 2 Gt a⁻¹ (or 12.5%) after the large front retreat starting in 2005.

Figure 3. Ice discharge (top graph) and cumulative ice mass change (bottom graph) for RCE (red) and WE (blue) between 1986 and 2019, with mean in solid line and the shading include 95% of the ensemble members, against different observation: Mouginot (+), Mankoff (Y) and King (×). On the right, histograms of ice Discharge and ice mass change in 2019.

We find a better agreement between RCE and the discharge reconstructions during the initial period : the modelled discharges are close to Mankoff and others (2019) and Mouginot and others (2019) while King and others (2018) give slightly larger values. After 2005, the discharge of the ensemble mean matches (King and others, 2018; Mankoff and others, 2019), but is lower (Mouginot and others, 2019) which differs greatly from the two others after 2010. The reason for the difference is that King and others (2018); Mankoff and others (2019) used BedMachine thicknesses as input, leading to a better match, while Mouginot and others (2019) used radar flight lines with a deeper bed.

These differences are reflected in the mass change, for which RCE is in good agreement with the different datasets, showing that UI has lost ~220 Gt between 1985 and 2019. WE significantly overestimates mass loss due to excessively high ice flow before 2005.

4.1.2 Initial period dynamics (1985–2005)

We now discuss the results obtained for the surface velocity and surface elevation during the initial period from 1985 to 2005 where observations show relatively stable conditions as no large front retreat has been observed for the three major branches, UI-N, UI-C and UI-S. The two ensembles WE and RCE are evaluated by scoring the ensemble against the velocity and surface elevation time series.

The bias_u, the MAE_u and the CRPS_u obtained for the velocity and computed using all the observations for this initial time period are shown in Figure 4. In the slow flowing areas, WE agrees well with the observations with absolute metric values between 0 m and 50 m a⁻¹. In the fast sections of the streams, a positive bias exceeding 1000 m a⁻¹ is evident in the final 20 km of UI-N, where observed velocities reach 3000 m a⁻¹. This influence notably impacts both the MAE_u and the CRPS_u. UI-C has also a smaller positive bias of ~300 m a⁻¹, respectively a relative error of 10% while UI-S and UI-SS have a negative bias of ~−300 m a⁻¹ (i.e. relative errors of 10 and 50% respectively).

Figure 4. Surface velocity bias (top), MAE (middle) of the ensemble mean and CRPS (bottom) for WE (left) and RCE (right) during the period 1985–2005. Points A (+) and B (×) are used in (Figs 6, 7) as representative of UI-N and UI-S respectively. The grey and black lines in the first row are the 200 and 1000 m a⁻¹ velocity contours computed from RCE 1985–2005 average.

RCE appears to perform better than WE in modelling flow velocities in fast flowing areas (Fig. 4, right column), while performing equally well for slow flowing areas. This observation aligns with the understanding that Coulomb's law converges towards Weertman's law for these lower velocities (‖u‖ ≪ u ₀). Nevertheless, RCE still presents negative biases of −250 m a⁻¹ on UI-S and UI-SS, and a positive bias of 250 m a⁻¹ on UI-N, while no bias is visible for UI-C.

Both ensembles have larger MAE_u than CRPS_u in the streams, meaning that the forecast uncertainty included in the CRPS provides more information than the ensemble mean alone. However, the CRPS_u in the shear margins, at the transition between slow and fast flow, remain very high, especially near the front, with values ~700 m a⁻¹ for UI-N while they are lower than 300 m a⁻¹ in the interior of the stream. The high bias observed in the shear margins could be due to precisely matching large velocity gradients with the model, possibly due to neglected specific physical processes in those areas; or to uncertainties in the observations as it is difficult to accurately capture steep transitions in surface flow velocity using image cross-correlation on satellite imagery (Millan and others, 2019). High-resolution imagery (e.g. CNES's Pléiades or SPOT) would be needed to better map this challenging region. Outside the streams, the differences between MAE_u and CRPS_u are very small because of the low velocities and the small differences between the different members.

For the surface elevation (Fig. 5), WE shows a negative bias between −25 and −50 m on the upper part of UI-N and UI-C, i.e. the model underestimates the surface elevation. This bias is more pronounced inside the fast flowing regions and in the upper catchment of UI-N. This underestimation is associated with the positive bias on the velocity (Fig. 4). Similarly, a clear positive bias (~50 m) is visible on the lower part of UI-S in relation to a negative bias for the velocity.

Figure 5. Same as Figure 4 for surface elevation.

RCE has a lower overall bias than WE and does not exhibit a positive bias on UI-S. However, we still find an underestimation of the elevation (−75 m) close to the UI-N front.

For both ensembles, kilometre-scale undulations are visible for all scores associated with the surface elevation. As there is no such pattern in the observations, we associate it with uncertainties in the bed elevations that have not been taken into account, leading to an overconfidence of the ensembles with similar biases for all members. Indeed, additional examinations conducted on the RCE by introducing minor perturbations to the bed configuration revealed that these undulations ceased for the CRPS, although they persisted with the MAE. Otherwise, there is no major difference between the MAE_zs of the ensemble mean and the CRPS_zs for both ensembles, meaning that the ensemble does not provide any particular added value when evaluating elevations, unlike velocities.

4.1.3 Front retreat dynamics (1985–2019)

We are now interested in the temporal evolution and in particular in the ability of the two ensembles to reproduce the changes in velocity and elevation following the retreat of the front between 2005 and 2010. Figures 6 and 7 shows the evolution of the velocity, basal shear stress and surface elevation for two points A and B located 5.5 km upstream of the 2019 front position on UI-N and UI-S, respectively, as these two branches have different histories and behaviours. Due to the challenge of properly determining errors in the annual average observations, we have not taken them into account in the scores (Eqns (5)–(8)). However, for (Figs 6, 7), we calculated an estimate of the error for surface velocities and surface elevations (details are provided in the Supplementary materials).

Figure 6. Surface velocity (top), basal shear stress (middle) and surface elevation (bottom) at point A (see Figs 4, 5) of WE (in blue), RCE ensemble (in red) and observations (black dots with an estimate of the uncertainty in grey). For WE and RCE the mean is represented in solid line and the shading include 95% of the ensemble members.

Figure 7. Same as Figure 6 for point B (see Figs 4, 5).

For point A on UI-N, WE has a significant positive bias for the velocity with respect to the pre-retreat observations, i.e. for the period 1985–2005, in agreement with what is shown in Figure 4. The increase in modelled velocity for WE between 2005 and 2012 is ~1500 m a⁻¹ (from ~3500 to ~5000 m a⁻¹ for the ensemble mean), while the observations show a larger increase over the same period (from 2000 to 4200 m a⁻¹). RCE generally performs better than WE, showing mean velocities of ~2300 and 4200 m a⁻¹ before and after the retreat, respectively. This aligns more closely with observations, indicating RCE's effectiveness in capturing the velocity increase. We mainly attribute the difference in initial velocities to local differences in the basal shear stress, τ _b, due to the parameterisation of the friction coefficient with respect to the distance to the front in RCE. For WE, τ _b is very small between 0 and 40 kPa and corresponds to a large spread in velocity between 2500 and 5000 m a⁻¹. For RCE, the parameterisation Eqn (3) has led to an increase of β _RC due to the more advanced front during the initial period, as a result τ _b is slightly higher between 10 and 70 kPa. It corresponds to lower velocities in better agreement with the observations and with a lower spread. After 2007, the speed-up is a response to the loss of buttressing at the front as it retreats. For WE, τ _b increases with the velocity in agreement with the friction law Eqn (1). For RCE, the parameterisation decreases β _RC as the front retreats, leading to a decrease of τ _b, and a higher speed-up. The dispersion of RCE also increases during the retreat, and we attribute this to a higher sensitivity to τ _b as it decreases. As expected, because the friction coefficient has been calibrated to reproduce recent observations, both ensembles lead to a relatively good agreement with the observations after 2012, RCE leading to slightly lower velocities.

For the surface elevation, both ensembles start with a similar surface elevation in 1985 and underestimate by ~100 m the first available observations in 2003–04. In relation to the higher velocities, WE exhibits slightly larger thinning rates than RCE during this initial period. Both are unable to fully reproduce the amplitude of the observed thinning, which is over 160 m between 2005 and 2015 at this location (or 8% of the thickness), against 80 m for RCE and 100 m for WE (or 5% of the thickness). As a result, WE is in better agreement with the observations than RCE after 2015, while velocity changes have been smaller. In contrast to velocity, the initial elevations of WE are equally dispersed as those of RCE. We can observe that both dispersions are smaller after the front retreat. This means that members of RCE and WE with higher initial surface elevation, have a larger decrease than members with a lower initial elevation, so that the spread decreases. Expressed in terms of differences in thickness at the end of the period, the differences are relatively small. At this location, the ice is ~1100 m thick, so a difference of 50 m corresponds to an error of only 5%, as does the relative error in velocity.

Point B on UI-S shows different results (Fig. 7). The front retreats less than 1 km between 1985 and 2019 and the velocity shows a globally decreasing trend from ~1900 to 1700 m a⁻¹ with some inter-annual variability. Contrary to point A, RCE leads to velocities higher than WE initially despite a higher basal shear stress τ _b due to the initially more advanced front position. With RCE, there is a decrease of τ _b from 20 to 10 kPa due to the effect of the parameterisation, which reduces friction when the front retreats. However, there is no real trend on the mean velocity which is ~1600 m a⁻¹, and thus has a negative bias during the initial period, in agreement with (Fig. 4). WE shows a much larger velocity bias at the beginning and the mean velocity increases as the front retreats from 1200 m a⁻¹ in 1985 to 1600 m a⁻¹ in 2005. We attribute this result to the difference in the extrapolation of the friction coefficient in areas that were not glaciated during the inversion. This extrapolation leads to higher τ _b near the front which slows down the whole stream. As the front retreats, the velocity increases and return to values close to the present-day observations similar to RCE. WE shows also more inter-annual variability, certainly associated with the changing conditions near the front. For surface elevation, the initial disparity between WE and RCE surpasses 70 m, indicating a rapid adjustment during the additional 2 year relaxation for RCE. This adjustment is likely associated with significant velocity changes induced by differences near the front. RCE shows a decreasing trend between 1985 and 2019 and slightly underestimates the observed surface elevation. As depicted in Figure 5, within the initial 10 km near the front, WE displays a notable positive bias in surface elevation. However, with the slight retreat of the front, a downward trend emerges, accompanied by an increase in velocity. Consequently, the mean surface elevation becomes relatively close to the observations in 2019. Finally, in terms of the initial spread, WE exhibits a considerably larger range (~110 m) compared to RCE (40 m).

Finally, to summarise, Figure 8 highlights the difference in performance between RCE and WE in terms of velocity and elevation. This figure clearly shows that the RMSEs on velocity obtained by RCE are better than those of WE. Although some WE members seem to have reasonable scores (RMSE_u between 170 and 390 m a⁻¹ and RMSE_zs between 24 and 52 m for elevation), RCE members have better overall scores (RMSE_u between 120 and 160 m a⁻¹ and RMSE_zs between 24 and 34 m). For the CRPSs computed over the whole validation area and period, we always find better scores for RCE (CRPS_u of 43.1 m a⁻¹ and CRPS_zs of 16.7 m) than for WE (CRPS_u of 51.8 m a⁻¹ and CRPS_zs of 19.6 m).

Figure 8. Distribution of WE and RCE RMSE for the surface velocity (top) and the surface elevation (bottom) in the validation area (Fig. 1) over the full period (1985–2019). The vertical line shows the values of the ensemble mean.

4.2 Sensitivity analysis

4.2.1 Global variables

Here, we further the analysis by assessing the sensitivity of the results to input parameters, using the Sobol indices described in Section 3.4. We focus only on RCE which gives the best fit to the observations. Given that Sobol indices are derived from a small ensemble, potential under-sampling issues may arise in their computation. To assess the sensitivity of the results, we recalculated these indices using random samples of 100 members from the total of 120. The absolute values exhibited variations within the range of 0–3%, yet the magnitudes and conclusions remained consistent.

Figure 9 shows the evolution of the Sobol indices between 1985 and 2019 computed for the volume, the total mass loss and the ice discharge. For the volume, during the first 5 years, the dominant parameters to explain its variability, i.e. corresponding to the highest Sobol indices, are the regularisation coefficient λ _reg and the relaxation time t _relax. By definition, our initialisation procedure does not control the volume drift during the relaxation period. Because t _relax has a large influence on the initial volume compared to the other parameters, this means that the variability in the drift is smaller than the effect of integrating the model forward in time few more years. However, its importance on the volume decreases through time and is very small for the total mass loss, meaning that the response of the system to the front perturbations is not very sensitive to its initial volume. As for λ _reg, it controls the quality of the fit to the observed velocity, so we found natural that it has a large influence on the results. Too much regularisation will tend to smooth the velocity gradients and thus decrease the largest velocity and discharge. For the moment, this regularisation is imposed under the form of a smoothness constraint by penalising the first derivatives and there are no real physical justifications for this particular choice. Improving our understanding of the physical processes controlling the friction, should help to better define priors, and thus reduce the uncertainty in the initial friction fields.

Figure 9. Evolution of Sobol index of RCE through time for the volume (a), the total mass loss since 1985 (b) and the ice discharge (c).

After 1990, the importance of parameters that influence the initial surface elevation weakens, with increasing sensitivity to parameters influencing the ice dynamics such as the exponent of the friction law m (+25%), the enhancement factor E (+14%), the regularisation weight λ _div and the velocity field used for the friction inversion u _obs. This suggests that the variability of the volume, which is initially strongly related to the initial surface elevation, becomes relatively independent of this initial surface later in the simulation. It is then mainly influenced by the parameters directly affecting the ice dynamics (friction coefficient field, m and E) because of their impact on the ice discharge.

Observing the Sobol indices for ice discharge, the significance of the ice dynamics and friction coefficient field parameters becomes evident. Thus, the initial variability of the ice discharge is initially dominated by λ _reg, m and E, and to a lesser degree by λ _div and u _obs. This reflects the major influence of the friction on the overall ice dynamics. The increase in the influence of u _obs after 2008 shows the sensitivity of the model to the initial friction field when changes of front position are applied.

For the mass loss, we observe that the friction and rheology parameters m, E and λ _reg have a larger impact than the forcing parameter (SMB _model) or parameters influencing the initial surface elevation (t _relax, SMB _date). This indicates that the initial surface has little impact on the final mass loss, which is the most interesting variable in terms of sea-level rise. However, as we do not take into account the feedback between SMB and elevation in the model, this sensitivity to SMB could be underestimated.

4.2.2 Performance metrics

Analysing the Sobol indices for RMSE_u and RMSE_zs allows us to identify the parameters that significantly influence the variability in the model's capacity to reconstruct observations. As shown in Figure 10, the enhancement factor E is the parameter that has the most influence on RMSE_u, explaining alone 33% of its variance. Figure 11 shows the distribution of the RMSE_u as a function of E for RCE. It highlights that the maximum values are found for E <1 or >3.5. On the contrary, members with E between 1 and 3.5 perform significantly better. Given the dominant effect of E, no other significant optima are found in the remaining parameters in relation to RMSE_u.

Figure 10. Sobol index of RCE for the RMSE_u (a) and the RMSE_zs (b) calculated on the full period and the full validation area (see Fig. 1).

Figure 11. RMSE_u as a function of E, with each red dot representing a member of RCE and boxplots of RMSE_u for 5 subgroups of E values: between 0.5 and 1.4, between 1.4 and 2.3, between 2.3 and 3.2, between 3.2 and 4.1, between 4.1 and 5.

The sensitivity analysis on RMSE_zs (Fig. 10) indicates that the choice of SMB _model is an important factor of the variability of RMSE_zs. We found that members forced with SMB _RACMO perform better, with an average RMSE_zs of 28.3 m than members forced with SMB _MAR, which have an average RMSE_zs of 31.3 m. The largest differences are found in slow-flowing areas (Fig. 12), where changes in elevation are mostly controlled by the SMB. RACMO has less melt in the slow flowing areas than MAR. Nevertheless, the choice of SMB _model does not seem to be a key point in the variability of the total volume or its evolution as we have seen previously (Fig. 9).

Figure 12. RMSE of the mean of RCE members using MAR (a) or RACMO (b) for each gridcell over the overall period (1985–2019). The grey contour indicates the ensemble mean RCE values, averaged over the period 1985–2019, for 30 m a⁻¹.

5. Discussion

This study validates an ensemble ice-sheet model initialised via inverse methods with prescribed front positions, employing an extensive dataset of surface elevation, ice velocity, ice discharge and mass loss time series. Two distinct ensembles, employing different friction laws, undergo a systematic comparison. The RCE, which parameterises the evolution of effective friction based on the distance to the front to capture lower effective pressure and basal stress closer to flotation, consistently outperforms the WE formulation in reproducing 2-D time series of surface velocity and surface elevation. This improved performance in terms of spatial fields is also reflected in the global variables such as the ice discharge or mass loss. Our ensemble modelling approach successfully alleviates overconfidence, with the RCE adeptly capturing observational data, including velocity, mass loss or ice discharge, along with their estimated variability, thanks to its spread (Figs 3, 6, 7). The sensitivity analysis highlights the dominant impact of friction and ice rheology within the RCE approach, whereas initial surface elevation's influence remains secondary (Figs 9, 10). A rational reduction of the parameters, centred on the significant parameters such as E, m, λ _reg and u _obs, is recommended, as the variance of these four parameters represents 90% of the variance of the total mass loss, while the parameters affecting only the initial conditions of the surface could be excluded due to their minimal impact. Moreover, u ₀ and λ _div display lower impacts compared to these four parameters affecting friction (E, m, λ _reg, u _obs). Subsequent studies stand to benefit from a streamlined parameter space, enhancing exploration and potentially revealing nuanced relationships between performance metrics, global variables and remaining parameters.

One major limitation in our results arises from the comparison with surface elevation. Both ensembles appear overconfident, exhibiting discrepancies with data that fall outside the modelled range (see Fig. 6 and the discussion below for more details). Thinning rates arise from the difference between SMB and flux divergence. Discrepancies between the model and observations may result from inadequately addressed uncertainties in either of these terms. The biases in the regional climate models, possibly linked to their use of a constant surface elevation (Fettweis and others, 2017; Noël and others, 2018), may contribute to discrepancies in SMB. This could result in an underestimation of melting as surface elevation decreases. This can also be partially explained by assuming the absence of basal melting, which, nonetheless, remains relatively low compared to the SMB, specifically <10% in grounded areas (Karlsson and others, 2021). Additionally, overconfidence in flux divergence may arise from uncertainty in bed topography, which is based on sparse measurements and reconstructed using mass conservation principles (Morlighem and others, 2017). For instance, BedMachine version 5 displays lower bed elevations, ~50 m less in specific sections of the UI-N stream compared to version 3. This deeper bed could lead to increased ice discharge, as evidenced by the disparity between Mouginot's data and that of King and Mankoff (Fig. 3), resulting in greater mass loss and a more significant elevation decrease. Finally, this overconfidence might be influenced by the absence of spatial variability in other model parameters.

Comparing our study's findings with previous research on UI history reveals similar trends. Observational studies such as Khan and others (2013) and Larsen and others (2016), or the model study by Haubner and others (2018) demonstrate an increase in UI's annual mass loss starting from 2005, attributed to the retreat of UI-N. While direct comparisons of mass loss values are hindered by differences in catchments, these studies also indicate that changes in elevation and mass loss are predominantly driven by dynamical processes rather than alterations in SMB. Consequently, the inadequacy of our mass change representation can be attributed to our ice-flow model, rather than to SMB or the absence of elevation feedback. This discrepancy may be associated with a deficient representation of the bed, such as a deeper bed for UI-N in that region, as indicated by aircraft flight lines or the newer version of BedMachine (v5) (Morlighem, 2022). This aligns with the outcomes of our sensitivity analysis, wherein we highlighted that the most influential parameters during the historical period are linked to dynamics rather than external forcings. In the study by Khan and others (2013), a correlation is established between calving events of UI-N and periods of warmer air and water. Instead of implicitly addressing processes at the front by forcing the position, a valuable future step in the historical case study of UI would involve incorporating a calving law and a parameterisation of melt dependent on ocean forcing. Evaluating these laws based on their ability to reproduce front position changes, as done in other studies (e.g. Bondzio and others, 2018), would further enhance the analysis. In this context, the four distinct front dynamics observed in UI provide a compelling case study for evaluating the model's skill in replicating these diverse behaviours, thereby strengthening the robustness of the results.

Among the ice dynamic parameters, the flow law, particularly the significance of E, stands out. The uncertainty in the initial temperature field is challenging to estimate due to the prolonged response time of the thermal state, requiring extensive spin-ups. Mechanical coupling introduces significant complexities, hindering the initialisation of a consistent thermo-dynamical steady state (Schäfer and others, 2014). Alternatively, an iterative procedure, as demonstrated in Zhao and others (2018), or an inversion for both friction and ice rigidity, as explored in Ranganathan and others (2021), could have been employed. However, the existence of multiple solutions renders the problem particularly ill-posed and dependent on prior information (Arthern and Gudmundsson, 2010). It can be noted that better constraining the uncertainties on the ice rigidity would also allow for more robust friction-only inversions as shown by Babaniyi and others (2021). Here we have tried to take some of this into account by performing an ensemble of inversions for different values of E. The model seems to reproduce the glacier's acceleration or 1985 state without requiring processes related to thermo-mechanical coupling. However, several processes that could impact viscosity, such as thermo-mechanical feedback, cryo-hydrologic warming, damage, etc., have been omitted. Nonetheless, this could explain the poor replication of the velocities in our study in the shear margins, as studies have shown that phenomena such as damage or shear heating may have contributed to the changes in dynamics observed in certain ice flows (Van Der Veen and others, 2011; Bondzio and others, 2017). Our findings reveal that the model encounters challenges in reproducing velocities within the shear margins. It might also stem from other uncertainties in the flow law (i.e. the spatial variation in the ice rigidity field), but due to the high sensitivity to friction we cannot conclude here on the importance of these processes. Similarly, this misrepresentation might be traced back to input data, particularly as the computation of ice-flow velocities within shear margins poses a great challenge in satellite imagery-based methods (Millan and others, 2019).

The important impact of friction in reproducing the observation of the recent past reinforces the relevance of defining an appropriate friction law. While prior investigations have delved into constraining the form of friction laws using multiple velocity fields (Joughin and others, 2004; Gillet-Chaulet and others, 2016; Choi and others, 2022), our ensemble modelling approach surpasses the extent of investigation previously conducted in this regard. On the one hand, Gillet-Chaulet and others (2016) and Hillebrand and others (2022) have shown that, when considering Weertman or Budd laws, the exponent m should surpass 5 to effectively replicate the dynamics in rapidly flowing regions. On the other, Joughin and others (2019) proposed that a regularised Coulomb law, with a friction decreasing towards 0 at the grounding line, could effectively reproduce observed velocity changes between 2002 and 2017 in Pine Island Glacier. Choi and others (2022) showed that friction laws that include a parameterised dependence on the effective pressure better reproduce the observed acceleration and mass loss of the past decade in northwest Greenland. Our results with WE on UI-N, where the front has retreated the most, is coherent finding with these previous studies: the magnitude of the velocity changes in response to the retreating front is underestimated even for the lowest exponents of the friction law m when using a Weertman friction law. These members with m lower than 1/3 leads to a larger acceleration, but this change seems far too small compared to a switch to a regularised Coulomb law associated with a parameterisation of the effective pressure as suggested by Joughin and others (2019).

The regularised-Coulomb friction law takes into account the effective pressure in a parameterised way by reducing the friction in the first 15 km upstream of the ice front fast flow areas. The results with RCE show that this parameterisation allows to reproduce much better the surface velocities observed before the front retreat, yet at point A it only corresponds to an increase of the basal shear stress of the order of 20 kPa (Fig. 6). Assuming a constant water pressure such a variation would correspond to a change in ice thickness of <2 m whereas the observations show changes in altitude of more than 150 m over the period. As a consequence, a parameterisation where the friction coefficient would change linearly with the thickness above flotation as done in Joughin and others (2019) would lead to too large variations. Furthermore, the UI-S results also show that the velocities can be quite sensitive to the conditions at the front, which makes the interpretation of the results relatively difficult, as velocity variations can be the result of changes that happen locally or at distances of several ice thicknesses.

Because the basal friction field has been calibrated with recent observations, our methodology inherently yields distinct initial states for the WE and RCE ensembles. Notably, RCE leads to initial velocity distributions more in line with observational data (Figs 4, 5). The observed performance disparities predominantly stem from the inadequate initial state of WE, as depicted in Figure 8. This poor performance in terms of initial velocity and elevation fields could partly be attributed to the different methods for the extrapolation of the friction coefficients in the areas that have not been constrained by the inversion. To confirm that the superior performance of the RCE can be attributed to the disparity in friction laws rather than differences in initial states, we conducted a sensitivity experiment. This involved an ensemble using Weertman's friction law with a fixed friction field closer to the 1985 observations, denoted as WE85. For this ensemble, the initial friction field is derived through a procedure identical to that of the RCE, as outlined in the Supplementary materials. The results are shown and compared to WE and RCE in Supplementary Figures S6–S9. Because WE and RCE have been calibrated using recent velocity observations, both ensembles had relatively similar performances in terms of velocity once the fronts have reached a position close to that of 2010 (Figs 6, 7). On the contrary, WE85 and RCE have similar velocities before the front retreat, but the difference between the two ensembles increases when the retreats initiate and WE85 underestimates the more recent velocity observations. The differences in terms of mass loss are less striking because the discharges predicted by the two ensembles really start to diverge after 2005. These results reinforce our conclusion that it is not possible to reproduce all the observations with a friction law that does not take into account changes in the basal friction. This is effectively achieved in RCE in a parameterised way introducing the distance to the front, and progress are needed in our understanding of the basal processes for a physical modelling of the evolution of the basal conditions.

As detailed in the Supplementary materials, our parameterisation is built on the observation that friction coefficients derived from inversion beneath the central flowlines of UI-N, UI-C and UI-S exhibit a linear correlation with the distance to the glacier front. This correlation is especially pronounced within the initial few kilometres (Supplementary Fig. S1). This correspondence can be linked to the linear association between the height above flotation (h _af) and the distance to the front along these specific flow paths. This finding aligns with the outcomes of Habermann and others (2013), who established a linear relationship between h _af and β _RC (referred to as τ _c in their work) for the most rapid flowing parts of Jakobshavn Isbræ. Their study also revealed analogous temporal variations in both β _RC and h _af, albeit with β _RC displaying more localised fluctuations within the flow. The proposed parameterisation framework adeptly captures these identified characteristics (Supplementary Fig. S2).

Nonetheless, we chose to maintain the reliance on the dependence on distance from the glacier front rather than adopting the more conventional approach of using ice elevation above flotation. Several studies have alternatively used variables like height above the flotation level (Vieli and others, 2001; Bondzio and others, 2017; Haubner and others, 2018) or relationships involving ice thickness and bed elevation (Downs and Johnson, 2022). However, these associations are notably susceptible to variations in surface and bed elevations. Especially within the realm of historical reconstructions, surface elevation data often presents constraints such as sparsity, incompleteness or inaccuracy, in contrast to the more abundant availability of front data (Wood and others, 2021). Consequently, our chosen parameterisation is particularly advantageous for historical inquiries, where data on the distance to the glacier front or the grounding line for ice shelves extend further into the past compared to surface elevation data, which is predominantly accessible from the 1980s onwards.

Looking ahead, the parameterisation we are proposing is likely to have a substantial influence on both predictions of future sea-level rise and confidence in these projections. By demonstrating the model's capacity, augmented by our parameterisation, to faithfully replicate a 35 year interval of velocity and elevation data, we underscore the potential of inverse method models to effectively emulate short-term data dynamics contingent upon an appropriate friction law. This aligns with the implications outlined in Barnes and Gudmundsson (2022), affirming the predictive skill of ice-sheet models within temporal spans, despite the inherent uncertainties associated with basal slip mechanisms. Moreover, Brondex and others (2019) underline the consistent underestimation of mass losses predicted by Weertman's law. Given this context, the adoption of our parameterisation to diminish friction near the glacier front, contrasts with the Weertman law characterised by a constant friction coefficient field over time. Specifically, as the glacier front retreats in future scenarios, friction coefficients within the glacier will decrease for RCE in comparison to WE, leading to an augmented glacier acceleration, higher discharge rates and ultimately an increased mass loss. Implementing such a parameterisation on a larger scale of the entire ice sheet therefore suggests the possibility of an amplified contribution of the GrIS to global sea-level rise.

Concluding on the aspect of validation, our investigation underscores the potential of incorporating spatio-temporal velocity data to validate polar ice-sheet models, thereby extending the temporal scope for model validation. While previous research has also aimed to validate ice-sheet models using satellite data (Aschwanden and others, 2013; Price and others, 2017; Nias and others, 2023), the use of altimetry or gravity data has constrained validation to approximately a decade or slightly more. By using observed velocities here, we greatly enlarge this window with data over a 35 year period. In comparison, the altimetry data collected are much more sparse in time and space, although still extends back to the 2000s (i.e. 23 years long time series, Supplementary Fig. S4). This may raise the validation issue raised by Goelzer and others (2018), i.e. that these techniques currently suffer the same limitations of observational satellite datasets with short time coverage. The use of declassified spy satellite observations from the 1970s and 1990s, such as the KH-9 Hexagon Satellite Images mission (Dehecq and others, 2020), has the potential to considerably extend the observation period used in our study (1985–2019) and enhance the validation of ice-sheet models.

Our results highlight the need to develop and use spatio-temporal velocity and elevation series at the scale of the GrIS to validate models using inverse methods. This validation then provides a better understanding of the behaviour of the global ice sheet and greater confidence in the projections made by this type of ice-sheet model with all its assumptions. Furthermore, our parameterisation should also be tested under different conditions. To achieve this, validation on other Greenland glaciers or at the scale of the ice sheet appears to be the most reliable, considering the impossibility of temporal validation for UI due to limited data availability before 1985.

6. Conclusion

In this study, we performed several ensemble simulations to model the dynamic evolution of UI by forcing the evolution of the front position. By evaluating the ensemble model against the time series of surface elevation and ice velocity, we found that the observations are sufficiently robust to help constrain the model. To reproduce these data observations, the model must include a reduction in friction in the stream and an appropriate way to extrapolate friction in front retreat areas as proposed here with RCE. In the case of UI, where the front was 6 km further ahead in 1985, this results in a 20 kPa increase in basal stress at the stream 5.5 km from the current front, representing a 200% increment (from 10 to 30 kPa). Our ensemble approach facilitates the assessment of a range of simulations, incorporating uncertainties in model parameters, in comparison with observations. This process enhances confidence in the model's capability to accurately replicate past and, consequently, future changes in ice dynamics. It is very likely that this choice of friction law has an impact on the projections of the future evolution of marine-terminated glaciers in Greenland. With a friction law that reduces near-front friction, it is likely that the future contribution of glacier dynamics to mass loss will be greater than for a friction law without reduction as has been done so far (Åkesson and others, 2021).

The sensitivity analysis, made possible by the ensemble approach, showed the main role of the initial friction field for a tidewater glacier where strong changes in dynamics and geometry occur. The other sources of uncertainty, related to the initial surface or the SMB, have on the contrary a lower influence on this historic reconstruction (ice discharge, ice mass loss, velocity and elevation).

These results should be reinforced in the future on this case study by performing a projection of UI with a similar approach to see its contribution to sea-level rise. It will be interesting to compare the impact of the sources of uncertainty identified here, such as the friction law or the flow law, with those related to atmospheric and oceanic forcings or shared socioeconomic pathways. Furthermore, extending this study to other study cases or to the scale of Greenland can overcome potential biases related to the study of this glacier and highlight whether the model performance remains the same.