Reducing uncertainties in projections of Antarctic ice mass loss

Abstract

Climate model projections are often aggregated into multi-model averages of all models participating in an intercomparison project, such as the Coupled Model Intercomparison Project (CMIP). The “multi-model” approach provides a sensitivity test to the models' structural choices and implicitly assumes that multiple models provide additional and more reliable information than a single model, with higher confidence being placed on results that are common to an ensemble. A first initiative of the ice sheet modeling community, SeaRISE, provided such multi-model average projections of polar ice sheets' contribution to sea-level rise. The SeaRISE Antarctic numerical experiments aggregated results from all models devoid of a priori selection, based on the capacity of such models to represent key ice-dynamical processes. Here, using the experimental setup proposed in SeaRISE, we demonstrate that correctly representing grounding line dynamics is essential to infer future Antarctic mass change. We further illustrate the significant impact on the ensemble mean and deviation of adding one model with a known bias in its ability of modeling grounding line dynamics. We show that this biased model can hardly be identified from the ensemble only based on its estimation of volume change, as ad hoc and untrustworthy parametrizations can force any modeled grounding line to retreat. However, tools are available to test parts of the response of marine ice sheet models to perturbations of climatic and/or oceanic origin (MISMIP, MISMIP3d). Based on recent projections of Pine Island Glacier mass loss, we further show that excluding ice sheet models that do not pass the MISMIP benchmarks decreases the mean contribution and standard deviation of the multi-model ensemble projection by an order of magnitude for that particular drainage basin.

Introduction

The contribution of the Antarctic ice sheet to sea-level rise (SLR) has steadily increased during the last two decades. In the early 1990s, the amount of snow falling over the ice sheet was more or less balanced by the total coastal discharge. Today, the ice sheet loses mass at a rate of $\approx$ 80 Gt yr${}^{-\mathrm{1}}$, equivalent to $\approx$ 0.2 mm yr${}^{-\mathrm{1}}$ of the global eustatic SLR . Proximal geological evidence shows that the western part of the Antarctic ice sheet may have collapsed during warm periods of the late Pleistocene . Such collapses have probably been driven by an unstable retreat of the marine-based regions (i.e., underlying bedrock below sea level), characterized by a retrograde bed slope. The underlying process, marine ice sheet instability (MISI), is supported by theoretical and numerical results . After MISI has been initiated, it could lead to a collapse of the contemporary west Antarctic ice sheet and have the potential to cause sea level to rise by $\approx$ 3.3 m , leading to a drastic impact on human societies . However, conditions for the initiation of such a collapse, and rate of retreat remains poorly known . The potential for MISI underscores the urgent need for reliable projections of Antarctic mass balance in order to conceive efficient regional and global adaptation strategies.

Current projections for mean sea-level rise in 2100 range from 0.28 to 0.98 m depending on the Representative Concentration Pathways (RCP) scenarios, and the contribution of ice sheets represent about a third of the total projected SLR . However, this likely range excludes the possibility of a collapse of west Antarctica. Since the latest IPCC Assessment Report (AR5), new modeling initiatives tend to show that both Pine Island and Thwaites glaciers may have initiated MISI . Significant progress in the ability of marine ice sheet models to reproduce observed dynamical changes in coastal regions has led to these novel results . Still, ice sheet models have not reached the level of development that other models, simulating other components of the climate system, have reached. Antarctic and Greenland ice sheet model ensembles, in particular, remain in their infancy. Only one attempt has been produced so far, namely the SeaRISE initiative, which has been extensively reported in three pivotal papers by and . Results of all participating models were aggregated into unweighted model averages to produce SLR projections. However, the confidence in related projections remains low because of the unproven ability of many participating models to cope with coastal dynamics .

Parallel to the SeaRISE initiative, specific model intercomparison exercises (MISMIP and MISMIP3d

MISMIP: Marine Ice Sheet Model Intercomparison Project.

) have been designed to improve our understanding of grounding line dynamics (i.e., dynamics of the limit between the grounded ice sheet and the downstream floating ice shelf). These initiatives have led to the formulation of requirements regarding physics and numerical approaches to adequately simulate the flow of coastal outlet glaciers in contact with the ocean . In this respect, proposed a multi-model intercomparison to evaluate the response of Pine Island Glacier (PIG) to changes at the grounding line, based on models meeting these MISMIP and MISMIP3d requirements.

In this paper, we assess the origin of uncertainty in recent ice sheet model projections of Antarctic sea-level contribution for PIG, based on SeaRISE and results due to , and , guided by MISMIP and MISMIP3d. We further evaluate the potential bias introduced by models that are limited by marine ice sheet physics and reassess SLR projections for Pine Island Glacier based on MISMIP-tested models. We clearly demonstrate the effect of abandoning the “one model, one vote” approach . A brief inventory of the physics implemented in common ice sheet models is presented in Sect. . In Sect. , we compare the SeaRISE sample with a simplified ice sheet model (SISM) and demonstrate that a proper representation of grounding line dynamics is essential in reducing uncertainties. Finally, a global ensemble analysis for PIG basin is presented.

Ice sheet models inventory/variety

Stokes equations and approximations

The basic problem in ice sheet modeling is to solve the gravity-driven flow of an incompressible and nonlinear viscous ice mass, further extended with a constitutive equation relating stresses to strain rates, i.e., $${\mathit{\tau }}_{ij}=\mathrm{2}\mathit{\eta }{D}_{ij},$$ where $\mathit{\tau }$ is the deviatoric stress tensor and $D_{ij}$ are the components of the strain rate tensor. The effective viscosity $\mathit{\eta }$ is then expressed as $$\mathit{\eta }=\frac{\mathrm{1}}{\mathrm{2}}E{A}^{-\mathrm{1}/n}{D}_e^{(\mathrm{1}-n)/n},$$ where $D_e$ is the strain-rate invariant. Models use a temperature-dependent coefficient $A$, and set $n$ $=$ 3, according to Glen's flow law. $A$ is possibly adjusted through an enhancement factor $E$ that classically accounts for the anisotropic behavior of ice . The velocity (and pressure field) of an ice body is computed by solving the Stokes problem: $$\begin{array}{}& & \text{div}\mathit{u}=\mathrm{0},& & \text{div}\mathit{\tau }-\text{grad}p+{\mathit{\rho }}_i\mathit{g}=\mathrm{0},\end{array}$$ where $p$ is the isotropic pressure and $\mathit{g}$ the gravitational acceleration.

Apart from the boundary conditions, which are discussed below, this model represents the most complete mathematical description of ice sheet dynamics and is commonly called a full-Stokes model. Owing to the considerable computational effort, approximations to these equations are often used, such as higher-order, shallow-shelf and shallow-ice approximations. These approximations involve dropping terms from the momentum balance equations as well as simplifying the strain rate definitions and boundary conditions. Higher-order Blatter–Pattyn-type models consider the hydrostatic approximation in the vertical direction by neglecting vertical resistive stresses . A particular case of this type of model is a depth-integrated hybrid model, combining both membrane and vertical shear stress, and of comparable accuracy to the Blatter–Pattyn model . Vertical shearing terms are included in the calculation of the effective viscosity, but the force balance is simplified. A further approximation, known as the shallow-shelf approximation (SSA), is obtained by neglecting vertical shear .

However, the earliest and most common approximation in large-scale ice dynamics simulations is the shallow-ice approximation (SIA). This approximation incorporates only vertical shear stress gradients opposing the gravitation drive, which is valid for an ice mass with a small aspect ratio (i.e., thickness scale much smaller than length scale) in combination with a significant traction at the bedrock. Its main advantage is that all stress and velocity components are locally determined. The approximation is not valid for key areas such as ice divides and grounding lines , since it excludes membrane stress transfer across the grounding line . The non-validity of SIA is overcome by some models through the use of grounding line flux or grounding line migration parametrizations based on solutions derived from matched asymptotics .

Boundary conditions

We will not list all boundary conditions of thermomechanically coupled ice sheet models, but focus on those that are of importance for grounding line migration. These pertain to the initialization of the ice sheet and conditions at the contact of the ice sheet with the ocean boundary.

Initialization

Initialization of ice sheet models to reproduce the current ice sheet state is commonly done through long-term paleo simulations (paleo spin-up). This has the advantage of establishing a reasonable temperature regime within the ice column . However, reproducing current ice sheet geometry and velocities remains of limited accuracy. To circumvent this pitfall, inverse methods have recently been introduced. Basal drag (or ice viscosity) is inferred by minimizing the misfit between observed and modeled surface velocity e.g., or observed surface elevation . For thermomechanically coupled simulations, a steady-state temperature field is used , which nevertheless has a limiting impact on short-term projections .

The marine boundary

It has long been hypothesized that grounding line migration may provoke unstable behavior when the ice sheet rests on a retrograde bed slope below sea level , leading to a potential collapse of marine-based areas . However, despite major developments in numerical ice sheet modeling, the majority of state-of-the-art models in the 1990s and 2000s did not exhibit significant grounding line retreat when ocean forcing (through ice shelf melting) was applied. The response of Antarctic ice sheet models to a warmer climate led to a higher volume due to atmospheric precipitation increase, but not to a mass loss due to inland migration . The incapacity of numerical ice sheet models to cope with grounding line migration has been seriously challenged by remote sensing observations at the end of the 1990s showing grounding line retreat and substantial mass loss in the Amundsen Sea sector e.g.,. However, ice sheet models were still unable to cope with such rapid changes when IPCC AR4 was released .

A verification of ice sheet models became feasible due to a boundary layer theory developed by , who showed that in the absence of lateral buttressing, grounding line positions are unique and stable on a downward-sloping bedrock. The theory also confirmed the MISI hypothesis in the case of a retrograde bed (unstable grounding line positions), although marine ice sheets are not unconditionally unstable in two horizontal dimensions as buttressing may stabilize grounding line over a retrograde bed . Two model intercomparison exercises have subsequently been organized to verify whether numerical ice sheet models produce results in agreement with the boundary layer theory MISMIP, MISMIP3d;. It further allowed requirements for numerical models to be identified, in order to cope with grounding line migration. Basic conclusions of both intercomparisons are that in order to resolve grounding lines, the inclusion of membrane stresses across the grounding line is required at a sufficiently small grid size ($<$ 500 m), even when a subgrid interpolation of the grounding line ($<$ 5 km) is preferred. An exception is the use of a grounding line parametrization based on the boundary layer theory due to Schoof, which has been successfully implemented in a series of models . This works well for coarse spatial resolutions, but the short-term transient response remains questionable when compared to other approaches , especially since the theory has been developed for the steady-state case.

In a more recent paper, further scrutinized the results from the ice2sea MISMIP3d intercomparison to demonstrate that a clear distinction in the response to marine forcing could be related to the complexity of the model physics. The study shows that at least higher-order approximations are necessary to accurately simulate the flow across the grounding line, as the presence of vertical shearing in the force budget softens the effective viscosity at the grounding line, leading to a faster response on short timescales. The result is clearly different for SSA models that are stiffer at the grounding line and seem to overestimate the contribution to SLR. This has also been confirmed by a model intercomparison of Pine Island Glacier .

Essential characteristics of Antarctic SeaRISE models together with SISM. More details on the numerics can be fund in .

Name	Resolution	Flow	Initialization
	(km)	equation
SICOPOLIS	10	L1Lx	spin-up
Potsdam	15	L1Lx	quasi-steady-state
PennState3D	20	L1Lx/Heuristic	spin-up
AIF	40	L1Lx	present-day geometry – enhancement factor adjustment
UMISM	20	SIA	spin-up
ISSM	3	L1Lx	data assimilation
SISM	15	SIA	quasi-steady-state

Description of SISM (simplified ice sheet model)

To demonstrate the importance of the proper inclusion of a marine boundary in large-scale ice sheet models, we developed a simple (but in terms of marine conditions – wrong) ice sheet model. The simplified ice sheet model (SISM) is a numerical ice sheet model based on the physics inherent to well-known ice sheet models e.g.,. It is a two-dimensional vertically integrated model, solving the Stokes equations according to SIA. The time-dependent evolution of the ice sheet is based on mass conservation and the model is initialized through 1000 (and 100) years of surface relaxation, starting from the BEDMAP2 data set . Details of the model and the model runs are given in Appendix A.

Contrary to the SeaRISE experiments, climate forcing is not applied and the present climate conditions are retained during the whole run (see Appendix A for details). Note that models with a similar degree of physical complexity in the description of ice flow have been included in the SeaRISE multi-model ensemble .

Grounding line dynamics are not explicitly included in SISM. However, melting at the grounding line is introduced by subtracting the amount of basal melt from the surface mass balance at the last grounded grid point. Ice thickness becomes zero when the ice thickness in that grid point – determined from ice advected from upstream and the local mass balance – becomes zero (or negative). Therefore, grounding line retreat is purely due to melting and not due to any physical process operating at the grounding line. Hence, a marine ice sheet instability (retreat of the grounding line on a retrograde slope in absence of melt perturbation and significant buttressing) is not simulated with this simplified model .

Using SISM, we perform a number of the Antarctic SeaRISE experiments and investigate the impact of including a model with a known bias on the ensemble projection.

Sea-level projections

SeaRISE and PIG model ensembles: approaches

SeaRISE ensemble

The SeaRISE initiative led to the first attempt to evaluate multi-ice sheet models ensembles. It is important to note that at the time the experiments were designed, circa 2008, SeaRISE's primary goal was to investigate the sensitivity of ice sheet models to external forcing. Its baseline hypothesis presumed that there was no “best” ice sheet model around and that ensemble modeling would potentially lead to a better understanding of ice sheet models . Six models participated in the SeaRISE modeling of the Antarctic ice sheet, with a large variety of approximations to the Stokes equations and different treatments with respect to implementing grounding line dynamics (see Table ). More details on the physics and numerics of SeaRISE models can be found in .

The SeaRISE experiments all start from an initial present-day ice sheet, which is built up using either a paleo spin-up or assimilation methods. Perturbations in boundary conditions are then imposed for 500 years and compared to a control run to remove the long-term drift. Climate forcing experiments refer to the ensemble mean of AR4 A1B changes in temperature and precipitation being imposed for 94 years and held constant at the values of year 94 for the remainder of the 500-year runs. An amplification factor of 1, 1.5 and 2, respectively, is applied in order to simulate warmer climate scenarios (experiments C1, C2 and C3, respectively). Subsequently, basal sliding perturbations are implemented through a uniform increase of basal sliding (amplified by a factor 2, 2.5 and 3, respectively, for experiments S1, S2 and S3, respectively) and the sensitivity of Antarctic ice shelves to sub-ice shelf melt was performed by applying a uniform melt rate at the base of floating ice (2, 20 and 200 m yr${}^{-\mathrm{1}}$, experiments M1, M2 and M3, respectively). These sensitivity experiments (or combinations of them) were further used to evaluate the dynamic contribution of Antarctica to sea-level rise for the 21st century under the various RCP scenarios , with estimated median contributions ranging from 0.07 m for the low-emission RCP2.6 scenario and 0.09 m for the strongest RCP8.5 . While the reliability of such projections has been questioned , this has not been further evaluated and discussed so far.

PIG ensemble

While SeaRISE focussed on modeling the whole Antarctic ice sheet, a number of studies have simulated the effect of ice shelf melting at the basin scale. Pine Island Glacier and Thwaites Glacier in particular have shown considerable contemporary grounding line retreat and thinning . present a first comprehensive modeling of Pine Island Glacier (PIG) based on a SSA model and using assimilation methods for initialization. Although this particular model did not participate in any MISMIP intercomparison, its physics (SSA) and spatial resolution around the grounding line (down to 140 m) make this model compliant with MISMIP recommendations.

propose a model intercomparison of PIG based on three models of varying complexity, i.e., an SSA model, a higher-order model and a full-Stokes model. These three models also took part in the ice2sea MISMIP3d intercomparison and produced verified results at the spatial resolution used in the PIG intercomparison. Their approach consists of computing an initial state as close as possible to the current geometry and surface velocities using assimilation methods. Melting and calving perturbations are further applied. They show that the response of PIG is mainly driven by the bedrock topography rather than the type and the amplitude of the perturbation and further conclude that PIG is probably already engaged in a MISI. The study finally estimates the contribution of PIG to SLR over the next 20 years ranging from 3.5 to 10 mm. Finally, use a higher-order model over PIG to simulate its dynamical response to marine forcing over the next 50 years using a higher-order model, with spatial resolutions down to 500 m at the grounding line. This particular model did not perform the MISMIP experiments, but as with the model, physics and numerical implementation are conform to MISMIP and MISMIP3d recommendations.

All models presented above will be compared below for the PIG basin. However, we start the analysis with an evaluation of the importance of marine processes on ice sheet response on a pan-Antarctic scale.

Change in grounded area vs. contribution to SLR after 200 years for all models participating in the SeaRISE Antarctic experiments. For reference, a loss of 4.0469 $\times$ 10${}^{14}$ m${}^{\mathrm{3}}$ of ice above flotation equates 1 m rise in mean global sea level. Blue shading corresponds to the sliding experiments with increasing darkness related to an increase in forcing (S1, S2 and S3). Similarly, the lightest red to the darkest one corresponds to melting experiments M1, M2 and M3.

[Figure omitted. See PDF]

Evolution of the Antarctic grounded area as computed by the five models which participated in the SeaRISE experiment M3 (a–e) and similar results obtained by SISM (f). Colors correspond to the time of ungrounding.

[Figure omitted. See PDF]

Grounding line migration in the SeaRISE ensemble

Figure  displays the contribution to SLR after 200 years as a function of the change in grounded area. It clearly highlights the fact that large contributions to SLR always go along with significant changes in the extension of the grounded ice sheet. In other words, having models able to cope with grounding line dynamics is a prerequisite before establishing projections of upper bound dynamic contribution of the Antarctic ice sheet to SLR. The evolution of the grounded area as a response of SeaRISE experiment M3 is presented in Fig. . Despite the drastic and unrealistic perturbation (200 m yr${}^{-\mathrm{1}}$ melt rate, designed to approximate a sudden collapse of all ice shelves), the response of the participating models varies widely, from a limited grounding line retreat to almost a complete collapse of all the marine sectors within a period of 200 years. Moreover, amongst the models presenting a significant retreat, the impacted regions are different with significant differences in grounding line retreat rates.

Large differences in model response are essentially due to two factors: models that correctly implement melting under the ice shelves will fail to produce a significant retreat if the grounding line area is not properly sampled (spatial resolution below 500 m), when using a physical approximation based on SSA, or lacking a parametrization of grounding line dynamics based on the boundary layer theory due to . This failure has been clearly illustrated by and . On the other hand, models that implement melting at the grounding line, i.e., the last grounded grid point, melt grounded ice away, thereby mimicking grounding line retreat. The result is unphysical in both implementation (since melting occurs under the ice shelves) and reaction (spatial resolution and/or physical model are unappropriate). The SISM model illustrates this perfectly, as grounding line retreat in this model is not due to ice-dynamical processes at the grounding line, but due to ice being melted away at the grounded line. The retreat rates produced by this model are within the range produced by SeaRISE, due to the fact that several models within SeaRISE implement grounding line melt in a similar fashion.

Impact of SISM on the SeaRISE ensemble

Since the number of models participating in the Antarctic SeaRISE experiments is rather limited, we may expect that adding a model (e.g., SISM) to the sample will significantly impact on the ensemble mean projections, thereby questioning its relevance. Its effect is illustrated in Fig.  and Table .

Global sea-level increase (cm) projected by SISM after a 1000-year spin-up, together with mean and standard deviation projected by the SeaRISE models extended with SISM (described in Sect. ) or SeaRISE models alone (in brackets) for each experiment at 100 and 200 years. Note that similar numbers are presented in but a different sea-level conversion factor was used (4.0 $\times$ 10${}^{14}$ m${}^{\mathrm{3}}$ equivalent to 1 m mean sea level).

	100		Standard	200		Standard
	SISM	Mean	deviation	SISM	Mean	deviation
S1	11.4	17.1 (18.1)	6.7 (6.8)	20.5	30.1 (31.6)	12.1 (12.4)
S2	16.9	22.2 (23.1)	8.4 (8.8)	30.1	38.5 (39.9)	15.0 (15.9)
S3	22.2	27.2 (28.0)	12.9 (14.0)	39.4	48.1 (49.6)	20.7 (20.9)
M1	17.4	8.6 (6.8)	5.6 (4.0)	28.9	17.4 (15.1)	9.5 (8.6)
M2	104.9	74.6 (68.5)	27.7 (26.1)	151.5	133.9 (130.4)	42.0 (46.0)
M3	465.6	363.9 (343.5)	301.7 (332.7)	758.3	575.1 (538.5)	446.0 (488.4)

Evolution of the contribution to SLR for all the models participating in SeaRISE experiments S1, S2, S3, M1 M2 and M3 (gray lines). SLR contribution computed by SISM for similar perturbations are presented in black, after a 1000- and 100-year spin-up (continuous and dashed line, respectively).

[Figure omitted. See PDF]

Compared to other models, the contribution to SLR with SISM is close to the SeaRISE ensemble mean for sliding experiments and is amongst the largest for melting perturbations, but it is not a striking outlier. As a reminder, SISM is based on simple model physics, isothermal and a surface mass balance set to present-day conditions and not evolving following any RCP scenario. Taking SISM into account in the ensemble unweighted mean leads to two distinct impacts. When considering melt perturbation, adding SISM to the ensemble usually increases both the mean and standard deviation of the ensemble projections. The increase in mean is substantial, up to 20 % for experiments M2 after 100 years. We can anticipate that adding a biased model which would present a limited capacity of grounding line retreat would lead to a decrease of the ensemble mean contribution to SLR together with an increase in the related standard deviation, as the sample size increases. The particular case of sliding experiments and experiment M3 is instructive: the projected contribution of SISM is fortuitously close to the SeaRISE ensemble mean. Including the SISM in the ensemble mean projections slightly affects the mean but also decreases the standard deviation. Ironically, in the particular situation where a biased model projection is coincidentally close to the ensemble mean, introducing such a model may be wrongly interpreted as improving the confidence in the ensemble projection.

The SeaRISE experiments were rerun with SISM, starting from a different spinup (100 years instead of 1000 years; Fig. ). Despite significant differences (several hundreds of meters in ice thickness in coastal areas) between both geometries, the prognostic runs are hardly affected in terms of SLR contribution over the next 200 years, indicating that an initial state close to the observed ice sheet geometry is a weak constraint when large perturbations are applied.

Ensemble analysis on PIG

In view of the small SeaRISE sample, we extended the sample with recent regional studies (basin-scale), focused on Pine Island Glacier . Since the most significant changes in grounding line position and mass balance are currently observed over PIG , this drainage basin appears to be the most appropriate region to evaluate the impact of model physics/numerics on SLR projections. Amongst these studies, argue that PIG is already experiencing MISI and that forthcoming mass change projected by models is relatively similar irrespective of the perturbation amplitude, even for an almost complete collapse of the current ice shelf. Therefore, this implies that comparison with SeaRISE experiments is feasible, despite the difference in melt perturbations between the various studies. However, this comparison remains qualitative as most of the models we compare use different setups that are certainly responsible for a part of the spread between models.

Figure  presents the evolution of the cumulated contribution of PIG drainage basin to SLR for the period 2000–2050 according to the SeaRISE M3 experiment and according to , and . As mentioned above, estimations from SeaRISE range from a very limited retreat of the grounded line (e.g., Fig. a) and relatively low contribution to SLR (below 5 mm cumulated in 2050) to an extremely high discharge of 3 mm yr${}^{-\mathrm{1}}$ and a collapse of the entire drainage basin within a century (e.g., Fig. d). As expected, SISM is amongst the models predicting the highest contribution for PIG, but this model result stays within the whole sample range.

A striking feature of Fig.  is that all projections due to , and occupy a limited range compared to the full range of the SeaRISE sample, with SLR contribution between 2.3 and 18.8 mm by 2040, compared to 2.8 and 146.4 mm, respectively.

Discussion

Most models in the SeaRISE sample have a coarse spatial grid size ($>$ 10 km, see Table ), which – despite the physical approximations – do not sample grounding line dynamics as stipulated in the MISMIP intercomparison . Only one model uses the parametrization based on the boundary layer theory due to , and has been tested against MISMIP, exhibiting a behavior similar to models that do capture grounding line dynamics at high spatial resolution. While such parametrized models are probably less reliable for transient effects, they capture the essence of grounding line migration and stability . The basin-scale simulations for PIG are performed with models that capture grounding line dynamics at the spatial resolution required and with appropriate physics. Most models did participate in the MISMIP intercomparison at the spatial resolution used in the PIG analysis.

Evolution of the cumulated contribution of Pine Island Glacier to SLR until 2050 as computed by models participating in SeaRISE experiment M3 (blue lines), SISM for the SeaRISE M3 forcing (purple line), together with estimations from (black lines), (red lines) and (orange lines) for the same region. Models are ordered according to their complexity: SIA (dotted lines with triangles), SSA (dashed lines with squares), heuristic approach (dashed lines with circles), L1Lx (dashed line) and full-Stokes (continuous). Starting time of SeaRISE experiments is 2004 . Starting times of experiments computed by and correspond to the acquisition year of the surface velocities used for inversion, 1996 and 2008, respectively. As detailed in , the starting time of their experiments corresponds to the last grounding line measurements available, i.e., completed in 2011 . Cumulated contribution was offset to zero in 2009. The light gray (dark gray) envelope encompasses the SeaRISE ensemble (MISMIP-verified ensemble, respectively).

[Figure omitted. See PDF]

Models that capture grounding line dynamics are within the dark gray envelope in Fig.  (associated mean of 8.6 mm and standard deviation of 4.9 mm in 2040), compared to the light-gray envelope that represents the SeaRISE sample (mean of 50.1 mm and standard deviation of 67.5 mm in 2040). Not only do the MISMIP-verified models occupy a smaller range than the full sample but a distinction between the physical representation of each of the models tends to appear. It is not expected that a model according to full Stokes would show the same results as a SSA model, since the physical model is different. Such a distinction clearly appears for models that capture grounding line migration as was already shown in based on an ideal ice sheet geometry. For instance, SSA grounding lines move faster in reaction of a perturbation, because they are stiffer at the grounding line (i.e., the viscosity of the grounding line is higher). Higher-order models (and in the limit full Stokes models) produce a slower response, as the viscosity at the grounding line is lower due to the inclusion of vertical shearing in the stress tensor. Such differentiation is not captured whenever the spatial resolution at the grounding line is too coarse and obliterates the characteristics of the physical model. Results computed by with a higher-order model are in agreement with the estimation computed with the higher-order model used in , despite differences in setup and perturbation. Only the estimations produced by with a SSA model are different to the response of the SSA model in . Finally, the model based on the parametrized approach has the highest contribution to SLR of the sample of models that capture grounding line migration. Although this comparison remains qualitative because boundary conditions and perturbations differ from one modeling experiment to the other, this result is in line with the results of MISMIPs, which tend to demonstrate that application to a “real” case seems to endorse the conclusions of “idealized” case simulations. However, it would probably deserve specific controlled experiments (i.e., similar forcing, domain, initial geometry …) to better apprehend the inter-model spread.

It is peculiar to note that the models due to exhibit marine ice instability (and presumably also the case for and ) and their response is to a large extent indifferent to the amplitude of the perturbation applied. Yet, their contribution to SLR on the timescales considered is smaller than the majority of the models that were used in SeaRISE that did not capture any MISI. This poses serious questions as to whether the inherent complexity of an ice sheet model (thermomechanics, sliding, surface mass balance) is decisive in the process of representing ice sheet response to marine forcing. This issue will definitely become important when ice sheet models will be fully coupled to ocean models at a spatial resolution that complies with both systems. Furthermore, coupling could exhibit a series of other feedbacks between the rate of sub-shelf melting and changes in the sub-shelf cavity shape, which are currently unknown.

Conclusions

The SeaRISE initiative has been the first multi-ice-sheet model ensemble projection to evaluate the future contribution of Antarctica to SLR. Results of all participating model results were taken into account, irrespective of the inherent difference in complexity between the models. A similar approach is used in AOGCMs (Atmosphere-Ocean General Circulation Models) community . However, it is probably simpler to evaluate the ability of ice sheet models to perform adequately when compared to AOGCMs, as fewer key processes are at play. In any case, and whatever the component of the modeled climate component, first order processes must be taken into account . As an example, it sounds particularly inappropriate to use an atmospheric model unable to compute a radiation balance to make any projections of mean surface temperature. Similarly, to compute projections of Antarctic contribution to SLR, ice flow models have to be evaluated on their ability of incorporating grounding line dynamics. If this process is not implemented within a model, any of its projections pertaining to coastal regions are unreliable, even on decadal timescales. Furthermore, solely based on the evolution of the modeled ice sheet, it may be hard to discriminate whether the projection is reasonable or not. Indeed, ad-hoc parametrizations can force any model to retreat, but the lack of physics makes any projection of the retreat and contribution to SLR untrustworthy. Owing to the small number of ice sheet models, including biased models has a strong effect on the mean, as well as on the dispersion of the results.

Benchmark experiments to evaluate the ability of models to cope with grounding line dynamics have been recently developed and others will emerge (MISMIP$+$). Ice sheet models should be evaluated using these benchmarks before being applied on actual cases. Such an approach has been followed by for PIG. Taking into account only a selection of models with appropriate physics and numerics to compute grounding line dynamics significantly reduces the spread of the projected contribution to SLR, reinforcing our confidence in projecting future glacier evolution. Initiatives to produce new multi-ensemble models will undoubtedly be launched in the near future (e.g., ISMIP6

http://www.climate-cryosphere.org/activities/targeted/ismip6

). Their ability to decrease uncertainties will strongly depend on whether inappropriate models (i.e., unvalidated grounding line dynamics) will be included or not.

SISM description and model setup

SISM (simplified ice sheet model) is an isothermal two-dimensional ice sheet model based on the shallow-ice approximation (SIA), using Glen's flow law as a constitutive equation. The vertical mean horizontal ice velocities in the grounded ice sheet are calculated from the local ice geometry e.g., $$\mathit{u}={\mathit{u}}_{\text{b}}-\frac{\mathrm{2}}{n+\mathrm{2}}A\left(\mathit{\rho }g{)}^n{H}^{n+\mathrm{1}}\right|\mathrm{\nabla }h{|}^{n-\mathrm{1}}\mathrm{\nabla }h,$$ where $\mathit{u}$ $=$ ($u$, $v$) are the horizontal flow velocities, ${\mathit{u}}_{\text{b}}$ is the velocity due to basal sliding, $n$ $=$ 3 is the flow-law exponent, $A$ $=$ 5 $\times$ 10${}^{-17}$ Pa${}^{-n}$ yr${}^{-\mathrm{1}}$ is the flow-rate factor, $\mathit{\rho }$ $=$ 910 kg m${}^{\mathrm{3}}$ is the ice density, $H$ (m) is the ice thickness and $h$ (m) is the surface elevation. Basal sliding follows a typical Weertman sliding law, and is defined as $${\mathit{u}}_{\text{b}}=A_{\text{s}}{\mathit{\tau }}_{\text{d}}^m,$$ where ${\mathit{\tau }}_{\text{d}}$ $=$ $-\mathit{\rho }gH\mathrm{\nabla }h$ is the driving stress (Pa), $A_{\text{s}}$ $=$ 3 $\times$ 10${}^{\mathrm{9}}$ (m yr${}^{-\mathrm{1}}$ Pa${}^{-m}$) a sliding factor and $m$ $=$ 2 . Ice-sheet evolution is based on mass conservation, written as a diffusive equation, i.e., $$\frac{\partial H}{\partial t}=\dot{a}-\mathrm{\nabla }\left(\mathit{u}H\right)=\dot{a}+\mathrm{\nabla }D\mathrm{\nabla }h,$$ where $\dot{a}$ is the surface mass balance (m yr${}^{-\mathrm{1}}$), and $D$ spatially varying diffusion coefficients, defined by $D$ $=$ $-\mathit{u}H/\mathrm{\nabla }h$.

The model is numerically solved on a finite-difference grid with a spatial resolution of 15 km. BEDMAP2 data are used as basic input for the model. The ice sheet equation (Eq. ) is solved as a diffusion equation and diffusion coefficients are calculated on a staggered Arakawa B-grid. Since it is an isothermal model, sliding was allowed on the whole domain, but given the very low driving stresses in the interior, its impact on overall ice dynamics remains limited.

Initialization of the model is based on a relaxation of the surface elevation, starting from the BEDMAP2 present-day ice sheet geometry, for a period of 1000 and 100 years, respectively. Surface mass balance is obtained from and , based on the output of a regional atmospheric climate model for the period 1980 to 2004, and calibrated using observed mass balance rates. During these runs, Eq. () was solved for the grounded ice mask. Starting from this initialization, the model was run forward another 200 years to determine the model drift for the forcings.

Similar to the SeaRISE experiments, the model is forced according to three basal sliding and three basal melt scenarios and for a period of 200 years, during which the basal sliding factor is multiplied by 2, 2.5 and 3, respectively, and basal melting at the grounding line is set to 2, 20, and 200 m yr${}^{-\mathrm{1}}$, respectively. Since the model does not have ice shelves, melting was applied at the last grounded grid point that is in contact with the ocean. Whenever ice thickness in that particular grid point becomes afloat or zero for a bedrock below sea level, the grid point becomes ocean and is not taken up in the grounded mask anymore. This simple, and profoundly wrong, mechanism allows the grounding line to retreat. Grounding line advance is, however, not possible and not anticipated given the significant melting factors applied. Similar mechanisms of grounding line retreat were also applied for some of the models within the SeaRISE intercomparison.

Acknowledgements

This study was funded by the French National 960 Research Agency (ANR) under the SUMER (Blanc SIMI 6) 2012 project referenced as ANR-12-BS06-0018. This research is also a contribution to the Belgian Federal Research Project “IceCon” (Constraining Antarctic ice mass change) BELSPO, contract number SD/CA/06A. Edited by: G. H. Gudmundsson