Assimilation of surface observations in a transient marine ice sheet model using an ensemble Kalman filter

Marine-based sectors of the Antarctic Ice Sheet are increasingly contributing to sea level rise. The basal conditions exert an important control on the ice dynamics and can be propitious to instabilities in the grounding line position. Because the force balance is non-inertial, most ice flow models are now equipped with time-independent inverse methods to constrain the basal conditions from observed surface velocities. However, transient simulations starting from this initial state usually suffer from inconsistencies and are not able to reproduce observed trends. Here, using a synthetic flow line experiment, we assess the performance of an ensemble Kalman filter for the assimilation of transient observations of surface elevation and velocities in a marine ice sheet model. The model solves the shallow shelf equation for the force balance and the continuity equation for ice thickness evolution. The position of the grounding line is determined by the floatation criterion. The filter analysis estimates both the state of the model, represented by the surface elevation, and the basal conditions, with the simultaneous inversion of the basal friction and topography. The idealised experiment reproduces a marine ice sheet that is in the early stage of an unstable retreat. Using observation frequencies and uncertainties consistent with current observing systems, we find that the filter allows the accurate recovery of both the basal friction and topography after few assimilation cycles with relatively small ensemble sizes. In addition it is found that assimilating the surface observations has a positive impact on constraining the evolution of the grounding line during the assimilation window. Using the initialised state to perform century-scale forecast simulations, we show that grounding line retreat rates are in agreement with the reference; however remaining uncertainties in the basal conditions may lead to significant delays in the initiation of the unstable retreat. These results are encouraging for the application to real glacial systems.

1 Introduction

Despite recent significant improvements in ice sheet models, the projected magnitude and rate of the Antarctic and Greenland ice sheets' contribution to 21st century sea-level rise (SLR) remains poorly constrained (Church et al., 2013). Improving our ability to model the century-scale magnitude and rates of mass loss from marine ice sheets remains a key scientific objective (Scambos et al., 2017).

Improving SLR estimates requires, amongst other things, correctly modelling the dynamics of the grounding line (GL), i.e. the location where the ice detaches from its underlying bed and goes afloat on the ocean (Durand and Pattyn, 2015). In the GL vicinity, the stress regime changes from a regime dominated by vertical shearing in the grounded part to a buoyancy-driven flow dominated by longitudinal stretching and lateral shearing (Pattyn et al., 2006; Schoof, 2007). Because this transition occurs on horizontal dimensions that are smaller than the typical grid size of large-scale ice sheet models, many studies have focussed on the ability of the numerical model to properly simulate grounding line migration using synthetic experiments (e.g. Vieli and Payne, 2005; Durand et al., 2009; Gladstone et al., 2012; Seroussi et al., 2014). Two Marine Ice Sheet Model Intercomparison Projects (MISMIP) have allowed the identification of the minimum requirements to properly resolve GL motion: (i) inclusion of membrane stresses and (ii) a sufficiently small grid size or a subgrid interpolation of the GL (Pattyn et al., 2012, 2013). These results suggest that, in realistic applications, the numerical error could be reduced below the errors associated with uncertainties in the initial model state, in the model parameters, and in the forcings from the atmosphere and ocean.

For obvious reasons of inaccessibility, the basal conditions (topography and friction) are an important source of uncertainties. Because of the intrinsic instability of marine ice sheets resting over a seaward up‐sloping bed, the resolution of the bed topography in the coastal regions can significantly affect short-term ice sheet forecasts (Durand et al., 2011; Enderlin et al., 2013). Analytical developments have shown that the flux at the grounding line depends on the friction law and its coefficients (Schoof, 2007; Tsai et al., 2015). The sensitivity of model projections to the basal friction has been confirmed by several numerical studies on both synthetic and real applications (Joughin et al., 2010; Ritz et al., 2015; Brondex et al., 2017, 2019). In particular, Brondex et al. (2017) have shown that, for unbuttressed ice sheets, spatially varying friction coefficients can also lead to stable GL positions in up‐sloping bed regions.

Uncertainties in the model state and parameters can be reduced by data assimilation (DA). The objective of formal DA methods is to update the model using observations in a framework consistent with the model, the data and their associated uncertainties (Bannister, 2017). Most ice flow models are now equipped with variational methods to constrain the basal conditions from surface observations (e.g. MacAyeal, 1993; Vieli and Payne, 2003; Larour et al., 2012; Gillet-Chaulet et al., 2012). However most studies perform “snapshot” calibrations, where the inversion is performed at a unique initial time step. The state of the model produced from this calibration is therefore sensitive to inconsistencies between the different datasets. The resulting transient artefacts are usually dissipated during a relaxation period where the model drifts from the observations.

Because historic remote sensing data collections are spatially incomplete as well as temporally sparse, most distributed maps are mosaicked, stacked or averaged to maximise the spatial coverage at the expense of the temporal information (Mouginot et al., 2012). However, in the last few years, the development of space-borne ice sheet observations has entered a new era with the launch of new satellite missions, considerably increasing the spatial and temporal resolution of surface observations. Because they require linearised versions of the forecast model and of the observation operator, extending the existing variational methods implies important numerical developments (e.g. Goldberg et al., 2016; Larour et al., 2016; Hascoët and Morlighem, 2018). In Goldberg and Heimbach (2013), a time-dependent adjoint ice flow model is derived using a source-to-source algorithmic differentiation software combined with analytical methods. The DA capabilities are illustrated with a suite of synthetic experiments, including the simultaneous inversion of the basal topography and friction from surface observations and the assimilation of transient surface elevations to retrieve initial ice thicknesses. In a real-world application to a region of West Antarctica, they show that assimilating annually resolved observations of surface height and velocities between 2002 and 2011 allows the improvement of the initial model state, giving better confidences in projected committed mass losses (Goldberg et al., 2015). Because of the complexity of the code, Larour et al. (2014) use an operator-overloading approach to generate the adjoint and assimilate surface altimetry observations from 2003 to 2009 to constrain the temporal evolution of the basal friction and surface mass balance of the Northeast Greenland Ice Stream.

Ensemble DA methods, based on the ensemble Kalman filter (EnKF), have been successful in solving DA problems with large and non-linear geophysical models. Comparative discussions of the performances and advantages of variational and ensemble DA methods can be found in, e.g. Kalnay et al. (2007), Bannister (2017) and Carrassi et al. (2018). As they aim at solving similar problems, a recent tendency is to combine both methods to benefit from their respective advantages.

EnKF approximates the state and the error covariance matrix of a system using an ensemble that is propagated forward in time with the model, avoiding the computation of the covariance matrices and the use of linearised or adjoint models. Contrary to time-dependent variational methods where the objective is to find the model trajectory that minimises the difference with all the observations within an assimilation window, EnKF assimilates the observations sequentially in time as they become available using the analysis step of the Kalman filter, as illustrated in Fig. 1. The model trajectory is then discontinuous and, at a given analysis, the model is only informed by past and present observations. For the retrospective analysis of a time period in the past, i.e. a reanalysis, ensemble filters can easily be extended to smoothers to provide analyses that are informed by all past, present and future observations (Evensen and van Leeuwen, 2000; Li and Navon, 2001; Cosme et al., 2012; Nerger et al., 2014). Since the first version introduced by Evensen (1994), many variants have been developed, mainly differing in the way the Kalman filter analysis is rewritten and the analysed error covariance matrix is resampled (e.g. Burgers et al., 1998; Houtekamer and Mitchell, 1998; Pham et al., 1998; Bishop et al., 2001; Nerger et al., 2012). A review of the most popular EnKFs using common notations can be found in Vetra-Carvalho et al. (2018). Efficient and parallel algorithms have been developed, and because they are independent of the forward model, several open-source toolboxes that implement various EnKFs are now available, e.g. OpenDA (https://www.openda.org, last access: 25 June 2018) and PDAF (http://pdaf.awi.de, last access: 25 June 2018).

Figure 1 Principle of data assimilation (adapted from Carrassi et al., 2018). Having a physical model able to forecast the evolution of a system from time t=t₀ to time t=T_f (cyan curve), the aim of DA is to use available observations (blue triangles) to correct the model projections and get closer to the (unknown) truth (dotted line). In EnKFs, the initial system state and its uncertainty (green square and ellipsoid) are represented by N_e members. The members are propagated forward in time during n₁ model time steps dt to t=T₁ where observations are available (forecast phase, orange dashed lines). At T=t₁ the analysis uses the observations and their uncertainty (blue triangle and ellipsoid) to produce a new system state that is closer to the observations and with a lower uncertainty (red square and ellipsoid). A new forecast is issued from the analysed state and this procedure is repeated until the end of the assimilation window at t=T_f. The model state should get closer to the truth and with lower uncertainty as more observations are assimilated. Time-dependent variational methods (4D-Var) iterate over the assimilation window to find the trajectory that minimises the misfit (J₀) between the model and all observations available from t₀ to T_f (violet curve). For linear dynamics, Gaussian errors and infinite ensemble sizes, the states produced at the end of the assimilation window by the two methods should be equivalent (Li and Navon, 2001).

As Monte Carlo methods, EnKFs suffer from under-sampling issues as often the size of the ensemble is much smaller than the size of the system to estimate. Localisation and inflation are popular methods to counteract these issues and to increase the stability of the filtering. Because they are based on the original Kalman filter equations, EnKFs are optimal only for Gaussian distributions and linear models. However, the many applications in geoscience with large and non-linear models have shown that the method remains robust in general and EnKFs are used in several operational centres with atmosphere, ocean and hydrology models (e.g. Sakov et al., 2012; Houtekamer et al., 2009; Hendricks Franssen et al., 2011). While firstly developed for numerical weather and ocean prediction where the forecasts are very sensitive to the model initial state, the method is also widely used, e.g. in hydrology, for joint state and parameter estimations (Sun et al., 2014).

In the context of ice sheet modelling, encouraging results have been obtained by Bonan et al. (2014) for the estimation of the state and basal conditions of an ice sheet model using the ensemble transform Kalman filter (ETKF; Bishop et al., 2001; Hunt et al., 2007). They study the performance of the method using idealised twin experiments where perturbed observations generated from a model run are used in the DA framework to retrieve the true model states and parameters. Using a flow line shallow ice model, they show that both the basal topography and basal friction can be retrieved with good accuracy from surface observations with realistic noise levels, even for relatively small ensembles. The method has been further developed to assimilate the margin position in a shallow ice model that explicitly tracks the boundaries with a moving mesh method (Bonan et al., 2017).

The purpose of this paper is to explore the performance of ensemble Kalman filtering for the initialisation of a marine ice sheet model that includes GL migration. In particular, we want to address (i) the quality of the analysis for the simultaneous estimation of the basal topography and friction in the context of a marine ice sheet that is undergoing an unstable GL retreat and (ii) the effects of the remaining uncertainties for the predictability of GL retreat. The ice flow model and the EnKF used in this study are described in Sect. 2. To test the DA framework, we define a twin experiment in Sect. 3. Section 4 presents the results for both the transient assimilation and the forecasts. Finally, perspectives and challenges for real applications are discussed in Sect. 5, before concluding remarks.

2 Methods

2.1 Ice flow model

The gravity-driven free surface flow of ice is solved using the finite-element ice flow model Elmer/Ice (Gagliardini et al., 2013).

For the force balance, we solve the shelfy stream approximation (SSA) equation (MacAyeal, 1989) in one horizontal dimension. This is a vertically integrated model that derives from the Stokes equations for small aspect ratio and basal friction. In 1D, this leads to the following non-linear partial differential equation for the horizontal velocity field u:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}H{\displaystyle \frac{\partial u}{\partial x}}\right)-{\mathit{\tau }}_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial x}},\end{array}$$

with ρ_i the ice density, g the gravity norm, H=z_s–z_b the ice thickness, and z_s and z_b the top and bottom surface elevations, respectively. Using Glen's constitutive flow law, the vertically averaged effective viscosity $\stackrel{\mathrm{‾}}{\mathit{\eta }}$ is given by

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{‾}}{\mathit{\eta }}={\displaystyle \frac{\mathrm{1}}{H}}\underset{z_{\mathrm{b}}}{\overset{z_{\mathrm{s}}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}{D}_{\mathrm{e}}^{(\mathrm{1}-n)/n}\text{d}z,\end{array}$$

where D_e is the second invariant of the strain-rate tensor, here equal to $D_{\mathrm{e}}^{\mathrm{2}}=(\partial u/\partial x{)}^{\mathrm{2}}$, A is the rate factor and n is the creep exponent, taken equal to the usual value n=3 in the following. The basal friction τ_b is null under floating ice and is represented by the non-linear Weertman friction law for grounded ice:

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{b}}=C{u}^m,\end{array}$$

with C and m the friction coefficient and exponent, respectively. In the following, we use the classical power law with $m=\mathrm{1}/n=\mathrm{1}/\mathrm{3}$. When in contact with the ocean, the ice is assumed to be in hydrostatic equilibrium. The floating condition is evaluated directly at the integration points, and τ_b in Eq. (1) is set to 0 wherever ice is floating (Seroussi et al., 2014).

The time dependency is introduced by the evolution of the top and bottom free surfaces. Because of the hydrostatic equilibrium, the ice sheet topography is fully defined by the bed elevation b and only one prognostic variable. Equation (1) is then coupled with the vertically integrated mass conservation equation for the evolution of the ice thickness H:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(uH\right)}{\partial x}}=a_{\mathrm{s}}-a_{\mathrm{b}},\end{array}$$

with a_s the surface accumulation rate and a_b the basal melt rate. The free surfaces z_s and z_b are obtained from the floating condition which, for z_s, using a constant sea level z_sl=0, gives

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{ll}{z}_{\mathrm{s}}=b+H& \text{for}H\ge -b\frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}\\ z_{\mathrm{s}}=H\left(\mathrm{1}-\frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}\right)& \text{otherwise}\end{array}\right.\end{array}$$

with ρ_w the sea water density.

2.2 Data assimilation

2.2.1 Filter algorithm

For the assimilation, we use the error subspace ensemble transform Kalman filter (ESTKF; Nerger et al., 2012). Originally derived from the singular evolutive interpolated Kalman filter (SEIK; Pham et al., 1998), ESTKF leads to the same ensemble transformations as the ETKF but at a slightly lower computational cost. In practice we use the local version of the filter implemented in PDAF (http://pdaf.awi.de, last access: 25 June 2018; Nerger et al., 2005b) and coupled to Elmer/Ice in an offline mode. This section outlines the ESTKF algorithm.

As an EnKF, ESTKF approximates the state x^k and the error covariance matrix P_k of a system at time t_k using an ensemble of N_e realisations ${\mathit{x}}_i^k$, $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{e}}$. The state vector, of size N_x, contains the prognostic variables and model parameters to be estimated and is approximated by the ensemble mean,

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{‾}}{\mathit{x}}}^k={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}{\mathit{x}}_i^k,\end{array}$$

while the error covariance matrix is approximated by

$$\begin{array}{}\text{(7)}& {\mathbf{P}}_k={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}-\mathrm{1}}}{{\mathbf{X}}^{\prime }}_k{{\mathbf{X}}^{\prime }}_k^T,\end{array}$$

where ${{\mathbf{X}}^{\prime }}_k=({\mathit{x}}_{\mathrm{1}}^k-{\stackrel{\mathrm{‾}}{\mathit{x}}}^k,\mathrm{\dots },{\mathit{x}}_{N_{\mathrm{e}}}^k-{\stackrel{\mathrm{‾}}{\mathit{x}}}^k)\in {\mathbb{R}}^{N_x\times N_{\mathrm{e}}}$ is the ensemble perturbation matrix.

The algorithm can be decomposed in two steps, the forecast and the analysis. Superscripts f (a) denote quantities related to each step. The forecast propagates the state and the error covariance matrix of the system forward in time, from a previous analysis at $t=t_{k-\mathrm{1}}$ to the next observation time t=t_k. For this, the numerical model ℳ_k, assumed perfect in the sequel, is used to propagate each ensemble member individually during n_dt model time steps:

$$\begin{array}{}\text{(8)}& {\mathit{x}}_i^{f,k}={\mathcal{M}}_k\left({\mathit{x}}_i^{a,k-\mathrm{1}}\right).\end{array}$$

At t=t_k, a vector of observations ${\mathit{y}}_{\mathrm{o}}^k$ of size N_y (usually with N_y≪N_x) is available. ${\mathit{y}}_{\mathrm{o}}^k$ is related to the true system state x^f by ${\mathit{y}}_{\mathrm{o}}^k=\mathcal{H}\left({\mathit{x}}^f\right)+{\mathit{ϵ}}^k$, where the observation error ϵ^k is assumed to be a white Gaussian distributed process with known covariance matrix R^k, and ℋ is the observation operator that relates the state variables to the observations. When ${\mathit{y}}_{\mathrm{o}}^k$ is the observed surface velocities, the relation between the observations and the system state, i.e., the ice sheet geometry, and parameters, i.e. the boundary conditions, is given by the force balance Eq. (1); thus ℋ is a non-linear elliptic partial differential equation.

The analysis provides a new estimation of the system state by combining the information from the forecast and the observations. In the following we will omit the time index k in the notations as the entire analysis is performed at t=t_k. As other EnKFs, ESTKF uses the Kalman filter update equations to compute the analysed system state ${\stackrel{\mathrm{‾}}{\mathit{x}}}^a$ and covariance matrix P^a from the forecast, the observations and their uncertainties:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathit{x}}}^a={\stackrel{\mathrm{‾}}{\mathit{x}}}^f+\mathbf{K}\mathit{d}\\ {\mathbf{P}}^a=(\mathbf{I}-\mathbf{KH}){\mathbf{P}}^f,\end{array}\right.\end{array}$$

where $\mathit{d}={\mathit{y}}_{\mathrm{o}}-\mathcal{H}\left({\stackrel{\mathrm{‾}}{\mathit{x}}}^f\right)$ is the innovation and K is the Kalman gain given by

$$\begin{array}{}\text{(10)}& \mathbf{K}={\mathbf{P}}^f{\mathbf{H}}^T({\mathbf{HP}}^f{\mathbf{H}}^T+\mathbf{R}{)}^{-\mathrm{1}}.\end{array}$$

Here, H is the linearised observation operator at the forecast mean. However, in practice H does not need to be computed as it always acts as an operator to project the ensemble members in the observation space. Defining the forecast ensemble projected in the observation space by ${\mathit{y}}_i^f=\mathcal{H}\left({\mathit{x}}_{\mathrm{1}}^f\right)$, $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{e}}$ with ${\stackrel{\mathrm{‾}}{\mathit{y}}}^f$ the ensemble mean, we make the linear approximation

$$\begin{array}{}\text{(11)}& {\mathbf{Y}}^f={\mathbf{HX}}^f,\end{array}$$

with ${\mathbf{X}}^f=({\mathit{x}}_{\mathrm{1}}^f,\mathrm{\dots },{\mathit{x}}_{N_{\mathrm{e}}}^f)\in {\mathbb{R}}^{N_x\times N_{\mathrm{e}}}$ the forecast ensemble matrix and ${\mathbf{Y}}^f=({\mathit{y}}_{\mathrm{1}}^f,\mathrm{\dots },{\mathit{y}}_{N_{\mathrm{e}}}^f)\in {\mathbb{R}}^{N_y\times N_{\mathrm{e}}}$ its equivalent in the observation space.

In practice, with large models (N_x>>1), the covariance matrices P^f and P^a of size N_x×N_x can not be formed, so that, to be implemented, the analysis (Eq. 9) needs to be reformulated. Moreover, the sample covariance matrix approximated with an ensemble of size N_e (Eq. 7) is only a low-rank approximation of the true covariance matrix and its rank is at most N_e−1. ESTKF uses this property to write the analysis in a (N_e−1)-dimensional subspace spanned by the ensemble and referred to as the error subspace (Nerger et al., 2005a). The forecast covariance matrix P^f is then rewritten as

$$\begin{array}{}\text{(12)}& {\mathbf{P}}^f={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}-\mathrm{1}}}{\mathbf{LL}}^T,\end{array}$$

where $\mathbf{L}\in {\mathbb{R}}^{N_x\times N_{\mathrm{e}}-\mathrm{1}}$ is given by

$$\begin{array}{}\text{(13)}& \mathbf{L}={\mathbf{X}}^f\mathbf{\Omega }.\end{array}$$

The matrix $\mathbf{\Omega }\in {\mathbb{R}}^{N_{\mathrm{e}}\times N_{\mathrm{e}}-\mathrm{1}}$ defined as

$$\begin{array}{}\text{(14)}& {\mathrm{\Omega }}_{ij}=\left\{\begin{array}{ll}\mathrm{1}-\frac{\mathrm{1}}{N_{\mathrm{e}}}\frac{\mathrm{1}}{\frac{\mathrm{1}}{\sqrt{N_{\mathrm{e}}}}+\mathrm{1}}& \text{for}i=j,i<N_{\mathrm{e}}\\ -\frac{\mathrm{1}}{N_{\mathrm{e}}}\frac{\mathrm{1}}{\frac{\mathrm{1}}{\sqrt{N_{\mathrm{e}}}}+\mathrm{1}}& \text{for}i\ne j,i<N_{\mathrm{e}}\\ -{\displaystyle \frac{\mathrm{1}}{\sqrt{N_{\mathrm{e}}}}}& \text{for}i=N_{\mathrm{e}}\end{array}\right.\end{array}$$

projects the ensemble matrix X^f onto the error subspace. The multiplication with X^f subtracts the ensemble mean and a fraction of the last column of the ensemble perturbation matrix ${{\mathbf{X}}^{\prime }}^f$ from all other columns.

After some algebra using Eqs. (12) and (9), P^a can be written as a transformation of L,

$$\begin{array}{}\text{(15)}& {\mathbf{P}}^a={\mathbf{LAL}}^T,\end{array}$$

with the transform matrix $\mathbf{A}\in {\mathbb{R}}^{N_{\mathrm{e}}-\mathrm{1}\times N_{\mathrm{e}}-\mathrm{1}}$ given by

$$\begin{array}{}\text{(16)}& {\mathbf{A}}^{-\mathrm{1}}=\mathit{\rho }(N_{\mathrm{e}}-\mathrm{1})\mathbf{I}+({\mathbf{Y}}^f\mathbf{\Omega }{)}^T{\mathbf{R}}^{-\mathrm{1}}\mathbf{Y}\mathbf{\Omega },\end{array}$$

where $\mathit{\rho }\in [\mathrm{0},\mathrm{1}]$ is the forgetting factor discussed in Sect. 2.2.2.

Finally, the update step is obtained as a single equation for the transformation of the forecast ensemble X^f to the analysed ensemble X^a as

$$\begin{array}{}\text{(17)}& {\mathbf{X}}^a={\stackrel{\mathrm{‾}}{\mathbf{X}}}^f+{\mathbf{X}}^f\mathbf{\Omega }(\stackrel{\mathrm{‾}}{\mathbf{W}}+\mathbf{W}),\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathbf{X}}}^f$ is the matrix where the columns are given by the forecast ensemble mean, $\stackrel{\mathrm{‾}}{\mathbf{W}}$ is a matrix where the columns are given by the vector

$$\begin{array}{}\text{(18)}& \stackrel{\mathrm{‾}}{\mathit{w}}=\mathbf{A}\left(\mathbf{Y}\mathbf{\Omega }{)}^T{\mathbf{R}}^{-\mathrm{1}}\right({\mathit{y}}_{\mathrm{o}}-{\stackrel{\mathrm{‾}}{\mathit{y}}}^f),\end{array}$$

and W is given by

$$\begin{array}{}\text{(19)}& \mathbf{W}=\sqrt{N_{\mathrm{e}}-\mathrm{1}}\mathbf{C}{\mathbf{\Omega }}^T,\end{array}$$

where C is the symmetric square root of A obtained by singular value decomposition.

Finally, the analysed ensemble X^a is used as the initial ensemble for the next forecast, and so on up to the end of the data assimilation window.

We draw attention to several remarks on the algorithm.

• To compute the innovation d, we have made the same linear approximation $\mathcal{H}\left({\stackrel{\mathrm{‾}}{\mathit{x}}}^f\right)={\stackrel{\mathrm{‾}}{\mathit{y}}}^f$ as Hunt et al. (2007). This choice is consistent with the computation of the covariance matrices P^fH^T and HP^fH^T in Eq. (10) using the linear approximation Eq. (11) (Houtekamer and Mitchell, 2001).

• Several ensembles can have the same mean and covariance matrix, which is why several EnKFs exactly satisfy Eq. (9) but lead to different ensemble transformations and thus different analysed ensembles (Vetra-Carvalho et al., 2018). With the same arguments several variants of ESTKF can be introduced, e.g. by replacing Ω in Eq. (19) by a random matrix with the same properties or using a Cholesky decomposition to compute C.

• As written here, the ESTKF leads to the same ensemble transformation as the ETKF. However, as the computations are not performed in the same subspace, tiny differences due to the finite precision of the computations may grow, leading to slight differences at the end of the assimilation window (Nerger et al., 2012).

• The leading computational cost of the ensemble transformation in ESTKF is $\mathcal{O}\left(N_y(N_{\mathrm{e}}-\mathrm{1}{)}^{\mathrm{2}}+N_{\mathrm{e}}(N_{\mathrm{e}}-\mathrm{1}{)}^{\mathrm{2}}+N_x{N}_{\mathrm{e}}(N_{\mathrm{e}}-\mathrm{1})\right)$, so it scales linearly with N_x and N_y (Nerger et al., 2012). Naturally, increasing N_e also requires an increase in the number of model runs and, in general, the objective is to get the ensemble size as small as possible. The performance of the algorithm also depends on the evaluation of the product of R⁻¹ with some vectors, which can become more expensive when the observation errors are spatially correlated.

2.2.2 Filter stabilisation: inflation and localisation

In practice for large-scale problems, EnKFs as Monte Carlo methods suffer from under-sampling issues. First, because of the rank deficiency of the covariance matrix P^f, the analysis adjusts the model state only in the error subspace, ignoring error directions not accounted for by the ensemble (Hunt et al., 2007). This can result in an analysis that is overconfident and underestimates the true variances. In the long run, the ensemble spread will become too small and the analysis will give too much weight on the forecast, finally disregarding the observations and diverging from the true trajectory. A common simple ad hoc remedy is to inflate the forecast covariance matrix with a multiplicative factor (Pham et al., 1998; Anderson and Anderson, 1999). Here, inflation has been introduced in Eq. (16) using the forgetting factor $\mathit{\rho }\in [\mathrm{0},\mathrm{1}]$ with ρ=1 corresponding to no inflation (Pham et al., 1998). It is the inverse of the inflation factor used by Bonan et al. (2014).

Second, the rank deficiency of P^f leads to the appearance of spurious correlations between parts of the system that are far away. As these correlations are usually small, a common remedy is to damp these correlations with a procedure called localisation. In covariance localisation, localisation is applied by using an ensemble covariance matrix that results from the Schur product of P^f with an ad hoc correlation matrix that drops long-range correlations (Hamill et al., 2001; Houtekamer and Mitchell, 2001). However, this localisation technique is not practical for square-root filters where P^f is never explicitly computed. Here, as in Bonan et al. (2014), we use a localisation algorithm based on domain localisation and observation localisation (Ott et al., 2004; Hunt et al., 2007). Both methods are illustrated in Sakov and Bertino (2011), who conclude that they yield similar results. Domain localisation assumes that observations far from a given location have negligible influence. In practice, the state vector in each single mesh node is updated independently during a loop through the nodes that can easily be parallelised for numerical efficiency. For each local analysis, only the observations within a given radius r from the current node are used. In addition, to avoid an abrupt cut-off, the observation error covariance matrix R is modified so that the inverse observation variance decreases to zero with the distance from the node using a fifth-order polynomial function which mimics a Gaussian function but has compact support (Gaspari and Cohn, 1999). Because it drops spurious long-range correlations and allows the local analyses to choose different linear combinations of the ensemble members in different regions, localisation implicitly increases the rank of the covariance matrix, leading to a larger dimension of the error subspace, implicitly increasing the effective ensemble size and the filter stability (Nerger et al., 2006; Hunt et al., 2007). However, it has been reported that localisation could produce imbalanced solutions (Mitchell et al., 2002). Here, because the force balance is non-inertial and the SSA assumes that the ice shelves are in hydrostatic equilibrium, this should not be an issue. Another disadvantage is that, when long-range correlations truly exist, the analysis will ignore useful information that could have been used from distant observations.

Here, the forgetting factor ρ and the localisation radius r will be used as tuning parameters of the filter. Improving the theoretical understanding of these ad hoc procedures and developing an adaptive scheme are active research areas and interested readers can refer to review articles (e.g. Bannister, 2017; Carrassi et al., 2018; Vetra-Carvalho et al., 2018).

3 Experimental design

To evaluate the performance of the DA framework we perform a twin experiment. In this section we first describe the synthetic reference simulation that will be used to assess the performance of the DA framework. From this reference, we generate a set of synthetic noisy observations that will be used by the assimilation scheme. Finally, we describe the initial ensemble constructed using a priori or background information.

3.1 Reference simulation

We start by building an initial steady marine ice sheet. The domain extends from x=0 km, where we apply a symmetry condition, u=0 in Eq. (1), to x=800 km where we have a fixed calving front. We use 1D linear elements with a uniform mesh resolution of 200 m, leading to 4001 mesh nodes.

Following Durand et al. (2011), we generate a synthetic bed geometry that reproduces a typical large-scale overdeepening with some small-scale roughness. The bed $b=b_{\mathrm{trend}}+b_{\mathrm{r}}$ is the sum of a general trend b_trend defined as

$$\begin{array}{}\text{(20)}& b_{\mathrm{trend}}=\left\{\begin{array}{ll}-\mathrm{1100}+x& \text{for }\times \le \mathrm{450}\mathrm{km}\\ -\mathrm{650}-\mathrm{5}(x-\mathrm{450})& \text{for }\times >\mathrm{450}\mathrm{km},\end{array}\right.\end{array}$$

and a roughness signal b_r that is computed at 200 m resolution using a random midpoint displacement method (Fournier et al., 1982). This is a classical algorithm for artificial landscape generation. In 1D, the algorithm recursively subdivides a segment, and a random value drawn from a normal distribution 𝒩(0,σ²) is added to the elevation of the midpoint. The standard deviation σ is decreased by a factor of 2^h between two recursions. Here we have used 12 recursions using an initial standard deviation σ=500 m and a roughness h=0.7. The resulting bed is shown in Fig. 2.

Figure 2 (a) Reference ice sheet topography every 10 years from t=0 to t=200  a (black to grey). (b)–(d) GL position as a function of the simulation time for the reference (black line), for the ensemble (grey lines), and for the deterministic forecast (magenta line) (b) without assimilation, (c) with assimilation up to t=20  a and (d) with assimilation up to t=35  a. The right column shows a zoom on the first 40 years. In (c)–(d), the horizontal dashed line shows the end of the assimilation window.

For the basal friction, we use a synthetic sinusoidal function with two wavelengths for C ($\mathrm{MPa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$)

$$\begin{array}{}\text{(21)}& C=\mathrm{0.020}+\mathrm{0.015}\text{sin}\left(\mathrm{5}{\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{L}}\right)\text{sin}\left(\mathrm{100}{\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{L}}\right),\end{array}$$

with L=800 km (Fig. 3).

Figure 3 (a) RMSE between the reference and the analysed ensemble mean for the bed and friction coefficient. (b) For the bed and friction coefficient, the reference is shown in black, the synthetic bed measurements in the top panel are shown as green dots, the ensemble mean before assimilation is in blue and at t=35 a in red. The shading shows the ensemble spread between the minimum and maximum values, before assimilation (blue) and at t=35 a (red). The dashed vertical lines show the GL position at t=0 and t=35 a.

While not tuned to match any specific glacier, this synthetic design compares relatively well to the conditions found in Thwaites Glacier (Antarctica). Thwaites has been the focus of many recent studies as it is undergoing rapid ice loss and, connected to deep marine-based basins, its retreat could trigger a large-scale collapse of the West Antarctic Ice Sheet over the next centuries (Scambos et al., 2017). In Fig. 4, C and b are compared with model results from Brondex et al. (2019) along three streamlines. In Brondex et al. (2019), C has been inferred from the observed surface velocities using a time-independent control inverse method and a SSA model. We can see that our synthetic design is realistic in terms of both amplitude and spatial variations. As the other characteristics (geometry, small flow divergence and convergence) are also similar, the model velocities have a good order of magnitude.

Figure 4 Thwaites Glacier (Antarctica). Model results from Brondex et al. (2019): model velocities (top) and friction coefficient C and bed elevation b extracted along three streamlines (same colour code). Synthetic values used in this study are shown with black dashed lines. Note that the mesh resolution varies from ∽200 m close to the GL, shown in yellow in the top panel, to ∽10 km at the upstream end of the streamlines.

Using a uniform ice rigidity $B=(\mathrm{2}A{)}^{-\mathrm{1}/n}=\mathrm{0.4}\mathrm{MPa}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$, we grow an ice sheet to steady-state using a uniform surface accumulation $a_{\mathrm{s}}=\mathrm{0.5}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ and no basal melting a_b=0. The steady state GL is located at x=440 km, just downstream of the region of overdeepening (Fig. 2).

In Jenkins et al. (2018), observed ice flow accelerations in the Amundsen Sea sector have been attributed to the decadal oceanic variability, where warm phases associated with increased basal melt induce a thinning of the ice shelves, reducing their buttressing effect and initiating short-lived periods of unstable retreat of the most vulnerable GLs. In a flow line experiment the ice shelf does not exert any buttressing effect. Using a suite of melting and calving perturbation experiments for Pine Island Glacier, Favier et al. (2014) have shown that, when initiated, the dynamics of the unstable retreat are fairly independent of the type and magnitude of the perturbation. Here, to trigger the initial acceleration, we instantaneously decrease the ice rigidity to $B=\mathrm{0.3}\mathrm{MPa}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$ at t=0, keeping all the other parameters constant.

This initial perturbation induces an acceleration, a thinning and a retreat of the GL. The model is then run for 200 years with a time step d$t=\mathrm{5}{\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$. After a short stabilisation at x=437.2 km between t=13 a and t=32 a, the GL retreats at a rate of approximately 1 km a⁻¹ during the following 100 years, and then the rate decreases as the GL enters an area of down-slopping bed (Fig. 2). The retreat rate shows small variations associated with spatial variations in the topography and basal friction.

3.2 Synthetic observations

From the reference run, we generate synthetic noisy observations that are typical of the resolution and performance of actual observing systems.

For the bed, we mimic an airborne radar survey conducted perpendicular to the ice flow with an along-flow resolution of approximately 15 km. For this, we randomly select 54 locations between x=0 and x=800 km and then linearly interpolate the true bed and add a random uncorrelated Gaussian noise with a standard deviation ${\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{obs}}=\mathrm{20}\mathrm{m}$ (Fig. 3).

We assume that the surface elevation and velocities are observed at an annual resolution at each mesh node. We then add an uncorrelated Gaussian noise with a standard deviation ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}=\mathrm{10}\mathrm{m}$ for the surface elevation and ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}=\mathrm{20}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ for the velocity. The most recent velocity products are now posted with a monthly to annual resolution (Mouginot et al., 2017; Joughin et al., 2018). The reported uncertainty for individual velocity estimates using the 6 and 12 d image pairs from the Sentinel-1A/B satellites is 6.2 and 17.5 m a⁻¹ for the two horizontal velocity components in stable conditions; however this could be underestimated in the coastal areas. For the surface elevation, the spatial and temporal resolution as well as the coverage and uncertainty will depend on the sensors. The ArcticDEM (http://arcticdem.org, last access: 31 January 2019) is a collection of openly available digital surface models derived from satellite imagery and posted at 2 m spatial resolution. After co-registration, a standard deviation ranging from 2 to 4 m has been reported for the uncertainty of elevation difference between two individual models of static surfaces (Dai and Howat, 2017). Using the same satellites, Greenland digital elevation models are now posted with a 3-month temporal resolution (https://nsidc.org/data/nsidc-0715, last access: 28 August 2019).

3.3 Assimilation setup

We recall that our aim is to initialise the model using the DA framework to estimate the state together with the basal conditions. As a simplification to realistic experiments, we assume in the following that the ice rheological properties (represented by the Glen flow law and its parameters) and the forcing (represented by the surface and basal mass balances in Eq. 4) are perfectly known. In addition, we assume that the form of the basal friction follows Eq. (3) with $m=\mathrm{1}/\mathrm{3}$, so that only the spatially varying friction coefficient C is uncertain.

In our model, as the force balance Eq. (1) contains no time derivative, the velocity is a diagnostic variable. Because of the flotation condition, the topography can be represented by only one prognostic variable. The state vector x is then given by the free surface elevation z_s at every mesh node, and we use the floatation Eq. (5) for the mapping between the ice thickness H and z_s. The state vector is augmented by the two parameters to be estimated, the bedrock topography b and the basal friction coefficient C. For the parameters we assume a persistence model, i.e. no time evolution, during the forecast step (Eq. 8). Because the velocities are insensitive to the basal conditions where ice is floating, these two parameters are included in the state vector only for the nodes where at least one member is grounded. In addition, to insure that C remains positive, we use the following change of variable for the assimilation C=α². Although it does not insure uniqueness of the estimation as α and −α would lead to the same C, this change of variable is classical (MacAyeal, 1993) and was chosen as the reference friction coefficient spans only 1 order of magnitude. Similar performances were found using the other classical change of variable C=10^α as in Gillet-Chaulet et al. (2012).

Because both z_s and b are included in the state vector, the analysis does not conserve the ice sheet volume, for either the ensemble mean or the individual members. However, as illustrated in Fig. 1, the estimation of z_s and b, and thus of the ice thickness, should be improved at each analysis as more data are assimilated, and the final state is the best estimation provided by the filter knowing the model, all the observations during the assimilation window and their uncertainties. As mentioned in the introduction, if the main interest is an analysis of past volume changes, a smoother or a variational method might be more appropriate. The smoother extension of the ESTKF can be found in Nerger et al. (2014). Note however that, interestingly, if we expect that the filter will improve the estimation of the ice thickness, there is no guaranty in general that it will provide a better estimate of the total volume as an a priori state with a totally different thickness distribution could lead, by compensation of the errors, to a perfect estimate of the true volume.

Kalman-based filters are based on the hypothesis of the independence between the background, i.e. the initial ensemble, and the observations that are used during the assimilation. As the synthetic bed observations will be used to construct the initial ensemble (see next section), we assimilate only the surface elevation and velocity observations, every year from t=1 a up to t=35 a. The observation operator ℋ is a simple mapping for the surface elevation and is given by the non-linear SSA equation (Eq. 1) for the surface velocities.

Finally, to illustrate the effect of the transient assimilation on model projections on timescales relevant for sea level projections, the analysed states at t=20 a and t=35 a are used to run deterministic and ensemble forecasts up to t=200 a. The deterministic forecast uses the ensemble mean produced by the analysis while the ensemble forecast propagates the full ensemble.

3.4 Initial ensemble

For atmosphere and ocean models, the initial state is usually sampled from a climatology, either observed or from a model run. This method can not be used for the parameters and the initial ensemble must reflect the background and the estimation of its uncertainty, available a priori before the assimilation. Following previous studies (Gudmundsson and Raymond, 2008; Pralong and Gudmundsson, 2011; Bonan et al., 2014; Brinkerhoff et al., 2016), we assume that the initial distributions for b and C are Gaussian with a given mean and a prescribed covariance model. Furthermore we assume no cross-correlation between the initial b, C and z_s, and we draw the initial ensembles independently.

For b and C, the initial samples are drawn using the R package gstat (Pebesma and Wesseling, 1998). As is classical in geostatistics, the covariance model is prescribed using a variogram γ(d) that is half the variance of the difference between field values as a function of their separation d. It is usually defined by two parameters, the sill s that defines the semi-variance at large distances and the range r_a which, for asymptotic functions, is defined as the distance where the γ(r_a)=0.95s. The package gstat allows one to directly draw simulations, i.e. random realisations of the field, from the prescribed spatial moments (Pebesma and Wesseling, 1998).

For the bed we use an exponential function,

$$\begin{array}{}\text{(22)}& \mathit{\gamma }\left(d\right)=s(\mathrm{1}-e^{-\frac{\mathrm{3}d}{r_{\mathrm{a}}}}),\end{array}$$

with r_a=50 km and s=4000 m². We also add a nugget model defined by

$$\begin{array}{}\text{(23)}& \mathit{\gamma }\left(d\right)=\left\{\begin{array}{ll}\mathrm{0}& d=\mathrm{0}\\ \mathrm{nug}& d>\mathrm{0}\end{array}\right.,\end{array}$$

with nug=200 m². This model is meant to represent the bed measurement error. To draw the initial ensemble, the simulations are conditioned with the bed observations. This procedure gives an initial ensemble that is drawn from the posterior probability distribution that would be obtained using ordinary kriging with the same observations and variograms. The ensemble mean and spread for a 50-member ensemble are shown in Fig. 3 and the first three members are shown in Fig. 5. As expected, the ensemble spread increases with the distance from the observations. At the observation locations, the spread is controlled by the nugget. For the individual members, the nugget controls the small-scale variability, resulting in a roughness larger than the reference. When averaged this roughness disappears, and the ensemble mean has a much smoother topography.

Figure 5 For the initial ensemble for the bed and friction coefficient, the ensemble mean is the dashed blue curve, and the shading shows the ensemble spread. Coloured solid lines show the first three members. The reference is shown in black and the synthetic bed measurements are shown as green triangles.

For the friction coefficient, we assume that we know the mean value $C_{\mathrm{mean}}=\mathrm{0.020}\mathrm{MPa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$ and draw unconditional simulations. For the spatial dependence, we use a Gaussian function $\mathit{\gamma }\left(d\right)=s(\mathrm{1}-e^{-\mathrm{3}(\frac{d}{r_{\mathrm{a}}}{)}^{\mathrm{2}}})$ for the variogram using a range r_a=2.5 km and a sill $s={\mathrm{8.10}}^{-\mathrm{5}}{\mathrm{MPa}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}/\mathrm{3}}{\mathrm{a}}^{\mathrm{2}/\mathrm{3}}$. This results in initial ensemble members that have approximately the same maximal amplitude as the reference, as shown in Fig. 5.

For the free surface, we initialise all the members using the observed (noisy) free surface at t=0. Doing so, we implicitly assume that the spread of the ensemble induced by the uncertain initial conditions at the first analysis is small compared to the spread induced by the uncertain parameters. This is motivated by the fact that divergence anomalies induced by uncertainties in model parameters can typically reach tens to hundreds of metres per year in fast-flowing areas (Seroussi et al., 2011).

4 Results

4.1 Assimilation

To assess the performance of the DA in retrieving the basal conditions, we compute the root-mean-square error (RMSE) between the analysed ensemble mean and the reference for both the bed and the friction coefficient, RMSE_b and RMSE_C. After each analysis, the RMSE is computed using all the nodes where the basal conditions have been updated by the assimilation, i.e. at least one member is grounded, and where x≥ 300 km. The last value mentioned is close to the position reached by the grounding line after 200 years in the reference simulation; moreover, during the assimilation window, the reference velocity at this location is close to 80 m a⁻¹ (Fig. 6), so that the relative noise is ∼25 % and we do not expect too much improvement from the DA upstream as the velocity tends to 0.

Figure 6 Velocity, u, at t=1 a and t=35 a. The reference is in black; the ensemble mean before and after the analysis is in blue and red, respectively. The shading shows the ensemble spread between the minimum and maximum. The dashed vertical black line indicates grounding line position.

Here the size of the state vector x, N_x, is approximately 8400, i.e. z_s at every node and the basal conditions, b and C, in the grounded part. To test the performances of DA in conditions that would be numerically affordable for real applications, we run the assimilation with relatively small ensemble sizes N_e=30, N_e=50 and N_e=100. In this case, inflation and localisation are required to counteract the effects of undersampling and we test a range of forgetting factors ρ and localisation radius r. The errors obtained at t=20 a relative to the errors from the initial ensemble mean are shown in Fig. 7. The performances of the assimilation for N_e=50 and N_e=100 are very similar. The filter diverges and produces errors larger than the initial errors for a localisation radius r≤4 km. However, for larger localisation radii, the assimilation is relatively robust for a wide range of r and ρ, with errors reduced by ∼30 % for b and ∼40 % for C. Decreasing the ensemble sizes reduces the filter performance but there is still a reduction of the errors by ∼20 % and ∼30 %, respectively, with N_e=30. For the two smallest ensembles, there is an optimal value for r, and increasing r above this value decreases the filter performance. In general, this optimal value for r increases as ρ decreases, because the ensemble spread reduction induced by assimilating more observations is counterbalanced by the inflation.

Figure 7 RMSE at t=20 a, relative to the initial (before assimilation) RMSE as a function of the forgetting factor ρ and the localisation radius r for different ensemble sizes N_e. (a) For the bed and (b) for the friction coefficient. Black lines show isovalues spaced by 5 %.

In the sequel we discuss the results obtained with an ensemble size N_e=50. As a compromise between the performances in retrieving b and C, we choose a forgetting factor ρ=0.92 and a localisation radius r=8 km. The evolution of the RMSEs as a function of assimilation time together with the initial and final ensembles are shown in Fig. 3. RMSE_b decreases steadily from ∼25 m for the initial ensemble at t=0 to ∼12 m at t=35 a. For the basal friction, RMSE_C is decreased by a factor of 1.75 during the first 10 years; then there is still a slight but much smaller improvement as new observations are assimilated.

At the end of the assimilation, for both fields, the spatial variations are well reproduced by the ensemble mean, and, compared to the initial ensemble, the difference from the reference is decreased everywhere except between 300 and 325 km for C. The reduction in the error is also accompanied by a diminution of the ensemble spread, represented by the minimum and maximum values in Fig. 3. This reduction is the most important just upstream of the grounding line where the relative noise for the velocity is the smallest. For the first 100 km upstream of the grounding line, the ensemble standard deviation increases by a factor of 4, from approximately 4 to 17 m for b and from $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{MPa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$ for C. Downstream of the GL where all members are floating, the model is insensitive to the basal conditions and the initial ensemble is unchanged.

We expect that uncertainties in the ice sheet interior should not affect the short-term forecast of the coastal regions (Durand et al., 2011); however for completeness we also show the results for the first 300 km in Fig. 8. For the bed there is only a small improvement of the ensemble mean with an RMSE decreasing from 50 to 45 m after 35 years. Because the relative observation error on the velocity is very high in the first kilometres, the reduction of the ensemble spread due to the assimilation of new observations is very small and eventually outperformed by the inflation, leading to an ensemble spread that becomes larger than before the assimilation. The model seems more sensitive to the basal friction and this effect is less pronounced for C with a continuous decrease in the RMSE and a small reduction of the ensemble spread everywhere.

Figure 8 Same as Fig. 3 but with the RMSE computed for $x\in [\mathrm{0},\mathrm{300}]\mathrm{km}$.

Figure 2 shows that some members undergo a fast GL retreat of a few kilometres before assimilation at the end of the first year. Interestingly, as the assimilation updates both the thickness and the bed, it also corrects the GL position, which never departs by more than a few nodes from the reference for the rest of the assimilation period.

As in realistic simulations the true bed and friction are not available to assess the performance of the DA, we also look at the variables assimilated by the model. Figure 9 shows the RMSEs between the ensemble mean and the reference for the velocity u (RMSE_u) and the free surface z_s (${\text{RMSE}}_{z_{\mathrm{s}}}$), computed for the entire domain ($\mathrm{0}\le x\le \mathrm{800}\mathrm{km}$). We also report the evolution of the ensemble spread, computed as the square root of the averaged ensemble variance. The velocities before and after the analysis at t=1 and t=35 a are shown in Fig. 6. The RMSEs are largely decreased during the first few years, especially for the velocity with an error of more than 300 m a⁻¹ before the first assimilation to approximately the noise level 20 m a⁻¹, at t=20 a. For z_s, ${\text{RMSE}}_{z_{\mathrm{s}}}$ is already below the noise level before the first analysis and decreases relatively steadily to reach ∼2 m after 35 years. RMSE_u increases at the end of the period when the reference GL leaves the stable region. As can be shown in Fig. 6, the error is dominated by the larger difference over the ice shelf due to the few members that still have their GL at the stable location, largely affecting the ensemble mean.

Figure 9 RMSE (solid lines) and square root of the averaged ensemble variance (dashed lines) during the assimilation window for (a) the velocity, u, and (b) the free surface, z_s. Each year, the blue triangle and the red square are the RMSEs before and after the analysis, respectively. Each segment represents a 1-year forecast step.

In general, during the first and last years of the assimilation period, the error and the ensemble spread increase during the forecast step. The analysis step reduces both the error and the ensemble spread (Fig. 9). With the stabilisation of the grounding line, both the error and the spread remain relatively stable during the forecast, and as RMSE_u and ${\text{RMSE}}_{z_{\mathrm{s}}}$ have already reached levels comparable to the observation noise, there is no much improvement during the analysis. After a few assimilation steps, as expected for a reliable ensemble, the error and the spread have similar values.

Similar conclusions are drawn if the assimilation is pursued up to t=50 a. Because of the sensitivity of the ice shelf velocities to the grounding line position, RMSE_u shows a higher variability, but, with a few exceptions, stays close to the noise level. RMSE_b and RMSE_c stagnate as the reconstruction continually improves mostly in the first few tens of kilometres upstream of the GL where the relative noise on u is the smallest.

To assess the influence of the observation uncertainties in the performance of the DA, we repeat the experiment with the same localisation and inflation but different levels for the uncertainties on the observed surface velocity (${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}$) and surface elevation (${\mathit{\sigma }}_{z_{\mathrm{s}}}^{\mathrm{obs}}$) (see Sect. 3.2). We recall that these uncertainties are not correlated spatially and temporally. As shown in Fig. 10, the performance of the DA to retrieve both b and C increases when the uncertainty on the velocity observation ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}$ decreases. However, when looking at the model velocities and surface elevation, this improvement is not significant as the RMSEs were already below the noise level. As shown in Fig. 11, as expected, decreasing ${\mathit{\sigma }}_{z_{\mathrm{s}}}^{\mathrm{obs}}$ improves the analysis for the surface elevation. However, it does not necessarily reflect on the basal conditions and, on the contrary, reducing ${\mathit{\sigma }}_{z_{\mathrm{s}}}^{\mathrm{obs}}$ below 10 m leads to an increase in RMSE_C from 0.004 to $\mathrm{0.005}\mathrm{MPa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{a}}^{\mathrm{1}/\mathrm{3}}$. However, again, this effect does not reflect on the model velocities that are retrieved with the same accuracy.

Figure 10 Sensitivity to the surface velocity observation error ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}$ : RMSEs after each analysis, computed only for x≥300 km for b and C. The thick red lines correspond to the results with ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}=\mathrm{20}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ and ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}=\mathrm{10}\mathrm{m}$ shown in Figs. 3 and 9. The horizontal dashed lines correspond to the observation errors ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}$ and the results are presented with solid lines using the same colour code. ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}=\mathrm{10}\mathrm{m}$ for all the experiments.

Figure 11 Sensitivity to the surface elevation observation error ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}$: RMSEs after each analysis, computed only for x≥300 km for b and C. The thick red lines correspond to the results with ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}=\mathrm{20}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ and ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}=\mathrm{10}\mathrm{m}$ shown in Figs. 3 and 9. The horizontal dashed lines correspond to the observation errors ${\mathit{\sigma }}_{z_{\mathrm{s}}^{\mathrm{obs}}}$ and the results are presented with solid lines using the same colour code. ${\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{obs}}=\mathrm{20}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ for all the experiments.

4.2 Forecast simulations

We now discuss model projections from the initial state to t=200 a.

Without assimilation, the deterministic forecast, i.e. using the ensemble mean basal conditions, rapidly leads to the fastest GL retreat, and after a few years the GL position is no longer included within the previsions from the ensemble (Fig. 2b). This is due to the fact that the ensemble mean is smoother than the reference and any of the ensemble members. The reference GL position is included in the ensemble, and at the end of the simulations most of the members are within ±25 km from the reference. However, for a few members the GL remains very stable near its initial position for tens to hundreds of years, eventually never switching to an unstable regime during the duration of the simulation. Retreat rates are relatively variable from one member to the other, depending on the basal conditions.

With assimilation, the ensemble mean is improved and the difference from the reference reduced. The deterministic forecast cannot be distinguished from the ensemble members any more (Fig. 2c–d). Retreat rates are closer to the reference, with the previsions from all the ensemble members being more or less parallel to the reference. We note however that when the forecast starts after an assimilation window of 20 years, i.e. during a period of stable GL position for the reference, the deterministic forecast leaves the stable position with a delay of approximately 25 years, and a few members remain stable for the entire simulation. On average between t=13 and t=32 a, the thinning rate at the GL in the reference simulation is approximately 0.6 m a⁻¹, decreasing to 0.25 m a⁻¹ during the last 2 years. The total thinning between two analyses is then much lower than the noise in the observed surface elevation and cannot be captured accurately by the DA. In addition, at the GL, the difference between the minimum and maximum bed elevation given by the ensemble is approximately 20 m. This remaining uncertainty induces a difference of more than 2 m for the floatation surface and combined with the small thinning rates explains the delays in the initiation of the instability.

Extending the assimilation window up to t=35 a when the reference has switched in a fast retreat allows the forcing of all the members in the unstable retreat. There is a very good agreement between the reference and the deterministic forecast up to t=110 a. This is also true for the ensemble, and after that the spread is larger and the predicted GLs are less retreated than for the reference.

These results can be summarised by looking at the distribution of the ensemble forecasts for the grounding line position and volume above floatation (VAF) at t=100 a in Fig. 12 where the relative VAF change is computed as $({\mathrm{VAF}}_{t=\mathrm{100}}-{\mathrm{VAF}}_{t=\mathrm{0}}^{\mathrm{ref}})/{\mathrm{VAF}}_{t=\mathrm{0}}^{\mathrm{ref}}$, with ${\mathrm{VAF}}_{t=\mathrm{0}}^{\mathrm{ref}}$ the reference VAF at t=0. As expected there is a clear correlation between grounding line retreat and mass loss, with higher retreat leading to higher mass loss. The distributions are clearly non-Gaussian; however, even without assimilation there is already a mode close to the reference. The mode is more pronounced, with more members close to the reference as observations are assimilated. As discussed before, with no assimilation or a short assimilation up to t=20 a before the unstable retreat, the deterministic forecast can be very different from the mode of the ensemble forecast. However they are very similar if the assimilation is pursued up to t=35 a, within 1 % for the relative volume loss or 5 km for the GL position.

Figure 12 Ensemble forecast at t=100 a: (a)–(c) relative change of volume above floatation (VAF) and (bottom) GL position with (a, b) no assimilation, (b, e) assimilation up to t=20 a and (c, f) assimilation up to t=35 a. The red circle correspond to the reference run and the magenta square to the deterministic forecast.

5 Discussion

Here, we have tested an ensemble Kalman filter to assimilate annually observed surface velocities and surface elevation in a marine ice sheet model. Similar to previous studies, we have shown that, in fast-flowing regions, it is possible to accurately separate and recover both the basal topography and basal friction from surface observations (Gudmundsson and Raymond, 2008; Goldberg and Heimbach, 2013; Bonan et al., 2014; Mosbeux et al., 2016). In view of our results, because the synthetic bed observations were already used once to generate the initial ensemble, it seems unnecessary to assimilate these same observations again during each analysis as in Bonan et al. (2014).

Using a scheme that assimilates time-dependent observations provides a model state consistent with transient changes and that can directly serve as an optimal initial condition to run forecast simulations without the need of an additional relaxation (Goldberg and Heimbach, 2013; Goldberg et al., 2015). Interestingly the position of the grounding line is also corrected during the analysis step, and the ensemble quickly converges within a few grid nodes from the reference. In addition, the ensemble framework naturally allows the estimation and propagation of the uncertainty of the estimated parameters. Each assimilation of new data improves the reconstruction of the basal conditions (Fig. 3), and the first 10 years are very efficient in reducing error and the spread of the model surface velocities and elevation (Fig. 9). Furthermore, we have shown that the remaining uncertainties in the basal conditions do not significantly affect GL retreat rates once the unstable retreat is engaged. However, they can lead to considerable delays in the initiation of the instability. If the assimilation is pursued up to the beginning of the instability (35 years in our experiment) all the members exhibit the unstable retreat, and centennial-scale model projections converge to the reference (Fig. 12).

Good results have been obtained with relatively small ensembles (50 to 100 members) for a state vector of size N_x≈8400 and N_y=8002 observations. Similar to Bonan et al. (2014), we still see an improvement with a 30-member ensemble but the performances to retrieve the basal conditions are not as good. Running 2D plane view simulations with such ensemble sizes is largely possible as demonstrated by Ritz et al. (2015), who, using hybrid shallow ice–shallow shelf model, have run a 200-year ensemble forecast of the whole Antarctic Ice Sheet using 3000 members.

We have used inflation and localisation to stabilise the filter. The inflation giving the best results in Bonan et al. (2014) (ρ=0.87–1.02) is similar to the values tested in this study. For the localisation radius r we have used values between 4 and 16 km, while they range from 80 to 120 km in Bonan et al. (2014). While this seems counter-intuitive as the velocities depend only on the local conditions with the shallow ice approximation used by Bonan et al. (2014), in fact, because we use a different grid size (dx=0.2 km compared to dx=5 km in Bonan et al., 2014), for each node we assimilate twice as many observations. Our results are in agreement with the adaptive localisation radius proposed by Kirchgessner et al. (2014). Using three different models, Kirchgessner et al. (2014) have shown that good performances are obtained when r is such that the effective local observation dimension, defined as the sum of the weights attributed to each observation during the local assimilation, is equal to the ensemble size. Here the observation weights decrease with the distance to the local assimilation domain following a fifth-order polynomial function mimicking a Gaussian function (Gaspari and Cohn, 1999). The value r=8 km used for the 50-member ensemble gives an effective observation dimension of 56. Future studies should investigate if this result can be transposed to realistic 2D simulations with unstructured meshes.

In the experiments presented above, we have used a depth-integrated model for the force balance equations where GL migration is implemented through a hydrostatic floatation condition. This allows a full description of the ice topography with only one prognostic variable. Adaptation of the framework to a full-Stokes model requires minimum adaptations as these models do not rely on the floatation condition and solve a proper contact problem for the grounding line migration (Durand et al., 2009); this implies incorporating the two prognostic free surfaces z_b and z_s in the state vector. These models might be more sensitive to unbalanced geometries that could result from the analyses, especially when localisation is used (Cohn et al., 1998; Houtekamer and Mitchell, 2001). However, the ESTKF, as the ETKF, induces a minimal transformation of the ensemble members and thus has better chances to preserve balance (Nerger et al., 2012).

Before generalising such methods to real glacial systems, several points must be taken into consideration. They are independent of the DA method but they will eventually be treated differently in a variational or in an ensemble framework.

First, if the implementation is not an issue, the computational cost implied by running a full-Stokes model might remain a limiting factor. Compared to the Stokes solution, the SSA is known to overestimate the effects of bed topography perturbations on the surface profile for wavelengths less than a few ice thicknesses (Gudmundsson, 2008). How this issue can affect the reconstruction of the basal properties has never been quantified; however snapshot basal friction inversions have shown that the solution is sensitive to the force balance approximation (Morlighem et al., 2010). In addition, the MISMIP experiments have shown that the GL position and its response to a perturbation depend on the force balance solved by the models (Pattyn et al., 2012, 2013). In real applications, the performance of DA can be improved by explicitly taking into account the model error. Several strategies have been developed to account for this error, one approach with EnKFs being to use different versions of the model for different ensemble members (Houtekamer et al., 2009). Further studies could investigate the potential benefits of using ensembles that combine several force balance approximations.

Second, the quality of the analysis and the accuracy of the error estimates depends on the observation error covariance matrix R. It is then important to provide meaningful error estimates. Recent velocity maps provide an error estimate reported as the 1σ value for each individual location (Mouginot et al., 2012; Joughin et al., 2018). In general, this value agrees well with independent estimates; however care must be taken when the maps result from a composite of different sensors or different periods, and in general it might be difficult to properly estimate R.

A review paper by Tandeo et al. (2020) illustrates the impacts of badly calibrated observation and model error covariance matrices in a sequential DA framework and discusses available methods and challenges for their joint estimation. For the question of the impact of systematic errors, i.e. bias, either in the model or in the observations, and their correction by augmenting the system state in variational and ensemble DA, interested readers are referred to Dee (2005).

Third, the results depend on prior assumptions on the control variables and their variability, represented here by the initial ensemble. For the basal topography, current reference maps provide local error estimates (Fretwell et al., 2013; Morlighem et al., 2017); however they do not provide information about spatial correlations so that generating initial ensembles with the correct statistics might be problematic. In addition, the gridding can result in a loss of information for some regions of dense measurements, or it can lead to too smooth terrains in sparsely sampled areas. With the aim of generating terrains that have the correct high-resolution roughness, Graham et al. (2017) propose a synthetic 100 m resolution Antarctic bed elevation that combines the reference topography of Bedmap2 (Fretwell et al., 2013) with an unconditional simulation where the spatial correlation is fitted from dense radar measurements. This method could be used to generate initial ensembles but requires having access to the initial high-resolution measurements. Generating initial ensembles for the basal friction might be more problematic as there is in general no independent a priori information about the magnitude and spatial variability of the basal friction. If there is a correlation between the basal drag and the seismic observations of the bed conditions at a large scale, a proper physical theory is still missing to quantitatively incorporate such information in the models (Kyrke-Smith et al., 2017). It could be interesting to investigate how the existing multi-model basal friction reconstructions, based on snapshot inversions, could be used to derive initial uncertainty statistics and reduce the initial ensemble spread.

Finally, in our synthetic applications, we have not accounted for all potential sources of uncertainty which are, for example, as follows.

• The ice flow law. The ice viscosity depends on the englacial temperature, which itself is a function of the ice sheet history and the boundary forcing, including the geothermal heat flux (e.g. Van Liefferinge and Pattyn, 2013). Several other processes also affect the ice viscosity, including damage and strain-induced mechanical anisotropy (e.g. Pimienta et al., 1987; Schulson and Duval, 2009; Borstad et al., 2013). For the stress exponent, if the value n=3 is used by most models, published values ranges between 1 and 5 (e.g. Gillet-Chaulet et al., 2011).

• The friction law. More and more direct or indirect evidence shows that the friction under fast ice streams is at least partially controlled by the presence of sediments, leading to a Coulomb-type friction law (e.g. Tulaczyk et al., 2000; Murray, 1997; Joughin et al., 2010; Gillet-Chaulet et al., 2016). For hard beds, the development of subglacial cavities also implies deviations for the classical Weertman friction law (Schoof, 2005; Gagliardini et al., 2007).

• The density. The firn layer is not accounted for in most models; however its depth and density affect the floatation condition and thus the GL position (e.g. Griggs and Bamber, 2011). Directly assimilating the GL position, using, for example, the moving mesh approach developed by Bonan et al. (2017), would certainly be beneficial in realistic applications to reduce the discrepancy between the modelled and observed GL (Goldberg et al., 2015).

• The external forcings from the atmosphere and the ocean. Increasing mass loss rates from the ice sheets, in a large portion, can be attributed to a response to oceanic forcing, but multiple challenges remain for a proper assessment of their magnitude (Joughin et al., 2012).

Realistic simulations with ice flow models cover a wide range of spatial and temporal scales, and the relative importance of these uncertainties as well as their representation in the models will certainly have to be evaluated partly on a case-by-case basis, requiring the development of a robust framework for a variety of applications.

6 Conclusions

Developing model initialisation strategies that properly reproduce the ice sheet dynamical mass losses observed over the last decades requires developing transient assimilation frameworks that are able to account for the growing availability of dense time series, especially from space observations. Here, we presented a synthetic twin experiment demonstrating the possibility of calibrating a marine ice model using an ensemble Kalman filter which requires fewer numerical developments than variational methods.

Using resolutions and noise levels consistent with current observing systems, good performances are obtained to recover both the basal friction and basal topography with an ensemble of at least 50 members. Localisation and inflation have been tuned manually; however the results are consistent over relatively wide ranges. Future studies should investigate how these values can be transposed to realistic applications. Nevertheless, there is an abundant and growing literature in other geophysical fields to overcome problems that we might be facing in future studies.

Once the GL enters an unstable region, retreat rates largely depend on the basal conditions; thus using DA to reduce the associated uncertainties largely increases the skill of the model to predict rates and magnitude of GL retreat for timescales relevant for sea level rise projections. In our simplified application, the assimilation of the surface observations was sufficient to capture the GL migration during the assimilation window, without explicitly assimilating the observed position. However, for the GL to enter an irreversible retreat, the thickness must reach a tipping point, i.e. the thickness at the GL must reach floatation. This can seriously impact the predictability of the system as, for small perturbations, remaining uncertainties on the basal conditions can lead to an uncertainty on the residence time of the GL on stabilisation points, which can be similar to the simulation timescale. However, if the assimilation is pursued up to a time when the glacier is engaged in unstable retreat, all the members exhibit instability, and the spread of centennial-scale model projections, in terms of volume and grounding line position, is largely reduced.

Finally, we have discussed the main challenges to tackle before generalising transient DA in ice sheet modelling. This includes a better assessment of the uncertainties in the model and in the observations used for the background and for the assimilation.

Appendix A: Notations

Table A1 Notations and values used in this study associated with the ice flow model.

Table A2 Notations and values used in this study associated with the ensemble filter.