2.1 Ice sheet model ensemble

2.2 Observations

The calibration target is based on a compilation of five satellite altimeter datasets of surface elevation changes from 1992 to 2015 by Konrad et al. (2017). The synthesis involves fitting local empirical models over spatial and temporal extents of up to 10 km and 5 years, respectively, as developed by McMillan et al. (2014). The satellite missions show high agreement, with a median mismatch of 0.09 m a⁻¹. The dataset has a resolution of 10 km×10 km spatially and 6 months temporally. Only the last 7 years (beginning of 2008 to beginning of 2015) of the dataset are used here for calibration. The following satellite missions contributed to this period: ERS-2 (until 2011), Envisat (until 2012), ICESat (until 2009) and CryoSat-2 (2010 to 2015). All of these carry radar altimeters, with the only exception being ICESat, which had a laser altimeter (lidar) as payload.

There is no exact start date of the simulations, which makes a dating of the calibration period difficult. However, the ice flow observations from Rignot et al. (2011) used for the ice sheet initialisation are largely from a 3-year period centred around 2008, which is why this is the first year of surface elevation change observations we use. We do not correct for possible changes in firn thickness and directly convert surface elevation change rates of grounded ice into rates of ice thickness change. An average of all fourteen 6-month intervals is used for calibration; however for one calibration approach the averaging is performed in basis representation (see Sect. 3.2 for details).

3.1 Principal component decomposition

Let y(θ_i) be the m dimensional spatial model ice thickness change output for a parameter setting θ_i, where m is the number of horizontal grid cells, and the model ensemble has n members so that ${\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots }{\mathit{\theta }}_n=\mathbf{\Theta }$, $\mathbf{\Theta }\subset [\mathrm{0},\mathrm{1}{]}^d\subset R^d$ being the whole set of input parameters, spanning in our case the d=5 dimensional model input space. The m×n matrix $\stackrel{\mathrm{̃}}{\mathbf{Y}}$ is the row-centred combined model output of the whole Nias et al. (2016) ensemble with the ith column consisting of y(θ_i) minus the mean of all ensemble members, $\overline{\mathit{y}}$, and each row represents a single location. In the following we will assume n<m. A principal component decomposition is achieved by finding U, S and V so that

where the m×n rectangular diagonal matrix S contains the n positive singular values of $\stackrel{\mathrm{̃}}{\mathbf{Y}}$, and U and V^T are unitary. The rows of V^T are the orthonormal eigenvectors of ${\stackrel{\mathrm{̃}}{\mathbf{Y}}}^T\stackrel{\mathrm{̃}}{\mathbf{Y}}$, and the columns of U are the orthonormal eigenvectors of $\stackrel{\mathrm{̃}}{\mathbf{Y}}{\stackrel{\mathrm{̃}}{\mathbf{Y}}}^T$. In both cases the corresponding eigenvalues are given by diag(S)². By convention U, S and V^T are arranged so that the values of diag(S) are descending. We use B=US as shorthand for the new basis and call the ith column of B the ith principal component. The first five principal components have been normalised $\left(\frac{{\mathbf{B}}_i}{\left|{\mathbf{B}}_i\right|}\right)$ for Fig. 2 to show more detail of the spatial pattern.

The fraction of ensemble variance represented by a principal component is proportional to the corresponding eigenvalue of U, and typically there is a number k<n for which the first k principal components represent the whole ensemble sufficiently well. We choose k=5 so that 90 % of the model variance is captured (Fig. 2).

with B^′ and V^′ consisting of the first k columns of B and V.

This truncation limits the rank of $\stackrel{\mathrm{̃}}{\mathbf{Y}}$ to k=5. The PCs are by construction orthogonal to each other and can be treated as statistically independent.

Figure 2The first five normalised PCs of the model ice thickness change fields, building an orthogonal basis. They represent the main modes of variation in the model ensemble and are unitless since normalised. The lower left graph shows the fraction of total variance represented by each PC individually (grey) and cumulatively (red), based on squared singular values.

3.2 Observations in basis representation

One of the calibration approaches we investigate uses the PCs derived before for both the model and observations (see Sect. 3.4). For this we have to put the observations onto the same basis vectors (PCs) as the model data. Spatial m dimensional observations z_(xy) can be transformed to the basis representation by

for z_(xy) on the same spatial grid as the model output y(θ), which has the mean model output $\overline{\mathit{y}}$ subtracted for consistency.

We perform the transformation as in Eq. (3) for all of the biyearly observations over a 7-year period to get 14 different realisations of $\widehat{\mathit{z}}$. Due to the smooth temporal behaviour of the ice sheet on these timescales we use the observations as repeated observations of the same point in time to specify $\widehat{\mathit{z}}$ as the mean and use the variance among the 14 realisations of $\widehat{\mathit{z}}$ to define the observational uncertainty in the calibration model (Sect. 3.4).

Figure 3(a) Mean observed ice thickness change from 2008 to 2015 based on data from Konrad et al. (2017). (b) As (a) but projected to first five PCs and reprojected to spatial field.

Figure 3 shows that large parts of the observations can be represented by the first five PCs from Fig. 2. This is supported by the fact that the spatial variance (Var()) of the difference between the reprojected and original fields is substantially smaller than from z_(xy) alone:

It is only the part of the observations which can be represented by five PCs (right of Fig. 3) which will influence the calibration.

3.3 Emulation

For a probabilistic assessment we need to consider the probability density in the full, five-dimensional parameter space. This exploration can require very dense sampling of probabilities in the input space to ensure appropriate representation of all probable parameter combinations. This is especially the case if the calibration is favouring only small subsets of the original input space. In our case more than 90 % of the calibrated distribution would be based on just five BISICLES ensemble members. For computationally expensive models sufficient sampling can be achieved by statistical emulation, as laid out in the following.

A row of ${{\mathbf{V}}^{\prime }}^T$ can be understood as indices of how much of a particular principal component is present in every ice sheet model simulation. Emulation is done by replacing the discrete number of ice sheet model simulations by continuous functions or statistical models. We use each row of ${{\mathbf{V}}^{\prime }}^T$, combined with Θ, to train an independent statistical model where the mean of the random distribution at θ is denoted ω_i(θ). Here the training points are noise-free as the emulator is representing a deterministic ice sheet model, and therefore ${\mathit{\omega }}_i\left(\mathbf{\Theta }\right)=[{{\mathbf{V}}^{\prime }}^T{]}_i$ for principal components $i=\mathrm{1},\mathrm{\dots },k$. Each of those models can be used to interpolate (extrapolation should be avoided) between members of Θ to predict the ice sheet model behaviour and create surrogate ensemble members.

We use Gaussian process (GP) models, which are a common choice for their high level of flexibility and inherent emulation uncertainty representation (Kennedy and O'Hagan, 2001; O'Hagan, 2006; Higdon et al., 2008). The random distribution of a Gaussian process model with noise-free training data at a new set of input values θ_* is found by (e.g. Rasmussen and Williams, 2006)

where $N(\cdot ,\cdot )$ represents a multivariate normal distribution and the values of $K(\mathbf{\Theta },\mathbf{\Theta }{)}_{ij}=c({\mathit{\theta }}_i,{\mathit{\theta }}_{\mathit{j}})$ are derived from evaluations of the GP covariance function $c(\cdot ,\cdot )$. Equivalent definitions are used for $K({\mathit{\theta }}_{*},\mathbf{\Theta })$, $K(\mathbf{\Theta },{\mathit{\theta }}_{*})$ and $K({\mathit{\theta }}_{*},{\mathit{\theta }}_{*})$; note that $K({\mathit{\theta }}_{*},{\mathit{\theta }}_{*})$ is a 1×1 matrix if we emulate one new input set at a time. We use a Matern ($\frac{\mathrm{5}}{\mathrm{2}}$) type function for $c(\cdot ,\cdot )$ which describes the covariance based on the distance between input parameters. Coefficients for $c(\cdot ,\cdot )$ (also called hyper-parameters), including the correlation length scale, are optimised on the marginal likelihood of ω(Θ) given the GP. We refer to Rasmussen and Williams (2006) for an in-depth discussion and tutorial of Gaussian process emulators.

Due to the statistical independence of the principal components we can combine the k GPs to

The combined Ω is in the following called emulator and ω(θ), and the entries of the diagonal matrix Σ_ω(θ) follow from Eq. (4). We use the python module GPy for training (GPRegression()) and marginal likelihood optimisation (optimize_restarts()). In total we generate more than 119 000 emulated ensemble members. Emulator estimates of ice sheet model values in a leave-one-out cross-validation scheme are very precise with squared correlation coefficients for both emulators of R²>0.988 (see Supplement for more information).

3.4 Calibration model

Given the emulator in basis representation, a calibration can be performed either after reprojecting the emulator output back to the original spatial field (e.g. Chang et al., 2016a; Salter et al., 2018) or in the basis representation itself (e.g. Higdon et al., 2008; Chang et al., 2014). Here we will focus on the PC basis representation.

We assume the existence of a parameter configuration θ^* within the bounds of Θ (the investigated input space) which leads to an optimal model representation of the real world. To infer the probability of any θ to be θ^* we rely on the existence of observables, i.e. model quantities z for which corresponding measurements $\widehat{\mathit{z}}$ are available. We follow Bayes' theorem to update prior (uninformed) expectations about the optimal parameter configuration with the observations to find posterior (updated) estimates. The posterior probability of θ being the optimal θ^* given the observations is

where $L(\mathit{z}=\widehat{\mathit{z}}|\mathit{\theta })$ is the likelihood of the observables to be as they have been observed under the condition that θ is θ^*, and π(θ) is the prior (uninformed) probability that $\mathit{\theta }={\mathit{\theta }}^{*}$. Following Nias et al. (2016) we choose uniform prior distributions in the scaled parameter range [0,1] (see also Sect. 2 and Eq. 11 in Nias et al., 2016). The emulator output is related to the real state of the ice sheets in basis representation, γ, by the model discrepancy ε:

We assume the model discrepancy to be multivariate Gaussian distributed with zero mean; $\mathit{\epsilon }=N(\mathrm{0},{\mathbf{\Sigma }}_{\mathit{\epsilon }})$. The observables are in turn related to γ by

where e is the spatial observational error and the transformation $({{\mathbf{B}}^{\prime }}^T{\mathbf{B}}^{\prime }{)}^{-\mathrm{1}}{{\mathbf{B}}^{\prime }}^T$ follows from Eq. (3).

We simplify the probabilistic inference by assuming the model error/discrepancy ε, the model parameter values Θ, and observational error e to be mutually statistically independent and e to be spatially identically distributed with variance ${\mathit{\sigma }}_e^{\mathrm{2}}$, so that

The k×k matrix $({{\mathbf{B}}^{\prime }}^T{\mathbf{B}}^{\prime }{)}^{-\mathrm{1}}$ is diagonal with the element-wise inverse of diag$({\mathbf{S}}^{\prime }{)}_i^{\mathrm{2}}$ as diagonal values. We estimate ${\mathit{\sigma }}_e^{\mathrm{2}}$ from the variance among the 14 observational periods for the first principal component constituting $\widehat{z_{\mathrm{1}}}$; i.e.

Note that the existence of γ is an abstract concept, implying that it is only because of an error ε that we cannot create a numerical model which is equivalent to reality. However abstract, it is a useful, hence common, statistical concept allowing us to structure expectations of model and observational limitations (Kennedy and O'Hagan, 2001). Neglecting model discrepancy, whether explicitly by setting ε=0, or implicitly, would imply that an ice sheet model can make exact predictions of the future once the right parameter values are found. This expectation is hard to justify considering the assumptions which are made for the development of ice sheet models, including sub-resolution processes. Neglecting model discrepancy typically results in overconfidence and potentially biased results.

The inclusion of model discrepancy can at the same time lead to identifiability issues where the model signal cannot be distinguished from the imposed systematic model error. Constraints on the spatial shape of the discrepancy have been used to overcome such issues (Kennedy and O'Hagan, 2001; Higdon et al., 2008). An inherent problem with representing discrepancy is that its amplitude and spatial shape are in general unknown. If the discrepancy were well understood the model itself or its output could be easily corrected. Even if experts can specify regions or patterns which are likely to show inconsistent behaviour, it cannot be assumed that these regions or patterns are the only possible forms of discrepancy. If its representation is too flexible, it can however become numerically impossible in the calibration step to differentiate between discrepancy and model behaviour.

For these reasons we choose a rather heuristic method which considers the impact of discrepancy on the calibration directly and independently for each PC. Therefore Σ_ε is diagonal with diag$\left({\mathbf{\Sigma }}_{\mathit{\epsilon }}\right)=({\mathit{\sigma }}_{\mathit{\epsilon }\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{\sigma }}_{\mathit{\epsilon }k}^{\mathrm{2}}{)}^T$. The three-sigma rule states that at least 95 % of continuous unimodal density functions with finite variance lie within three standard deviations from the mean (Pukelsheim, 1994). For the ith PC we therefore find ${\mathit{\sigma }}_{i\mathrm{95}}^{\mathrm{2}}$ so that 95 % of the observational distribution $N(\widehat{z_i},{\mathit{\sigma }}_{ei}^{\mathrm{2}})$ lies within 3σ_i95 from the mean of ω_i(Θ), i.e. across the n ensemble members. We further note that we do not know the optimal model setup better than we know the real state of the ice sheet and set the minimum discrepancy to the observational uncertainty. Hence ${\mathit{\sigma }}_{\mathit{\epsilon }i}^{\mathrm{2}}=max({\mathit{\sigma }}_{i\mathrm{95}}^{\mathrm{2}},{\mathit{\sigma }}_{ei}^{\mathrm{2}})$.

We thereby force the observations to fulfil the three-sigma rule by considering them as part of the model distribution ω_i(Θ) while avoiding overconfidence in cases where observations and model runs coincide.

3.5 Calibration model test

In this section we test our calibration approach on synthetic observations to see whether our method is capable of finding known-correct parameter values. We select one member of the BISICLES model ensemble at a time and add 14 different realisations of noise to it. The noise is added to see how the calibration performs if the observations cannot be fully represented by the ice sheet model.

We use spatially independent, zero-mean, normally distributed, random noise with variance equal to the local variance from the 14 periods of satellite observations. This way the variance incorporates dynamic changes (acceleration/deceleration of the ice thickness change) and technical errors (e.g. measurement and sampling errors). For each selected model run we generate 14 noise fields and add them to the single model ice thickness change field. These 14 realisations replace the 14 periods of satellite observations for the synthetic model tests.

For Fig. 4 the model run with central parameter values (=0.5) for basal traction, viscosity and ocean melt scaling factors, non-linear friction, and modified bedrock has been selected, as indicated by black circles. This parameter set has been selected as it highlights the limitations of the calibration; the results of 11 other synthetic model tests are shown in the Supplement.

Figure 4Likelihood of parameter combinations of the synthetic test case (evaluations of Eq. 12). Upper right panels show likelihood values marginalised to pairs of parameters including basal melt factor (B. Melt) and friction law exponent (Fr. Law), normalised to the respective maximum for clarity. Lower left panel shows likelihood values marginalised to individual parameters for the three scalar parameters (line plots), as well as the friction law exponent and bedrock topography map (labelled as L(Lin. Sliding) and L(Mod. Bed) respectively), normalised to an integral of one, consistent with probability density functions. The central values for traction, viscosity, and ocean melt as well as non-linear friction and modified bedrock are used. The parameter values are also shown by the black circles, while the values of the set of parameters with highest likelihood are shown by green crosses.

Download

Figure 4 illustrates which parts of the model input space are most successful in reproducing the synthetic observations of ice thickness changes during the calibration period. For visualisation we collapse the five-dimensional space onto each combination of two parameters and show how they interact. For a likely (yellow) area in Fig. 4 it is not possible to see directly what values the other three parameters have, but very unlikely (black) areas indicate that no combination of the remaining parameter values results in model configurations consistent with observations.

As can be seen from Fig. 4, marginal likelihoods of our calibration approach can favour linear friction even if the synthetic observations use non-linear friction. In addition, the ocean melt parameter is often weakly constrained or, as in this case, biased towards small melt factors. In contrast, the basal traction coefficient and viscosity scaling factors have a strong mode at, or close to, the correct value of 0.5, and the correct bedrock map can always be identified (Fig. 4 and Supplement). Different values of basal traction and viscosity have been tested in combination with both bedrock maps and show similar performance (see Supplement). The fact that the parameter setup used for the test is attributed to the maximal likelihood (green cross on top of black circle) supports our confidence in the implementation as the real parameter set is identified correctly as the best fit. Relative ambiguity with respect to friction law and ocean melt overrules the weak constraints on these parameters in the marginalised likelihoods. The higher total likelihood of linear friction can be traced back to a higher density of central ensemble members for linear friction. Non-linear friction produces more extreme ice sheet simulations as simulations with high velocities will have reduced (compared to linear friction) basal drag and speed up even more (and vice versa for simulations with slow ice flows). The frequency distribution of total sea level contribution and basis representation are therefore wider for non-linear friction (see Supplement). The relative density of ensemble members around the mode of the frequency distribution can, as for this test case, cause a smaller marginal likelihood for non-linear friction compared to linear friction (28 % to 72 %). This can be considered a caveat of the model ensemble which might very well be present in other ensembles which perturb the friction law in combination with other parameters. If the friction law cannot be adequately constrained, as is the case for all calibration approaches tested here, the prior belief in the optimal friction law must be set very carefully.

The signal of friction law and ocean melt is not strong enough to adequately constrain the calibration, even though both parameters are known to have a strong impact on the ice sheet (Pritchard et al., 2012; Arthern and Williams, 2017; Jenkins et al., 2018; Joughin et al., 2019; Brondex et al., 2019). This is likely related to the slower impact of those parameters compared to the others. A change in bedrock, basal traction or viscosity has a much more immediate effect on the ice dynamics. For example, if the basal traction field is halved, the basal drag will be reduced by the same amount, leading to a speed-up of the ice at the next time step (via the solution of the stress balance). The perturbation of ocean melt from the start of the model period has to significantly change the ice shelf thickness before the ice dynamics upstream are affected. The initialisation of the ensemble has been performed for each friction law individually, which means that the initial speed of the ice is by design equivalent. It is only after the ice velocities change that the different degrees of linearity in the friction law have any impact on the simulations. This does not mean that the simulations are insensitive to the ocean melt forcing and friction law; in fact Fig. 4 shows that both parameters have some impact on the simulation in the calibration period. It just means that the much more immediate effects of basal traction and viscosity are likely to dominate the calibration on short timescales.

From this test we conclude that basal friction law and ocean melt scaling cannot be inferred with this calibration approach and calibration period. We will therefore only calibrate the bedrock as well as basal traction and viscosity scaling factors. Several studies used the observed dynamical changes of parts of the ASE to test different friction laws. Gillet-Chaulet et al. (2016) find a better fit to evolving changes of Pine Island Glacier surface velocities for smaller m, reaching a minimum of the cost function from around $m=\mathrm{1}/\mathrm{5}$ and smaller. This is supported by Joughin et al. (2019), who find $m=\mathrm{1}/\mathrm{8}$ to capture the PIG speed-up from 2002 to 2017 very well, matched only by a regularised Coulomb (Schoof) friction law. It further is understood that parts of the ASE bed consist of sediment-free, bare rocks for which a linear Weertman friction law is not appropriate (Joughin et al., 2009). We therefore select non-linear friction by expert judgement and use a uniform prior for the ocean melt scaling.

3.6 Comparison with other calibration approaches

To put the likelihood distribution from Fig. 4 into context, we try two other methodical choices. First we calibrate in the spatial domain after reprojecting from the emulator results.

where y^′(θ_i) are the reprojected ice sheet model results after truncation for parameter setup θ. We set the model discrepancy to twice the observational uncertainty ${\mathit{\sigma }}_e^{\mathrm{2}}$ so that the reprojected likelihood L_(xy) simplifies to

Another approach is to use the net yearly sea level contribution (SLC) from the observations SLC(z_(xy)) and model SLC(y^′(θ_i)) for calibration, as done in e.g. Ritz et al. (2015).

Again, we set the model discrepancy to twice the observational uncertainty which we find from the variance of the yearly sea level contributions for the 14 biyearly satellite intervals. ${\mathit{\sigma }}_{\mathrm{SLC}}^{\mathrm{2}}=\mathrm{Var}\left(\mathrm{SLC}\right({\mathit{z}}_{\left(xy\right)}\left)\right)={\mathrm{0.035}}^{\mathrm{2}}$ [mm SLE² a⁻²].

Figure 5Likelihood of parameter combinations of the synthetic test case for reprojected emulator estimates (top, a; Eq. 14) and sea level rise contribution calibration (bottom, b; Eq. 15). Upper right panels show likelihood values marginalised to pairs of parameters including basal melt factor (B. Melt) and friction law exponent (Fr. Law), normalised to the respective maximum for clarity. Lower left panel shows likelihood values marginalised to individual parameters for the three scalar parameters (line plots), as well as the friction law exponent and bedrock topography map (labelled as L(Lin. Sliding) and L(Mod. Bed) respectively), normalised to an integral of one, consistent with probability density functions. The central values for traction, viscosity, and ocean melt as well as non-linear friction and modified bedrock are used. The parameter values are also shown by the black circles, while the values of the set of parameters with highest likelihood are shown by green crosses.

Download