Assessing the Impact of Retreat Mechanisms in a Simple Antarctic Ice Sheet Model Using Bayesian Calibration

The response of the Antarctic ice sheet (AIS) to changing climate forcings is an important driver of sea-level changes. Anthropogenic climate change may drive a sizeable AIS tipping point response with subsequent increases in coastal flooding risks. Many studies analyzing flood risks use simple models to project the future responses of AIS and its sea-level contributions. These analyses have provided important new insights, but they are often silent on the effects of potentially important processes such as Marine Ice Sheet Instability (MISI) or Marine Ice Cliff Instability (MICI). These approximations can be well justified and result in more parsimonious and transparent model structures. This raises the question of how this approximation impacts hindcasts and projections. Here, we calibrate a previously published and relatively simple AIS model, which neglects the effects of MICI and regional characteristics, using a combination of observational constraints and a Bayesian inversion method. Specifically, we approximate the effects of missing MICI by comparing our results to those from expert assessments with more realistic models and quantify the bias during the last interglacial when MICI may have been triggered. Our results suggest that the model can approximate the process of MISI and reproduce the projected median melt from some previous expert assessments in the year 2100. Yet, our mean hindcast is roughly 3/4 of the observed data during the last interglacial period and our mean projection is roughly 1/6 and 1/10 of the mean from a model accounting for MICI in the year 2100. These results suggest that missing MICI and/or regional characteristics can lead to a low-bias during warming period AIS melting and hence a potential low-bias in projected sea levels and flood risks.


Introduction
Coastal areas are at risk to sea-level rise, and will be more so if the marine part of the Antarctic ice sheet (AIS) were to collapse. A disintegration of the marine part of the AIS would raise global mean sea level by more than three meters and hence drive an increase in flooding vulnerability in many coastal areas [1]. Basic physics, the paleo-record, as well as model simulations suggest that such a disintegration could follow a highly nonlinear, relatively abrupt, and hysteresis type threshold response [2][3][4][5]. Previous studies suggest that anthropogenic greenhouse gas emissions could trigger a disintegration of the marine part of the AIS [2,[6][7][8], and implies that assessing these potential increased risks is important.
Hindcasts and projections of AIS dynamics are deeply uncertain. This uncertainty stems in part from the complexity of the coupled AIS / ocean / atmosphere systems combined with relatively sparse data and the severely limited ability to resolve relevant processes in current models [3,5,9,10]. For example, consider the highly vulnerable marine part of the AIS. The marine part of the AIS is vulnerable to anthropogenic climate change because it is grounded below the sea level and in direct contact with the ocean via ice shelves [2,8,11]. Ocean-ice interactions can in part drive AIS dynamics, for example, through melting of ice shelves and calving cliffs [5,12]. For the part of the ice sheet grounded below sea level, this melting and calving can lead to a runaway retreat due to threshold behavior, known as Marine Ice Sheet Instability (MISI) [2,3,5]. MISI is a runaway retreat of the grounding line as warming temperatures increase both the water depth and ice flux at the grounding line [8]. Once started, this mechanism will continue (even in the absence of external forcings) until the grounding line reaches an upward-sloping bed or enough buttressing is exerted to stabilize the grounding line on a retrograde slope [13]. Previous studies suggest that the additional process of Marine Ice Cliff Instability (MICI) may facilitate MISIs and therefore ice sheet disintegration [2,3,5,8,12]. MICI is defined as a weakening or structural failure of shear ice cliffs as warming temperatures increase crevasses and reduce the maximum supported cliff heights [8].
Here, we calibrate a previously published and relatively simple model of AIS volume loss [14] to reproduce a hindcast period from the last interglacial to the present using a pre-calibration method and a Bayesian inversion. The Bayesian inversion accounts for the heteroskedastic nature of the data and the model accounts for MISI (Fig 1). However, the model neglects the effects of MICI and is unable to capture regional characteristics of the ice sheet. To approximate the effects of missing MICI, we compare our projections to physically more realistic models [8,[15][16][17][18] and quantify bias during the last interglacial period, when MICI was potentially triggered. In the results section, we show that our projections are comparable to some previous expert assessments [4,[15][16][17][18], yet unable to account for the sizable melt generated from a model considering MICI [8] and that our mean hindcast is unable to reproduce the estimated mean during the last interglacial. In our analysis, we deduce that missing MICI is one possible mechanism that can lead to underestimation of warming period (i.e., the last interglacial and the year 2100) melting in the AIS model.

Antarctic ice sheet model
We adopt the DAIS (Danish Center for Earth System Science Antarctic Ice Sheet) model as a starting point for our analysis [14]. This model has been previously described in detail (see for example, [14]), hence we provide just a brief overview. The DAIS model considers three important changes in the ice sheet: ice sheet volume, volume loss in sea-level equivalence (SLE), and ice sheet radius. In this model, the volume for the entire AIS changes under four mechanisms: (i) ice sheet growth from precipitation, (ii) accumulation from precipitation minus melt runoff, (iii) accumulation from precipitation minus basal melt, and (iv) ice sheet disintegration from precipitation minus melt runoff and basal melt. The simplicity of the model stems from neglecting key processes and threshold responses. For example, the model misses the link of warming temperatures to hydro-fracturing of buttressing ice shelves, collapse of ice cliffs, and abrupt changes in melting rates by surpassing tipping points. Additionally, the simplicity of the model treats the AIS as a whole rather than capturing the regional characteristics of glaciers and subglacial basins, which would help trigger the MISI (see the section "Potential explanation for the discrepancy" for more details). However, the model quantifies the ice flux (the rate at which the ice flows towards the sea) at the grounding line (location where the ice shelf begins to float; averaged for the entire AIS), the water depth at the grounding line (distance between the sea level and the grounding line; averaged for the entire AIS), and the relationship between ice sheet volume and ice sheet radius [14] (Fig 1). In the model, ice sheet volume relates to the radius based on an undisturbed bed profile referenced to present day volume and radius taking into account isostatic adjustment and its effect on the displacement of seawater by ice [11,14]. Despite the drawbacks, the simplicity and computational efficiency of the DAIS model enables us to implement a Bayesian data-model fusion and to quantify uncertainties with a focus on the tails of the distributions.

Data
Forcings. The ice sheet forcings in the hindcast period span from 240 kyr BP to the year 1997 [14]. These forcings include Antarctic temperature reduced to sea-level, global mean sea level, global mean sea-level rate, and high latitude subsurface ocean temperature.
For projections beyond 1997, we generate forcings with the CNRM CM5 RCP8.5 temperature scenario (as obtained from the CMIP5 model output archive, http://cmip-pcmdi.llnl.gov/ cmip5/). This model simulation was chosen because both surface air and ocean temperatures were projected to the year 2300 using the extended RCP8.5 scenario. The future forcings are created using details and instructions specified in a previous analysis [14]. The Antarctic temperature anomaly is adjusted to a present-day (1961-1990) mean temperature of -18˚C and projected over a 1˚lat./long. Antarctica land mask [14]. We generate future global mean sea level and sea-level rates with an empirical global mean sea-level model [19] calibrated using historical NASA GISS global mean temperatures [20,21] and global mean sea-level observations [14] by minimizing the sum of the absolute residuals using a global optimization method [22]. The high latitude subsurface ocean temperature is the volume averaged temperature from 200-800 m depth between 52-70˚S [14]. In addition to the business-as-usual temperature scenario, the ocean temperatures are projected with the Antarctic temperature forcing and adjusted to a present-day mean temperature of 0.72˚C.
Observational constraints. We use four observational constraints of volume loss in SLE to fit the model: (i) the last interglacial (LIG; *120 kyr BP), (ii) the last glacial maximum (LGM; *20 kyr BP), (iii) the mid-Holocene (MH; *6 kyr BP), and the instrumental period (1992-2011) (S1 Table). Each observational constraint is based on previously published work and interpreted as a 95% confidence interval. These ranges include: 1.8-6.0 m (LIG), -6.9--15.8 m (LGM), -1.25--4.0 m (MH), and 1-2.9 mm (in the year 2002). Note that the instrumental period constraint is relative to the year 1992, whereas the other constraints are relative to the 1961-1990 period. A detailed discussion of the constraints is provided in the Supporting Information (S1 Text and S1 Table).

Calibration methods
We use two separate calibration stages: (i) a pre-calibration (e.g., [23]) using Latin hypercube sampling (LHS) and (ii) a full Bayesian inversion based on a Markov chain Monte Carlo (MCMC) algorithm. We compare these two methods to assess the potential improvement in the hindcasts and parameter estimates associated with moving from the simpler pre-calibration to a full Bayesian inversion and to provide a check for the Bayesian inversion.
Pre-calibration. We first implement pre-calibration using a LHS technique [24,25]. We generate 1.3 × 10 3 samples from the 13-dimensional model parameter space. We use uniform prior distributions for each physical model parameter and an inverse gamma distribution for the statistical parameter σ 2 P [26,27]. The parameters σ 2 P and σ 2 I approximate the effects of observational errors, unresolved internal variability, and model errors for the paleo and instrumental period (discussed below). LHS divides each parameter distribution into 1.3 × 10 3 equally probable intervals and then draws 1.3 × 10 3 sample points. We employ a maxi-min LHS, which optimizes the sample by maximizing the minimum distance between parameter values [25,28]. Using this method, we generate 1.3 × 10 3 parameter combinations and determine which parameter sets satisfy individual observational constraints (Table 1; hindcasts and projections shown in S1 and S2 Figs).
Markov chain Monte Carlo (MCMC). We implement a Bayesian inversion technique using a MCMC algorithm [29][30][31]. We represent the observational data of AIS SLE by y t , the model output as f(θ, t), the unknown model parameters θ, an unknown residual term R t , and time t over a hindcast period of 240 kyr BP to 1997. We approximate the AIS volume loss observations as the sum of the model output plus an unknown residual term, The observations, model output, and residuals are vectors from 1 to N (N = 4 in this case, since there are four observational constraints) such that y = (y 1 , . . ., y N ) t , f(θ, t) = (f(θ, t 1 ), . . ., f(θ, t N )) t , and R = (δ 1 + 1 , . . ., δ N + N ) t . Because thousands of years separate the observational constraints, we neglect the temporal correlation. Hence, we model the residuals R t as drawn from an independent and identically distributed normal distribution δ t with an added observation error t , This is a simple model for the data model discrepancy since accounting for independent and identically distributed residuals ignores the potentially complex autocorrelation structure [32] of AIS trends. The autocorrelation structure can be modeled from year to year if more data are available [32]. In our case, our simple model seems appropriate, because there are four data points with large time spans between observations. We hence use models of white noise, δ t that are independent normally distributed with different variance for the paleo and instrumental record, for instance, d t $ Nð0; s 2 P Þ for t = 1, 2, 3 and d t $ Nð0; s 2 I Þ for t = 4. The observation error ε t is heteroskedastic consisting of a normal distribution with zero means and a time dependent variance, t $ Nð0; t 2 t Þ. The white noise δ = (δ 1 , . . ., δ N ) t approximates the model structural uncertainty and internal climate variability whereas the time dependent variance of ϵ = ( 1 , . . ., N ) t is implemented by substituting in the known measurement error of each observational constraint. This implementation is important because ignoring the potentially complex error structure can result in overconfident parameter estimates, hindcasts, and projections (e.g., [32]). The data assimilation technique follows Bayes' theorem [33], where the posterior probability π of observing the parameters Θ given the data y is proportional to the likelihood L of the data given the parameters times the prior probability distribution of the parameters, pðY j yÞ / Lðy j YÞ Â pðYÞ: After specifying the prior parameter distributions (S3 Fig) and the likelihood function, the MCMC algorithm samples from the posterior distribution [29][30][31]34]. We construct the likelihood function, assuming that the observational data y consists of the model output f(θ, t) plus an unknown residual term R (Eqs 1 and 2). The uncertain parameters Θ include both the unknown model parameters θ and the unknown statistical parameters σ 2 P and σ 2 I (S3 and S4  Figs). To derive the likelihood of the data given the parameters L(yjΘ), we find the residuals R * N(0, S); The likelihood functions is then: We adopt uniform priors for the model parameters and an inverse gamma prior for the statistical parameter σ 2 P (α = 2; β = 1) (S3 Fig). The inverse gamma prior for the paleo period σ 2 P is a vague prior with a heavy tail. Due to convergence difficulties, an inverse gamma prior could not be used for σ 2 I ; a uniform prior of 0-0.0004 was used instead (S3 Fig). The adaptive MCMC uses Metropolis-Hastings updates with joint normal proposal [29,30,34]. Following standard practice, we use a large number of samples from the MCMC chain with 8 × 10 5 iterations, minimize effects of the initial parameter guess in the MCMC samples by removing a 4% burn-in, and thin the chains for analysis [31]. Visual inspection and the potential scale reduction factor suggest the Markov chains are well mixed and converged [35]. We drive the calibrated model with future forcings to project smooth simulations f(θ, t) of AIS volume loss in SLE. Following Eq (1), the smooth simulations are superimposed with the process noise to account for the internal variability and heteroskedasticity not captured by the model. We superimpose the smooth simulations with the paleo process noise from 240 to 5 kyr BP and superimpose the instrumental process noise from 1961 to the year 2300. We apply a linear regression from the paleo to the instrumental process noise to simulate a smooth transition period from 5 kyr BP to the year 1960. The simulations with process noise represent a probabilistic distribution of AIS volume loss in SLE and the Markov chains represent the empirical estimate of the joint parameter posterior (S3 and S4 Figs).

A full Bayesian inversion improves hindcasts relative to pre-calibration
Moving from the pre-calibration method to a full Bayesian inversion improves the fit of the model to the observational constraints. As perhaps expected, the more sophisticated Bayesian calibration method increases the percentage of samples passing through each individual observational constraint (Table 1). More importantly, we assess the hindcast skill of the model fit with the root-mean-square error (RMSE) between the mean and median model fit and the observational constraints. The Bayesian inversion improves the hindcast skill by reducing the root-mean square error from roughly 5.5 and 3.8 m (mean and median pre-calibration RMSE) to 0.8 m (Bayesian inversion RMSE for both mean and median).

The calibrated model captures MISI during the last interglacial period
The mean hindcast produces enough melt, decrease in radius, and increase in ice flux to suggest the model triggers MISI during the LIG period. MISI occurs when oceanic and atmospheric warming causes the ice flux at the grounding line to increase and the grounding line to retreat unabated until the grounding line reaches an upsloping bed or temperatures cool enough to reform and stabilize the buttressing ice shelf (Fig 1a and 1b) [8,13]. During the LIG, the AIS melts by roughly 2.9 m in SLE (mean) causing a decrease in the ice sheet radius (Fig 1c  and 1d). As this occurs, the water depth and the ice flux at the grounding line increase suggesting retreat of the grounding line (Fig 1e and 1f). The combination of the increased melt, ice flux, and water depth at the grounding line not only demonstrates that the model captures MISI, but that the model can produce enough melt to simulate a loss of the marine ice sheet (i.e., a contribution to global mean sea-level rise of more than three meters; [1]).

The calibrated model partially passes the hindcast test, with evidence for structural discrepancies
The hindcast test reveals some ability to reproduce observational constraints as well as structural model discrepancies. The model does reasonably well at reproducing the width of observational constraint (95% confidence interval) during the LGM, MH, and the instrumental period (Figs 2 and 3). However, our model fit does not reproduce the mean and the top percentage of melting during the LIG. During the LIG, our calibrated model hindcasts roughly 2.9 m, or between -0.1 and 5.9 m (mean and 95% credible interval) of AIS contribution to sea level. These results indicate that our mean hindcast is roughly 26% (1 m) lower than the LIG observed data (Figs 2 and 3a).

Future sea-level contributions diverge from projections considering MICI
Our projection, in the year 2100, produces results that are comparable to previous expert assessments [15][16][17][18]. These previous expert assessments consider the sensitivity or vulnerability of glacier collapse (e.g. Pine Island and Thwaites) within the next century, project changes in surface mass balance and potential rapid ice sheet dynamics, or use a sub-grid interpolation of basal melting at the grounding lines [15][16][17][18]. Our calibrated model (in 2100) has an AIS contribution of roughly 0-0.3 m (90% credible interval). In comparison, the previous expert assessments [15][16][17][18] generally produce a similar or narrower uncertainty range for the 90% confidence/credible interval or low to high estimate (Fig 3e and Table 2). Overall, our median projection is comparable, differing by roughly 3-5 cm (Fig 3e and Table 2).
Additionally, we compare our projections in the year 2100 to some non-model expert assessments [4,9,36]. In contrast to the comparisons with the model expert assessments [15][16][17][18], the non-model expert assessments produce a wider [9,36] or a narrower [4] uncertainty range for the 90% confidence/credible interval. The differences between our results versus the non-model expert assessments are likely due to several reasons. For instance, the non-model expert assessments produce their estimate based on an interpretation of various model studies and scaled to account for dynamical effects, based on formal expert elicitation by conducting expert interviews, and by using maximal kinematic constraints [4,9,36]. None of the non- model expert assessments are calibrated to a RCP8.5 or similar scenario, making comparison to these assessments difficult.
Although our model calibration results in projections that are comparable to previous expert assessments based on models [15][16][17][18], our projection is unable to account for the sizable melt generated from a physically more realistic model considering hydro-fracturing and MICI [8] (Fig 4). Our calibrated model without MICI projects an AIS contribution of roughly 0.09 m (mean) and 0-0.2 (± 1σ) in the year 2100 (Fig 4 and Table 2). In comparison, our results, which neglect MICI, miss 96 and 100% of the ± 1σ projections produced in the physically more realistic model [8] (Fig 4). Furthermore, our mean projection, in 2100, is roughly 1/6 and 1/10 (0.55 and 0.96 m smaller) of the mean produced in the physically more realistic model [8] (Fig 4 and Table 2).

Potential explanations for the discrepancy
The missing key processes and characteristics in the model such as the effect of MICI, the capturing of regional characteristics, and effect of hydro-fracturing of buttressing ice-shelves are potential explanations for this discrepancy. The calibrated model produces discrepancies during periods warmer than today (i.e., LIG and 2100). During the LIG, global mean temperatures were roughly 2˚C above today and evidence indicates a potential collapse of the marine ice sheet [1,37,38]. This rise in temperature is consistent with projections over the next few centuries [37,38]. As discussed in the introduction, melting and calving from warming temperatures can trigger rapid runaway retreat of the ice sheet grounded below sea level, MISI (Fig 1a  and 1b). MISI is highly dependent on the regional characteristics of subglacial basins. For instance, MISI can occur for the West Antarctic ice sheet based on the stability of Pine Island Glacier, Thwaites Glacier, and other glaciers. Yet, high volume losses will be defined by thresholds for the marine part of the East Antarctic ice sheet. A potential mechanism to trigger such losses faster is the MICI mechanism, which can help facilitate MISI in such marine basins and therefore ice sheet disintegration. Despite the models ability to capture MISI, the model neglects the effects of MICI and the ability to capture regional characteristics of the ice sheet. Hence missing MICI as well as the inability to estimate regional characteristics are two potential causes of the reduced melting.

Discussion and caveats
This analysis adopts a relatively simple model structure and neglects several potentially important sources of uncertainty. The modeling choice is motivated by the goal to produce a transparent analysis, but this simplicity also points to caveats and future research needs. For example, we analyze the question whether and how much neglecting MICI can reduce AIS melt during warming periods by comparing our results to projections from a physically more realistic model [8], a model accounting for MICI. A more refined (but also much more involved) approach to quantify this bias would be to calibrate a model with and without MICI. Additionally, we make expert judgments of the prior range for each parameter given that there are only four data points, but 13 parameters. Imposing hard bounds potentially cuts off the tails of the posterior parameter densities, however, providing the model with far too much flexibility given the small number of data points could lead to physically unrealistic results. Furthermore, our analysis neglects the impacts of many other uncertainties (e.g., about the paleo and future forcings). For example, the sea-level curve preceding the LIG might be considerably revised/improved using more recent work (see for example, [39,40]). Despite these caveats, the relative simplicity of the analysis enables the approximation of the effect of missing MICI on AIS melt in the model.

Conclusion
We calibrate a simple AIS model (that does not include a cliff instability mechanism nor is able to capture regional characteristics) with observational constraints over the past 240,000 years using a Bayesian inversion considering the heteroskedastic nature of the data. Using the hindcasts and projections, we compare our results to those from a pre-calibration method and expert assessments with potentially more realistic models. We approximate how neglecting fast processes (i.e., the MICI mechanism) in an AIS model can lead to biases in the AIS hindcasts and projections during warming periods. For the specific example considered, we show how missing MICI produces a lower mean hindcast (roughly 26% or 1 m smaller) during the LIG, a period when the marine ice sheet is suggested to have deglaciated. Additionally, the model is unable to account for roughly 96 and 100% of future AIS contributions predicted by a physically more realistic model accounting for MISI, MICI, and hydro-fracturing yet reproduces the projected median melt in other expert assessments in the year 2100. Overall, accounting for retreat mechanisms can potentially increase warming period AIS melt and reduce model discrepancy.
Supporting Information S1 Text. Discussion of constraints used for calibrating the DAIS model. (PDF) S1