Articles | Volume 13, issue 10
Model description paper
15 Oct 2020
Model description paper |  | 15 Oct 2020

ISSM-SLPS: geodetically compliant Sea-Level Projection System for the Ice-sheet and Sea-level System Model v4.17

Eric Larour, Lambert Caron, Mathieu Morlighem, Surendra Adhikari, Thomas Frederikse, Nicole-Jeanne Schlegel, Erik Ivins, Benjamin Hamlington, Robert Kopp, and Sophie Nowicki

Understanding future impacts of sea-level rise at the local level is important for mitigating its effects. In particular, quantifying the range of sea-level rise outcomes in a probabilistic way enables coastal planners to better adapt strategies, depending on cost, timing and risk tolerance. For a time horizon of 100 years, frameworks have been developed that provide such projections by relying on sea-level fingerprints where contributions from different processes are sampled at each individual time step and summed up to create probability distributions of sea-level rise for each desired location. While advantageous, this method does not readily allow for including new physics developed in forward models of each component. For example, couplings and feedbacks between ice sheets, ocean circulation and solid-Earth uplift cannot easily be represented in such frameworks. Indeed, the main impediment to inclusion of more forward model physics in probabilistic sea-level frameworks is the availability of dynamically computed sea-level fingerprints that can be directly linked to local mass changes. Here, we demonstrate such an approach within the Ice-sheet and Sea-level System Model (ISSM), where we develop a probabilistic framework that can readily be coupled to forward process models such as those for ice sheets, glacial isostatic adjustment, hydrology and ocean circulation, among others. Through large-scale uncertainty quantification, we demonstrate how this approach enables inclusion of incremental improvements in all forward models and provides fidelity to time-correlated processes. The projection system may readily process input and output quantities that are geodetically consistent with space and terrestrial measurement systems. The approach can also account for numerous improvements in our understanding of sea-level processes.

1 Introduction

Reliable projections of local sea-level change, together with robust uncertainties, are a key quantity for stakeholders to shape adequate and cost-effective mitigation and adaptation measures to sea-level rise (Kopp et al.2019). Most regional sea-level projections use a process-based approach, in which all relevant processes are modeled separately and summed up together, including the individual estimates of error, with their spatial signature (Slangen et al.2012; Church et al.2013b; Kopp et al.2014; Jackson and Jevrejeva2016; Kopp et al.2017; Jevrejeva et al.2019). These projections are widely used by coastal planners and stakeholders, as is, for example, demonstrated by the impact of Kopp et al. (2014, 2017) on assessment reports across the United States (Gornitz et al.2019; City of Boston2016; Kopp et al.2016; Kaplan et al.2016; Callahan et al.2017; Dalton et al.2017; Griggs et al.2017; Miller et al.2018; Boesch et al.2018).

In their simplest form, these process based projections can be expressed straightforwardly. We generally refer to these as KOPP14 (Kopp et al.2014)) and write

(1) RSL total θ , ϕ , t = i = 1 n F GRD , i ( θ , ϕ ) B i ( t ) + RSL sterodynamic θ , ϕ , t + RSL GIA θ , ϕ , t .

where RSLtotalθ,ϕ,t is the total projected relative sea-level (RSL) change at time t, latitude θ and longitude ϕ. For all barystatic processes, or processes that change the total ocean mass, the effects of gravity, rotation and deformation (GRD) on local sea level are computed by multiplying the total barystatic contribution Bi(t) by the associated barystatic-GRD fingerprint (abbreviated by “fingerprint” from here on), or FGRD,i(θ,ϕ), which is computed a priori. This procedure is generally used to include the effects of glacier and ice-sheet mass loss, as well as for projected changes in terrestrial water storage (TWS). Note here that our definition of GRD is not completely in line with Gregory et al. (2019), as glacial isostatic adjustment (GIA) is considered as a separate contributor, and the GRD contribution does contribute to global-mean sea-level changes. It is rather in line with the definition of contemporary GRD in Gregory et al. (2019). The effects of sterodynamic sea-level change RSLsterodynamicθ,ϕ,t, which is the sum of global thermosteric expansion and local sea-level changes due to ocean dynamics, is generally included by directly using estimates from Earth system models (ESMs) and atmospheric–oceanic global circulation models (AOGCMs), such as the output of the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al.2009). Finally, the GIA term RSLGIAθ,ϕ,t is generally accounted for using output from a periodically updated global model.

To derive uncertainties for these local projections of sea-level change, the barystatic components Bi are often sampled from a probability distribution found in published probabilistic projections, for example, from expert elicitation projects (e.g., Bamber et al.2019), or other ice-sheet models (DeConto and Pollard2016). The sterodynamic contribution often uses the inter-model spread as a source of the uncertainties. While the basis of each probabilistic projection is similar, each group adds additional components and physics to Eq. (1). For example, in Kopp et al. (2014) and Kopp et al. (2017), a Gaussian process regression model, based on tide-gauge observations, is used to account for the effect of non-climatic vertical land motion. Or in Jackson and Jevrejeva (2016) and Kopp et al. (2017), the GRD effects of ocean dynamics (Richter et al.2013) are explicitly taken into account, with Kopp et al. (2017) computing these effects over the entire projection time series.

One of the key strengths of this approach is its relative simplicity and transparency, as the process from probabilistic estimates of the underlying processes into local sea-level changes is a simple multiplication operation with the respective barystatic-GRD fingerprint. It provides a framework that outputs a probability density function (PDF) for RSL at any desired location, from which the expected sea-level change and its confidence intervals can be derived. This provides both efficient calibration/validation quantities to projections and streamlines incrementally updated projections. In essence, each modular input may be improved separately, so updates are unencumbered by the queueing up of new modules for incorporation into more complex ESMs and AOGCMs.

Recently, a growing body of research indicates that additional processes should be considered in this process-based approach. Inclusion of such processes is critical to improving the quantification of uncertainties in local sea-level change predictions, but they are not directly feasible within the framework of Eq. (1). Below, we highlight some of the key contributors to uncertainty that, until now, have not been considered together in large-scale estimates of sea-level change.

First, in Eq. (1), the multiplication of a barystatic mass contributor Bi(t) with a fingerprint FGRD,i(θ,ϕ), assumes that the fingerprint is constant through time, which is not always the case (Mitrovica et al.2011). Instead, a fingerprint results from feedbacks between the geometry of sea-level components. For example, local sea level depends on the geometry of ice mass loss, so temporal changes in ice geometry will directly translate into local sea-level changes (e.g.,  Larour et al.2017; Mitrovica et al.2018). As a result, this temporal variability not only affects the expected local sea-level changes but also its uncertainties, as the uncertainty of the input mass loss also has a pronounced spatial pattern due to relative limitations in measurement and data interpretation. An example of the inadequacy of temporally constant fingerprints is shown in Fig. 1 for a projection of Greenland's contribution to RSL at 2016 versus 2045 and 2075. Normalized RSL patterns are clearly different between all three times, and the differences are not local to just Greenland but spill over into regions such as North Europe, Alaska and the Canadian arctic.

Figure 1Normalized fingerprints for Greenland in (a) 2016, (b) 2045 and (c) 2075, based on the JPL ISSM experiment 5 simulation that contributed to ISMIP6 (Goelzer et al.2020). The relative sea-level change patterns in the North Atlantic and Arctic oceans vary significantly in both space and time, resulting in very different contributions to local RSL in the years 2016, 2045 and 2075.

Second, covariance in time as well as the covariance between the individual processes is not always small, though it is often considered to be or is approximated by a simple relationship (e.g., Church et al.2013b). Indeed, assuming so could cause a significant misrepresentation of the estimated uncertainties in local sea-level change. For example, Le Bars (2018) showed that most driving factors of sea level are correlated with global-mean temperature changes, and ignoring this inter-process covariance can underestimate uncertainty in local sea-level change. Note that in addition to covariance between processes, the uncertainty in individual processes may also be correlated temporally. Propagating this full spatiotemporal covariance into projections and its uncertainties promotes a better understanding of the spatial and temporal coherence of uncertainties, which could, for example, allow us to assess the likelihood of reaching specific sea levels by 2100 given observed sea-level change during the next 20 years.

Thirdly, recent work on the Antarctic Ice Sheet (AIS) shows a strong coupling between GIA, elastic surface deformation and ice mass loss (Gomez et al.2018; Barletta et al.2018; Larour et al.2019). Such relationships between these processes suggest that any uncertainties in computed ice-sheet histories and solid-Earth properties that propagate into GIA projections (Caron et al.2018) can also feed back into ice-mass-loss projections; thus, considering these processes as independent ignores these couplings. Here, the main problem is that projection frameworks are articulated in terms of changes in mass, while most ice-sheet models, GIA models, TWS evolution models and glacier models are explicitly described in terms of local mass change evolution (or thickness changes, in meters per water equivalent per year; m w.e.yr−1). In order to be able to account for strong couplings, or to even be able to ingest recent modeling results, one needs to propagate the local mass changes and the associated uncertainties into regional sea-level projections. This is particularly relevant now given new modeling runs that have been carried out within large modeling intercomparison projects (MIPs) such CMIP5 and CMIP6, as well as ISMIP6 or GlacierMIP2.

Figure 2Contribution to uncertainty in 100-year extreme warming simulations of AIS and three subregions of the AIS, tested for four different model variables independently. Each probability distribution function represents an ensemble of 800 ISSM ice-sheet forward model runs, sampled using the DAKOTA uncertainty quantification framework (Schlegel et al.2018).


Similarly, additional strong positive feedbacks between ice sheet and ocean dynamics have been evidenced in work from Goldberg et al. (2012, 2018, 2019) and Seroussi et al. (2017), among others. Specifically, these studies suggest that strong coupling between sub-ice-shelf ocean circulation (in particular melt rates) and ice-flow dynamics (in particular, grounding line dynamics and mass transport resulting in modifications of an ice-shelf draft) results in significant retreat of ice streams such as Thwaites Glacier and Pine Island Glacier, as Antarctica's warm circumpolar deepwater is advected close to their grounding line. Other high-frequency processes (at the daily to monthly level) such as ocean tides, and in particular how tidal currents affect water mass properties at ice-sheet marine margins (Padman et al.2018), are critical in understanding how mass loss rates will evolve. This will significantly impact how melt-rate parameterizations are developed to quantify melt rates, especially in the West Antarctic Ice Sheet area (Seroussi et al.2017). Significant work remains in calibrating such melt-rate parameterizations to correctly account for all aforementioned effects. While more work is required in terms of constraining such parameterizations, the impacts of such ice–ocean feedbacks have not been assessed in probabilistic sea-level models (PSLMs).

Finally, in the past decade, extensive work has been carried out to probabilistically characterize components such as GIA (Whitehouse et al.2012; Gunter et al.2014; Caron et al.2018; Melini and Spada2019) or ice-sheet mass balance (Larour et al.2012b, a; Schlegel and Larour2019; Schlegel et al.2013, 2015, 2016, 2018; Smith-Johnsen et al.2020). Substantial understanding of the impact of rheological parameters and ice history on the distribution of bedrock uplift and rate of change in geoid rates has been generated through modeling of GIA. Similarly, for ice-sheet models, significant knowledge has been generated about how the mass balances of both the AIS and Greenland Ice Sheet (GIS) are impacted by surface mass balance (SMB), ice-shelf basal melt, ice–bedrock friction, geothermal heat flux or ice rheology (see, e.g., Fig. 2a). All of these advances need to be fully integrated into new probabilistic projections of sea-level change, and a new approach therefore needs to be envisioned that will allow for such new processes to be accurately modeled.

Indeed, moving from strategies where continental-scale mass changes are sampled and multiplied with the corresponding fingerprint to actually sampling upstream model inputs is important for improving the state of the art. In particular, there is a strong need to fully account for spatial patterns of mass change and their uncertainty (see Fig. 2b–d). This applies to, among others, SMB, basal friction or ice and solid-Earth rheological properties.

For example, Eq. (1) relies on masses that are aggregated at the basin/continental level. However, most ice-sheet models compute high-resolution thickness change patterns that are not aggregated. This aggregation greatly reduces the complexity in representation of model physics and uncertainty propagated at the interface between ice-sheet models and PSLMs. A more comprehensive approach that re-establishes interfaces between forward models and PSLMs is therefore necessary, where model outputs are not aggregated or simplified.

Here, we propose a new framework for sea-level projections that is able to account for all terms in Eq. (1). We improve the existing process-based approach by using the Ice-sheet and Sea-level System Model (ISSM; Larour et al.2012c), which allows for inclusion of forward model physics. It also improves the modeling and sampling of covariance between input processes, both temporally and spatially through the computation of high-resolution barystatic-GRD patterns. The latter feature builds the basis for a geodetically compliant projection system where GRD patterns and their computation are done systematically and which does not introduce biases in the projections.

2 Methods

2.1 Theory

Sterodynamic sea-level changes form a significant contributor to both global-mean sea-level rise and are responsible for large parts of the regional deviations from the global-mean projected changes (e.g., Slangen et al.2012; Church et al.2013b; Slangen et al.2017). Following the CMIP5 conventions, sterodynamic sea-level changes consist of a global-mean thermosteric contribution (variable name zostoga) and a local dynamic contribution (variable name zos) with a zero mean over the oceans. Generally, an ensemble of model runs, either based on multiple models (e.g., Church et al.2013b) or on large-ensemble experiments based on perturbing a single model (for example,  Little et al.2017), can be used to directly sample regional sea-level changes. An alternative approach to generate more samples than model ensemble members is to determine common modes of variability, for example, by extracting the largest empirical orthogonal functions from each model and perturbing the associated principal components (e.g., Thompson et al.2016, Fig. 3).

While sterodynamic effects do not change the total ocean mass, ocean dynamics can be coupled to redistribution of ocean mass, which manifests in ocean-bottom-pressure changes, particularly on shallow shelf seas (Landerer et al.2007). These bottom-pressure changes load the solid Earth below and thus result in GRD effects, which are often referred to as self-attraction and loading (SAL) effects (Ray1998; Stepanov and Hughes2004; Vinogradova et al.2015). These SAL effects could cause several centimeters of additional sea-level rise above the sterodynamic signal in century-scale sea-level projections made by atmosphere–ocean general circulation models (AOGCMs) (Richter et al.2013). By adding the ocean-bottom-pressure changes to the sea-level equation solver, this effect can be incorporated in regional sea-level projections.

As depicted in Eq. (1), in the classical approach, static sea-level fingerprints are computed a priori for each individual process, which typically include glaciers (GLA), the GIS and AIS, and TWS. These fingerprints are subsequently multiplied by the equivalent barystatic contribution, which is often sampled from a PDF and added, together with the sterodynamic and GIA contribution, to obtain local RSL changes and the associated confidence intervals. This method is both transparent and simple, while maintaining computational efficiency due to the fact that the fingerprints do not have to be computed for each sample or time step.

However, several issues arise from this approach, which can be mitigated using a different method. First, it is assumed that the spatial pattern of mass loss is known a priori and does not vary over time. A common approach is to assume that the mass loss is uniformly distributed over the ice sheet, or that it follows the spatial pattern derived over the Gravity Recovery and Climate Experiment (GRACE) period. Jackson and Jevrejeva (2016) quantified the errors induced by assuming a uniform mass loss and found that this bias could be up to 1 and  10 cm for sites distant from and close to centers of mass loss. Furthermore, the approximation of time-invariant fingerprints could lead to biases, when the spatial pattern of mass loss varies over time.

Figure 3Diagram of ISSM's Sea-Level Projection System (ISSM-SLPS) model. The system is driven by requirements from Eq. (2). ISSM-SESAW is the GRD core of the system (in pink). ISSM-SLPS is a combination of ISSM-SESAW and a layer (in green) that handles STR, DSL and GIA inputs, as well as all uncertainty quantification aspects.


In our approach (Fig. 3), ISSM Sea-Level Projection System (ISSM-SLPS) solves for RSL as follows:

(2) RSL ( θ , ϕ , t ) = RSL STR ( t ) + RSL DSL ( θ , ϕ , t ) + RSL GIA ( θ , ϕ , t ) + RSL GRD ( θ , ϕ , t ) .

The first two terms on the right-hand side, i.e., RSLSTR(t)+RSLDSL(θ,ϕ,t), together represent the sterodynamic sea-level change. STR represents the global-mean thermosteric expansion and DSL local sea-level changes due to ocean dynamics. These can be obtained from CMIP results. The GIA contribution to ongoing sea-level change, RSLGIA, is given, for example, by Caron et al. (2018). The last term, RSLGRD, refers to the component of sea-level change due to mass-induced contemporary GRD response of the solid Earth (Gregory et al., 2019), excluding the GIA processes. This implies that viscoelastic deformation is split between long-term timescales and short-term fast rebound of the bedrock uplift, such as observed in West Antarctica (Barletta et al.2018), acting essentially over timescales of 10–100 years. This includes mass transport between the land and the ocean, as well as that due to dynamic ocean circulation. The latter field is provided by CMIP as the ocean-bottom-pressure (OBP) products. Note that GRD associated with land–ocean mass transport is usually termed “sea level fingerprint” (e.g., Mitrovica et al.2009), while the GRD due to OBP variability is termed the “self-attraction and loading” phenomenon (e.g., Ray1998). As we shall see, we unify both of these elements of contemporary GRD sea level in Eq. (3). Note also that the global-mean OBP is removed from the ocean models, since ocean dynamics do not add or remove any mass from/to the ocean. In fact, our projections rely on CMIP5 and CMIP6 “zos” (the sea-surface height change above geoid) and “zostoga” fields for which the Greatbatch correction has been applied (Greatbatch1994).

We compute RSLGRD using ISSM's Solid Earth and Sea-Level Adjustment module (ISSM-SESAW; Adhikari et al.2016). Assuming that all of land ice–water mass change directly modulates the ocean mass, we define a global mass-conserving loading function, Mglobal(θ,ϕ,t), that describes the change in mass per unit area on the solid-Earth surface as follows:

(3) M global ( θ , ϕ , t ) = M land ( θ , ϕ , t ) [ 1 - O ( θ , ϕ ) ] + ρ o [ H OBP ( θ , ϕ , t ) + RSL GRD ( θ , ϕ , t ) ] O ( θ , ϕ ) ,

where the land loading function (with dimensions of mass per unit area) Mland(θ,ϕ,t) is given by Adhikari et al. (2016):

(4) M land ( θ , ϕ , t ) = ρ i [ H AIS ( θ , ϕ , t ) + H GIS ( θ , ϕ , t ) + H GLA ( θ , ϕ , t ) ] + ρ w H TWS ( θ , ϕ , t ) .

Here, ρi is the ice density, ρw is the freshwater density, ρo is the mean density of ocean water, and HOBP is the (ocean) water equivalent height of the ocean-bottom-pressure change. Similarly, HAIS, HGIS and HGLA are the ice height change in the respective cryospheric domains, and HTWS is the freshwater height change in the non-cryospheric land domain. Note that we invoke an ocean function 𝒪(θ,ϕ) in Eq. (3) to ensure mass conservation in the system. This function is equal to 1 over the oceans and 0 everywhere else.

The contemporary mass transport function Mglobal(θ,ϕ,t) loads the underlying solid Earth that is self-gravitating, rotating and viscoelastically compressible. The induced spatial pattern of RSLGRD(θ,ϕ,t) is dictated by the perturbation in Earth's gravitational and rotational potentials and associated viscoelastic deformation of the solid Earth (Farrell and Clark1976; Milne and Mitrovica1998). In the absence of dynamic sea level and meteorologically induced high-frequency signals, the sea-surface height mimics the spatial pattern of the geoid (Gregory et al.2019). Therefore, we may write

(5) RSL GRD ( θ , ϕ , t ) = C ( t ) + G GRD ( θ , ϕ , t ) - B GRD ( θ , ϕ , t ) ,

where GGRD(θ,ϕ,t) and BGRD(θ,ϕ,t) represent the change in geoid and bedrock elevation induced by the loading of the solid Earth (Eq. 3), respectively. Spatial invariant C(t) is invoked to ensure mass conservation in the Earth system, and it may be readily derived by inserting Eq. (5) into Eq. (3) and integrating it over the solid-Earth surface.

Both GGRD(θ,ϕ,t) and BGRD(θ,ϕ,t) appearing in Eq. (5) may be partitioned into two components each: those related to gravitational potential and those to rotational potential. These components can be computed by convolving Mglobal(θ,ϕ,t) with respective Green's functions. These may be defined in terms of surface harmonics with loading Love numbers as coefficients. Given the structure and viscoelastic properties of the solid Earth, these numbers characterize the axisymmetric deformational and gravitational response of Earth to the applied unit surface load. The rotational components depend upon tesseral second-degree loading and tidal Love numbers as well as on the perturbation in Earth's inertia tensor, which in turn depends on Mglobal(θ,ϕ,t). In order to solve for RGRD(θ,ϕ,t), we require an a priori knowledge of Mglobal(θ,ϕ,t), which in turn depends on RGRD(θ,ϕ,t) itself. The system of Eq. (3) and (5) is therefore solved iteratively until a desired solution accuracy is achieved. One key feature of this field is that as ice sheets lose mass, the near-field relative sea level drops, and far-field sea level rises at a much larger rate than the barystatic term for the sake of mass conservation. While theoretical/numerical treatments on the topic are found elsewhere (e.g., Farrell and Clark1976; Mitrovica and Peltier1991; Mitrovica and Milne2003; Spada and Stocchi2007; Adhikari et al.2019), version 1.0 of the SESAW algorithm where RSLGRD is solved for is presented in Adhikari et al. (2016).

Figure 4A 3-D plot of the boundaries used to mesh each continental area of the Earth's surface. Regions include South and North America, Australia, Eurasia and the Pacific, as well as Greenland and Antarctica. In this particular scenario, Greenland has been subdivided into 18 regions along the boundaries defined in Zwally et al. (2012).

2.2 Meshing

SESAW is a mesh-based convolution based on Eq. (2) in Farrell and Clark (1976). As such, it relies on an anisotropic unstructured mesh of the surface of the Earth which is refined according to specific metrics such as distance to the nearest coastline, presence of loads (such as changes in ice thickness or TWS) and the complexity of the coastline. Given the amount of inputs being sampled for in the SLPS system, a systematic approach to refining such a mesh needs to be developed. The main tool for such a refinement is the ISSM implementation of the Bidimensional Anisotropic Mesh Generator (BAMG) anisotropic mesh refiner (Hecht2006a, b). This is a 2-D-based anisotropic mesher which can refine a mesh according to several constraints at the same time: a metric to specify directions along which the mesh resolution needs to be improved, specific vertex or segment positions, in particular vertex positions of the region outlines, and specified mesh resolutions for user-defined locations. Combining these constraints, we develop an approach based on meshing of a set of 2-D continental areas of the Earth, projection of such 2-D meshes onto the 3-D Earth surface and then stitching of the resulting meshes into one seamless global 3-D mesh.

Figure 5A 2-D adaptive mesh of South America using GRACE observations of ice and hydrological mass change (in cm yr−1) from 2003 to 2016. Seismic effects (Richter et al.2019) are not removed in this rendering of Patagonian ice mass loss that was directly taken from Adhikari and Ivins (2016). The Global Self-consistent, Hierarchical, High-resolution Geography database (GHSSH) L1 coarse version coastline is shown in green. Segments of the triangular mesh are plotted in black. The color bar for the thickness changes was saturated at [−1,1] cm yr−1 in order to improve the contrast of the figure given the high mesh resolution.

A plot of the 2-D regions is given in Fig. 4, which include South America, North America, Australia, Eurasia and the Pacific. At the north and south, we have regions defined for Antarctica and Greenland. Greenland itself has been further refined into 18 regions drawn along its main ice divides, following Zwally et al. (2012, Fig. 3). The approach facilitates a direct linkage of models from the existing literature, or potentially from previous ISSM studies such as Seroussi et al. (2017); Schlegel et al. (2018), without having to remesh the entire Earth. This in turn allows for direct comparisons between uncertainty quantification projection results where only one specific region is modified, hence allowing an approach where control runs can be compared against specific variations of an uncertainty quantification projection run.

An example mesh of the South American continent is shown in Fig. 5. This mesh relies on defined vertices for the outline, which match the outline vertices for the Pacific, Antarctica and Eurasian meshes, so that the stitching within a larger 3-D mesh can be done without redundancy in vertices along continental boundaries. In addition, GRACE ice mass trends from 2003 to 2016 (Adhikari and Ivins2016) are provided as a metric to be used for refinement of the mesh, in particular around the Patagonian ice fields. The minimum mesh resolution attained for this mesh is 500 m, and the largest is 1400 km. Finally, Global Self-consistent, Hierarchical, High-resolution Geography database (GHSSH) L1 coarse version (Wessel and Smith1996) was used as a vertex constraint, so that the final mesh perfectly coincides with the coastline dataset (in black). This allows for the most optimum sea-level solution using the SESAW solver.

Once each region has been meshed in 2-D using BAMG, it is projected onto latitude and longitude, and concatenated together to create a 3-D mesh. This is possible because each 2-D mesh relies on the same set of boundaries as shown in Fig. 4. The resulting mesh is shown in Fig. 6 and comprises 38 944 surface elements for 19 486 vertices. For comparison, an equiangular 1×1 grid would require 64 800 vertices, which is 3 times as many as for a coarse grid resolution.

Figure 6A 3-D Earth mesh stitched from 3-D projections of 2-D regional meshes of the following regions: South and North America, Australia, Eurasia and the Pacific, as well as Antarctica and Greenland. GRACE observations of ice mass change (in cm yr−1) from 2003 to 2016 (Adhikari and Ivins2016) are laid over the mesh. The left frame azimuth is 30 with an elevation of 64. The right frame azimuth is 205 with an elevation of 23.

2.3 Sampling and partitioning

In order to sample variables at each time step, our approach is to use a geographical partitioning of the unstructured mesh. An example is shown in Fig. 7, where a range of values from 1 to 5 has been attributed for each vertex (and element) of the mesh, corresponding, respectively, to Antarctica (1), Greenland (2), glaciers (3), the ocean (4) and land (5). For each partition and for each variable that is probabilistically sampled, we define a probability density function (PDF). For normal distributions, for example, this will be done through a mean and standard deviation.

Figure 7Partition vector (values from 1 to 5: 1 for Antarctica, 2 for Greenland, 3 for glaciers around the world, 4 for the ocean and 5 for ice-free land). The partition vector is used to sample probabilistic variables in a geographically consistent way, with PDF distribution moments (mean and standard deviation) defined for each partition area.

The algorithm for sampling through SLPS is explained in Algotithm 1, in the generic case where spatial covariance is available between variables.

Here, t is the time variable (ranging from the years 2019 to 2100, at 1-year intervals), j is the counter for each sample, from 1 to nsamples (in our case, 10 000), VAR is the sampled variable (from one of the SESAW inputs, excluding RSLGIA which is deterministic in our framework), VERTEX is a counter for all vertices in the mesh MESH, PARTITION is the partition vector (for example, ranging from 1 to 5 in Fig. 7), PDF is the joint probability distribution of variables across all geographical locations, DAKOTA is the sample generator in ISSM (Eldred et al.2008; Larour et al.2012b), alpha is the jth sample matrix of scaling factors with size (number of variables, number of partitions), VAR0 the unmodified variable (stored in memory at the beginning of the model run), and SLPS is the sea-level solver, generating RSL for a specific sample of all the probabilistic variables. In this algorithm, the PDF distribution is built by specifying its nature and parameters; e.g., the “type” argument can indicate the choice of a multivariate Gaussian distribution and “pdfspec_arg” can specify the vector of means and covariance matrix of alpha between each other and between partitions.

For this application, we assume that each variable and each partition is independent, and we set the mean of all distributions to 1. This ensures that values of alpha behave as scaling parameters. We use them to directly scale a variable locally, according to which partition area this location geographically belongs. This method is therefore significantly different from the approach in KOPP14, where the entire mass within a certain partition (for example, GIS or AIS) is sampled. Here, the sampling is a scaling of a vectorial field, which therefore preserves the local geographical distribution of a given variable. This is shown in Fig. 8 for a scenario where thinning rates of the GIS are sampled using one geographical partition (corresponding to Fig. 7 partition value of 2, in blue). We display the average thinning rates μ, μ+3σ and μ−3σ (for an arbitrary value of the standard deviation σ=5 %). The structure of the thinning rate as it is sampled is kept intact, implying that the spatial covariance of the variable being sampled across the mesh is kept closely similar across samples and within any given partition.

Figure 8Random sampling of thinning rates across Greenland. (a) GRACE-generated thinning rate pattern at 2005 (in m yr−1). (b) Thinning rate along the AB profile (from a) (in red, representing the average of the PDF) and samples generated at −3σ (blue) and +3σ (yellow) from the average.


Figure 9ISSM-SLPS projections based on Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) Representative Concentration Pathway 8.5 (RCP8.5). For each time step, we sample (10 000 times with Latin hypercube sampling, or LHS) the following inputs: HGIS, HAIS, thermal expansion of the ocean (STR), HTWS and glacier contributions HGLA (see AR5 WG1 chap. 13; Church et al.2013b). Each input's PDF is calibrated using the AR5 5 %–95 % projection confidence interval, similar to Kopp et al. (2014). The resulting global mean sea-level (GMSL) probability distribution function (PDF) distribution is shown in panels (a) (in time) and (b) (at a subset of time steps). The 5 %–95 % confidence interval (likely range, following AR5 definition) is plotted in black in panel (a), along with the temporal mean. Each time step is fully decorrelated from the previous time steps, with this test being used to validate against existing an existing AR5 projection.


Figure 10AR5 calibrated projection of RSL for nine cities around the world from 2007 to 2100. Sampling was carried out for HAIS, HGIS, HGLA, STR and HTWS using mean and standard deviations from AR5 (Church et al.2013a). The patterns for ice thickness are from GRACE 2003–2016 trends (Adhikari and Ivins2016). DSL is fully deterministic, from the CMIP5 Norwegian Earth System Model (NorESM) runs (Bentsen et al.2013). RSLGIA was deterministically set to 0. Each time step was sampled for using 10 000 LHS samples.


2.4 Modularity

The advantage of the partition approach as implemented in SLPS is that various approaches to probabilistic projections can be executed with the same framework. First, as we will show in the next section, the KOPP14 approaches are fully compatible with the SLPS framework. Indeed, fingerprint patterns can be recomputed using local thickness change rate patterns that are spatially constant on the basis of only one partition, such as the entire Greenland or Antarctic ice sheets (contrary to Fig. 12 where Greenland is subdivided). Once several partitions are adopted, however, the refinement in the fingerprint patterns significantly departs from the KOPP14 approach. Second, existing probabilistic assessments for specific components (such as the impact of changes in surface mass balance or basal friction in Antarctica – as shown in Schlegel et al.2018 – on ice thickness changes) can be used directly, using model output (for example, for thickness change rates), or PDF distributions from such model outputs. If the uncertainty quantification was done using a Bayesian framework, the model output statistics can be reused directly (using some type of uniform discrete sampling of each model output), hence replicating a Bayesian-type exploration approach of SLPS without incurring any additional computational cost (meaning, not having to rerun the analysis carried out to compute such model outputs). Third, ISSM modules can be activated upstream of the SLPS solver to push further the boundaries of the uncertainty assessment. For example, an analysis of the impact of SMB variations in one specific region of Antarctica could be carried out using the ice-flow modeling core of ISSM, capable of delivering ice thickness changes directly to the SLPS core. Fourth, these modules can be activated while remaining coupled to other modules. For example, in Larour et al. (2019), it was demonstrated that over centennial timescales, coupling between the elastic uplift of the grounding line and ice-flow-related grounding line migration are key to controlling the retreat of Thwaites Glacier in West Antarctica. Assessing the uncertainty brought by such processes on sea-level rise (SLR) projections would require this coupling to be activated, which could be done (assuming computational costs are still realistic) without modifications to the SLPS framework. Finally, given how closely ISSM can be integrated within Web Server architectures using its native JavaScript interface (Larour et al.2017), SLPS is potentially fully compatible with open-source types of collaborative approaches where inputs from the community could be provided directly to web servers running ISSM in the background to generate model projections without significant investment in a computation core and/or an interface to the latter.

3 Results and discussion

SLPS probabilistic projections were validated using model inputs from the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) (Church et al.2013a). AR5 supplies several projection components in SLR equivalents: the “expansion” term (STR), the “glacier” term (which can be converted into an average thickness change rate for HGLA), “antnet” and “greennet” for net barystatic contribution from the Antarctica and Greenland ice sheets, which can also be converted into an average change rate for HAIS and HGIS, and the “landwater” term for TWS contribution to SLR (which can be converted into an average change rate for HTWS). For each of these terms, AR5 supplies the mean projection and the 5 %–95 % confidence interval. We can use this information to calibrate PDF distributions for thickness change rates at each time step, with the mean of each PDF corresponding to the AR5 mean and the standard deviation calibrated from the 5 %–95 % interval (corresponding to the −1.65σ to 1.65σ interval). Because AR5 does not supply spatial patterns, we rely on GRACE 2003–2016 thickness change rate patterns from Adhikari and Ivins (2016) for HGLA, HAIS and HGIS. For HTWS, we assume a uniform spatial distribution over all the spatial partitions. STR is also considered uniform over all the oceans. DSL is not sampled but rather deterministically set to the DSL term of the CMIP5 NorESM-ME runs (Bentsen et al.2013). GIA is independently sampled (from Caron et al.2018) and probabilistically added as an independent PDF. The sampling is carried out on the partitions described in Fig. 7 with the notable exception that the GIS is further divided into 18 different basins as defined in Zwally et al. (2012) and as plotted in Fig. 12. For each year between 2007 and 2100, 10 000 sample runs of SLPS are carried out (with full geodetic capabilities of the SESAW core). For each partition, samples for the corresponding inputs are generated using a Latin hypercube sampling (LHS) algorithm. The runs were carried out on the Pleiades cluster at the NASA Ames Research Center, on 20 Ivy nodes (20 cores per node for an equivalent 400 cores) over 7 h.

Figure 9 shows projection results for GMSL computed at each time step between 2007 and 2100, and histograms for several time snapshots. We match the mean and 5 %–95 % confidence intervals of AR5 (Fig. 9a) as expected. We also show the evolution of RSL in Fig. 10 for nine cities around the world. We provide the mean and standard deviation for each PDF and show how the sampling of ice-related thickness changes impacts mean and standard deviation. In particular, as expected from the AR5 inputs, we show a marked increase in PDF spreads as time evolves from 2007 to 2100.

Figure 11 shows the impact of using existing statistics from Caron et al. (2018) to include into SLPS. These statistics were evaluated using a Bayesian inversion method based on simulated annealing (Kirkpatrick et al.1983), a variation of the Monte Carlo with Markov chains (MCMC) method (Metropolis and Ulam1949; Metropolis et al.1953). They can be used directly in SLPS either during a standard probabilistic projection run or a posteriori as is the case here. These statistics reflect the statistical fitness to a global GIA dataset composed of paleo-RSL indicators and vertical GPS trends. The impact of the migrating Laurentide isostatic bulge on Norfolk, Virginia, is apparent in Fig. 11, with an offset of 16 cm in the average projection for the city.

Figure 11AR5 calibrated projection of RSL for Norfolk (in brown) versus same projection in which GIA statistics from Caron et al. (2018) are used to account for GIA-induced RSL (in blue). Note that both PDF distributions have standard deviations that are essentially identical within 0.4 % relative difference.


Figure 12 shows results for a different experiment, in which we quantify the impact of refining the amount of partitions used to sample the uncertainty in ice thickness change rates. For the area of Greenland, we use either one partition (blue boundary) or 18 boundaries (brown basins) from the Zwally et al. (2012) dataset. Each basin is delimited by ice divides and thus represents a dynamically coherent area, expected to behave (short of ice divides migrating actively) independently from one another. We rerun an SLPS projection using a similar AR5 setup and display the contribution of ice-related basins to SLR in New York and Hawaii for 1 and 18 partitions respectively. As expected, the mean in PDF distributions are identical for both 1 and 18 partitions. However, the tails are much larger for the one-basin scenario. The relative difference in standard deviations between 1 and 18 basins ranges from −23 % for New York to −34 % for Hawaii. This implies that current probabilistic RSL projections are significantly overestimating (by 20 %–30 %) the width of the “likely” (5 %–95 %) range in ice-melt contribution to RSL.

Figure 12Impact of subsampling the GIS mass (from 1 basin for the whole ice sheet to 18 basins) on barystatic sea-level rise in New York and Hawaii. The distributions are a result of SLPS, where HGIS, HAIS and HGLA were sampled 10 000 times using an LHS algorithm. The means in PDF distributions for both scenarios are identical; however, the tails are much larger for the one-basin scenario. The relative difference in standard deviations between 1 and 18 basins ranges from −23 % for New York to −34 % for Hawaii. This implies that current probabilistic RSL projections could significantly overestimate (20 %–30 %) the “likely” (5 %–95 %) range in ice-melt contribution from glaciers and ice sheets.


This is understandable because of the fact that in a one-partition scenario, variations of ice thickness are dictated by scaling of the local ice thickness change rate mean by an identical scalar for the entire partition, which leads to more extreme values for the contribution to RSL. With finer partitions, basins that have low thickness change rates do not impact RSL as much, and in aggregate the total contribution range varies less. This can be visualized better by taking the example of New York, where following Larour et al. (2017) contributions from South Greenland are almost negligible. This implies that all the basins (and corresponding GRD patterns) in South Greenland will contribute zero variance to the PDF for RSL in New York. This will therefore result in smaller tails for projections that rely on more refined basins. A very similar conclusion was found in Schlegel et al. (2018), where Antarctica had to be subdivided in spatially coherent areas, which were not obvious initially and did not mandatorily map into individual basins. The issue is that the error distribution in model inputs had a specific spatial coherence that had to be respected. Assuming this coherence extended to the entire ice sheet led to significantly larger and unrealistic uncertainty ranges in model outputs. Of course, given differing dynamics in each geographical basin, we cannot assume that the input scaling should be similar (same standard deviation). This will modify the results in Fig. 12. But our aim here is to point out the issue of sub-partitioning as being essential in quantifying the right range of spread in modeled statistical outputs.

This analysis also shows that using SLPS, it is possible to efficiently address the question of how to sample uncertainty in a manner that is consistent with the local behavior of separate basins, glaciers and ice sheets. In Jackson and Jevrejeva (2016), for example, it is shown that the impact of glacier ice thickness variations around the world is significantly different and that relying on one fingerprint alone can lead to significant differences in the projection of glacier contribution (up to several percent). Our approach in SLPS ensures that the GRD contribution is systematically reassessed for each sample, at each time step, and the partitioning of our sampling ensures that we correctly capture the specificity of each glacier/ice/hydrological area and their unique mass change trends. It is to be noted that a similar approach is currently implemented in new instantiations of the KOPP14 projection system based on sampling of glacier projections across the 19 Randolph Glacier Inventory (RGI) areas used in the GlacierMIP results (Hock et al.2019). However, these areas can be very large in spatial extent (such as the low latitudes or north Asian areas) and should be broken down. Our approach scales for any barystatic contributor at any spatial scale (for example, sub-basin or at the glacier level) required by the structure of the error distribution of model inputs.

4 Conclusions

ISSM SLPS is a new sea-level probabilistic projection system which relies on a new partitioning approach to sampling of boundary conditions, forcings and inputs. It is compatible with previous probabilistic frameworks but allows for a more robust integration of state-of-the-art results in the modeling of ice flow in ice sheets and glaciers, sterodynamic sea level, TWS evolution and GIA. It re-establishes temporal correlation in projections where they were previously lacking and allows for better constraints on spatial and temporal covariance in the model inputs. In particular, it is capable of systematically computing geodetically compliant patterns of sea level that are consistent with space and terrestrial measurement systems. The system relies heavily on the use of high-resolution anisotropic meshes and allows for a better interfacing with existing modeling frameworks which operate at higher resolutions and consistently generate changes in mass density patterns around the globe. SLPS has been validated against previous frameworks and is fully backwards compatible. Differences between SLPS and previous approaches have also been shown both in terms of integration of GIA statistics and integration of new high-resolution sampling of ice thickness change patterns in Greenland. This new approach offers a roadmap towards further increasing the complexity and realism of sea-level probabilistic projection frameworks.

Code availability

The ISSM code and its SLPS components are available at (last access: 29 September 2020). The instructions for the compilation of ISSM and SLPS modules are available at (last access: 29 September 2020). The public SVN repository for the ISSM code can also be found directly at (Larour et al.2020) and downloaded using user name “anon” and password “anon”. The version of the code for this study, corresponding to ISSM release 4.17, is SVN version tag number 24683.

Data availability

All datasets used in the projections are freely available in the public domain and are referenced in the text.

Author contributions

EL carried out all the simulations and implemented SLPS into ISSM. He wrote the bulk of the manuscript. LC contributed computations for GIA. MM contributed enhancements to the adaptive meshers in ISSM. All authors contributed to the manuscript in terms of text, figures and comments.

Competing interests

The authors declare that they have no conflict of interest.


We would like to acknowledge the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research for all the high-performance computing and corresponding support. We would also like to acknowledge Ed Zaron for his insightful comments on the manuscript.

Financial support

This research has been supported by the Jet Propulsion Laboratory, California Institute of Technology (prime contract no. 80NMO0018B0004), the NASA Sea-level Change Team (N-SLCT, WBS no. 105393), NASA Cryospheric Sciences (WBS no. 105393), NASA Modeling Analysis and Prediction (MAP, WBS no. 105479) and NASA Earth Surface and Interior (ESI, WBS 105526) programs, as well as the NASA GRACE and GRACE-FO science team programs (WBS nos. 32700 and 106749).

Review statement

This paper was edited by Steven Phipps and reviewed by three anonymous referees.


Adhikari, S. and Ivins, E. R.: Climate-driven polar motion: 2003–2015, Sci. Adv., 2, e1501693,, 2016. a, b, c, d, e

Adhikari, S., Ivins, E. R., and Larour, E.: ISSM-SESAW v1.0: mesh-based computation of gravitationally consistent sea-level and geodetic signatures caused by cryosphere and climate driven mass change, Geosci. Model Dev., 9, 1087–1109,, 2016. a, b, c

Adhikari, S., Ivins, E. R., Frederikse, T., Landerer, F. W., and Caron, L.: Sea-level fingerprints emergent from GRACE mission data, Earth Syst. Sci. Data, 11, 629–646,, 2019. a

Bamber, J. L., Oppenheimer, M., Kopp, R. E., Aspinall, W. P., and Cooke, R. M.: Ice Sheet Contributions to Future Sea-Level Rise from Structured Expert Judgment, P. Natl. Acad. Sci. USA, 116, 11195–11200,, 2019. a

Barletta, V. R., Bevis, M., Smith, B. E., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S. A., Rovira-Navarro, M., Dalziel, I., Smalley, R., Kendrick, E., Konfal, S., Caccamise, D. J., Aster, R. C., Nyblade, A., and Wiens, D. A.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability, Science, 360, 1335–1339, 2018. a, b

Bentsen, M., Bethke, I., Debernard, J. B., Iversen, T., Kirkevåg, A., Seland, Ø., Drange, H., Roelandt, C., Seierstad, I. A., Hoose, C., and Kristjánsson, J. E.: The Norwegian Earth System Model, NorESM1-M – Part 1: Description and basic evaluation of the physical climate, Geosci. Model Dev., 6, 687–720,, 2013. a, b

Boesch, D., Boicourt, W., Cullather, R., Ezer, T., G.E. Galloway, J., Johnson, Z., Kilbourne, K., Kirwan, M., Kopp, R., Land, S., Li, M., Nardin, W., Sommerfield, C., and Sweet, W.: Sea-level Rise: Projections for Maryland 2018, University of Maryland Center for Environmental Science, Cambridge, MD, available at: Rise Projections for Maryland 2018_0.pdf (last access: 29 September 2020), 2018. a

Callahan, J. A., Horton, B. P., Nikitina, D. L., Sommerfield, C. K., McKenna, T. E., and Swallow, D.: Recommendation of Sea-Level Rise Planning Scenarios for Delaware:Technical Report, prepared for Delaware Department of Natural Resources and Environmental Control (DNREC) Delaware Coastal Programs, available at: SLR Technical Report 2017.pdf (last access: 29 September 2020), 2017. a

Caron, L., Ivins, E. R., Larour, E., Adhikari, S., Nilsson, J., and Blewitt, G.: GIA Model Statistics for GRACE Hydrology, Cryosphere, and Ocean Science, Geophys. Res. Lett., 45, 2203–2212, 2018. a, b, c, d, e, f

Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., Merrifield, M., Milne, G., Nerem, R., Nunn, P., Payne, A., Pfeffer, W., Stammer, D., and Unnikrishnan, A.: Sea Level Change Supplementary Material, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, available at: (last access: 29 September 2020), 2013a. a, b

Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., Merrifield, M., Milne, G., Nerem, R., Nunn, P., Payne, A., Pfeffer, W., Stammer, D., and Unnikrishnan, A.: Sea Level Change, chap. 13, Cambridge University Press, Cambridge, UK, New York, NY, USA, 1137–1216,, 2013b. a, b, c, d, e

City of Boston: Boston Research Advisory Group (BRAG) Climate Projections Full Report, Prepared for Climate Ready Boston, available at: (last access: 29 September 2020), 2016. a

Dalton, M., Dello, K., Hawkins, L., Mote, P., and Rupp, D.: The Third Oregon Climate Assessment Report, Oregon Climate Change Research Institute, College of Earth, Ocean and Atmospheric Sciences, Oregon State University, Corvallis, OR, available at: (last access: 29 September 2020), 2017. a

DeConto, R. and Pollard, D.: Contribution of Antarctica to past and future sea-level rise, Nature, 531, 591–597,, 2016. a

Eldred, M. S., Adams, B. M., Gay, D. M., Swiler, L. P., Haskell, K., Bohnhoff, W. J., Eddy, J. P., Hart, W. E., Watson, J.-P., Hough, P. D., and Kolda, T. G.: DAKOTA, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis, Version 4.2 User's Manual, Technical Report SAND 2006-6337, Tech. rep., Sandia National Laboratories, PO Box 5800, Albuquerque, NM 87185, 2008. a

Farrell, W. E. and Clark, J. A.: On Postglacial Sea Level, Geophys. J. Int., 46, 647–667, 1976. a, b, c

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096,, 2020. a

Goldberg, D. N., Little, C. M., Sergienko, O. V., Gnanadesikan, A., Hallberg, R., and Oppenheimer, M.: Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 2. Sensitivity to external forcings, J. Geophys. Res., 117, 1–11,, 2012. a

Goldberg, D. N., Snow, K., Holland, P., Jordan, J. R., Campin, J. M., Heimbach, P., Arthern, R., and Jenkins, A.: Representing grounding line migration in synchronous coupling between a marine ice sheet model and a z-coordinate ocean model, Ocean Model., 125, 45–60,, 2018. a

Goldberg, D. N., Gourmelen, N., Kimura, S., Millan, R., and Snow, K.: How Accurately Should We Model Ice Shelf Melt Rates?, Geophys. Res. Lett., 46, 189–199,, 2019. a

Gomez, N., Latychev, K., and Pollard, D.: A Coupled Ice Sheet-Sea Level Model Incorporating 3D Earth Structure: Variations in Antarctica during the Last Deglacial Retreat, J. Climate, 31, 4041–4054, 2018. a

Gornitz, V., Oppenheimer, M., Kopp, R., Orton, P., Buchanan, M., Lin, N., Horton, R., and Bader, D.: New York City Panel on Climate Change 2019 Report Chapter 3: Sea Level Rise, Annals of the New York Academy of Sciences, 1439, 71–94, 2019. a

Greatbatch, R. J.: A note on the representation of steric sea level in models that conserve volume rather than mass, J. Geophys. Res., 99, 12767++12771,, 1994. a

Gregory, J., Griffies, S., Hughes, C., et al.: Concepts and Terminology for Sea Level: Mean, Variability and Change, Both Local and Global, Surv. Geophys., 40, 1251–1289, 2019. a, b, c

Griggs, G., Arvai, J., Cayan, D., DeConto, R., Fox, J., Fricker, H., Kopp, R., Tebaldi, C., and Whiteman, E.: Rising Seas in California: An Update on Sea-Level Rise Science, California Ocean Science Trust, available at: (last access: 29 September 2020), 2017. a

Gunter, B. C., Didova, O., Riva, R. E. M., Ligtenberg, S. R. M., Lenaerts, J. T. M., King, M. A., van den Broeke, M. R., and Urban, T.: Empirical estimation of present-day Antarctic glacial isostatic adjustment and ice mass change, The Cryosphere, 8, 743–760,, 2014. a

Hecht, F.: BAMG: Bi-dimensional Anisotropic Mesh Generator, Tech. rep., FreeFem++, available at: (last access: 29 September 2020), 2006a. a

Hecht, F.: Maillage 2d, 3d, adaptation, Tech. rep., Universite Pierre et Marie Curie, Paris, 2006b. a

Hock, R., Bliss, A., Marzeion, B., Giesen, R. H., Hirabayashi, Y., Huss, M., Radic, V., and Slangen, A. B. A.: GlacierMIP – A model intercomparison of global-scale glacier mass-balance models and projections, J. Glaciol., 65, 453–467, 2019. a

Jackson, L. P. and Jevrejeva, S.: A probabilistic approach to 21st century regional sea-level projections using RCP and High-end scenarios, Glob. Planet. Change, 146, 179–189, 2016. a, b, c, d

Jevrejeva, S., Frederikse, T., Kopp, R. E., Le Cozannet, G., Jackson, L. P., and van de Wal, R. S. W.: Probabilistic Sea Level Projections at the Coast by 2100, Surv. Geophys., 81, 1–24, 2019. a

Kaplan, M., Campo, M., Auermuller, L., and Herb, J.: Assessing New Jersey's Exposure to Sea-Level Rise and Coastal Storms: A Companion Report to the New Jersey Climate Adaptation Alliance Science and Technical Advisory Panel Report, Prepared for the New Jersey Climate Adaptation Alliance, Rutgers University, New Brunswick, NJ, available at: (last access: 29 September 2020), 2016. a

Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P.: Optimization by Simulated Annealing, Science, 220, 671–680,, 1983. a

Kopp, R., Broccoli, A., Horton, B., Kreeger, D., Leichenko, R., Miller, J., Miller, J., Orton, P., Parris, A., Robinson, D., C.P.Weaver, Campo, M., Kaplan, M., Buchanan, M., Herb, J., Auermuller, L., and Andrews, C.: Assessing New Jersey's Exposure to Sea-Level Rise and Coastal Storms: Report of the New Jersey Climate Adaptation Alliance Science and Technical Advisory Panel, Prepared for the New Jersey Climate Adaptation Alliance. New Brunswick, New Jersey, available at: (last access: 29 September 2020), 2016. a

Kopp, R. E., Horton, R. M., Little, C. M., Mitrovica, J. X., Oppenheimer, M., Rasmussen, D., Strauss, B. H., and Tebaldi, C.: Probabilistic 21st and 22nd century sea-level projections at a global network of tide-gauge sites, Earth's future, 2, 383–406, 2014. a, b, c, d, e

Kopp, R. E., DeConto, R. M., Bader, D. A., Hay, C. C., Horton, R. M., Kulp, S., Oppenheimer, M., Pollard, D., and Strauss, B. H.: Evolving understanding of Antarctic ice-sheet physics and ambiguity in probabilistic sea-level projections, Earth's Future, 5, 1217–1233, 2017. a, b, c, d, e

Kopp, R. E., Gilmore, E. A., Little, C. M., Lorenzo-Trueba, J., Ramenzoni, V. C., and Sweet, W. V.: Usable Science for Managing the Risks of Sea-Level Rise, Earth's Future, 7, 1235–1269,, 2019. a

Landerer, F. W., Jungclaus, J. H., and Marotzke, J.: Ocean Bottom Pressure Changes Lead to a Decreasing Length-of-Day in a Warming Climate, Geophys. Res. Lett., 34, L06307,, 2007. a

Larour, E., Morlighem, M., Seroussi, H., Schiermeier, J., and Rignot, E.: Ice flow sensitivity to geothermal heat flux of Pine Island Glacier, Antarctica, J. Geophys. Res.-Ea. Surf., 117, F04023,, 2012a. a

Larour, E., Schiermeier, J., Rignot, E., Seroussi, H., Morlighem, M., and Paden, J.: Sensitivity Analysis of Pine Island Glacier ice flow using ISSM and DAKOTA, J. Geophys. Res., 117, F02009,, 2012b. a, b

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, F01022,, 2012c. a

Larour, E., Ivins, E. R., and Adhikari, S.: Should coastal planners have concern over where land ice is melting?, Sci. Adv., 3, e1700537,, 2017. a, b

Larour, E., Cheng, D., Perez, G., Quinn, J., Morlighem, M., Duong, B., Nguyen, L., Petrie, K., Harounian, S., Halkides, D., and Hayes, W.: A JavaScript API for the Ice Sheet System Model (ISSM) 4.11: towards an online interactive model for the cryosphere community, Geosci. Model Dev., 10, 4393–4403,, 2017. a

Larour, E., Seroussi, H., Adhikari, S., Ivins, E., Caron, L., Morlighem, M., and Schlegel, N.: Slowdown in Antarctic mass loss from solid Earth and sea-level feedbacks, Science, 364, eaav7908,, 2019. a, b

Larour, E., Morlighem, M., and Seroussi, H.: Ice-Sheet and Sea-Level System Model, svn repository, available at:, last access: 29 September 2020. a

Le Bars, D.: Uncertainty in Sea Level Rise Projections Due to the Dependence Between Contributors, Earth's Future, 6, 1275–1291,, 2018. a

Little, C. M., Piecuch, C. G., and Ponte, R. M.: On the Relationship between the Meridional Overturning Circulation, Alongshore Wind Stress, and United States East Coast Sea Level in the Community Earth System Model Large Ensemble, J. Geophys. Res.-Oceans, 122, 4554–4568,, 2017. a

Melini, D. and Spada, G.: Some remarks on Glacial Isostatic Adjustment modelling uncertainties, Geophys. J. Int., 218, 401–413, 2019. a

Metropolis, N. and Ulam, S.: The Monte Carlo method, J. Amer. Stat. Associ., 44, 335–341, 1949. a

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E.: Equation of State Calculations by Fast Computing Machines, The J. Chem. Phys., 21, 1087–1092,, 1953. a

Miller, I., Morgan, H., Mauger, G., Newton, T., Weldon, R., Schmidt, D., Welch, M., and Grossman, E.: Projected Sea Level Rise for Washington State - A 2018 Assessment, A collaboration of Washington Sea Grant, University of Washington Climate Impacts Group, University of Oregon, University of Washington, and US Geological Survey, available at: (last access: 29 September 2020), 2018. a

Milne, G. and Mitrovica, J.: Postglacial sea-level change on a rotating Earth, Geophys. J. Int., 133, 1–19, 1998. a

Mitrovica, J. and Milne, G.: On post-glacial sea level: I. General theory, Geophys. J. Int., 154, 253–267, 2003. a

Mitrovica, J. and Peltier, W.: On post-glacial geoid subsidence over the equatorial ocean, J. Geophys. Res., 96, 20053–20071, 1991. a

Mitrovica, J. X., Gomez, N., and Clark, P. U.: The Sea-Level Fingerprint of West Antarctic Collapse, Science, 323, 753,, 2009. a

Mitrovica, J. X., Gomez, N., Morrow, E., Hay, C., Latychev, K., and Tamisiea, M. E.: On the Robustness of Predictions of Sea Level Fingerprints, Geophys. J. Int., 187, 729–742,, 2011. a

Mitrovica, J. X., Hay, C. C., Kopp, R. E., Harig, C., and Latychev, K.: Quantifying the Sensitivity of Sea Level Change in Coastal Localities to the Geometry of Polar Ice Mass Flux, J. Climate, 31, 3701–3709,, 2018. a

Padman, L., Siegfried, M. R., and Fricker, H. A.: Ocean Tide Influences on the Antarctic and Greenland Ice Sheets, Rev. Geophys., 56, 142–184,, 2018. a

Ray, R. D.: Ocean self‐attraction and loading in numerical tidal models, Mar. Geodes., 21, 181–192,, 1998. a, b

Richter, A., Groh, A., Horwath, M., Ivins, E., Marderwald, E., Hormaechea, J. L., Perdomo, R., and Dietrich, R.: The Rapid and Steady Mass Loss of the Patagonian Icefields throughout the GRACE Era: 2002–2017, Remote Sens., 11, 909,, 2019. a

Richter, K., Riva, R., and Drange, H.: Impact of Self-Attraction and Loading Effects Induced by Shelf Mass Loading on Projected Regional Sea Level Rise, Geophys. Res. Lett., 40, 1144–1148,, 2013. a, b

Schlegel, N.-J. and Larour, E. Y.: Quantification of surface forcing requirements for a Greenland Ice Sheet model using uncertainty analyses, Geophys. Res. Lett., 46, 9700–9709,, 2019. a

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M., and Box, J. E.: Decadal-scale sensitivity of Northeast Greenland ice flow to errors in surface mass balance using ISSM, J. Geophys. Res.-Ea. Surf., 118, 667–680,, 2013. a

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M., and Box, J. E.: Ice discharge uncertainties in Northeast Greenland from boundary conditions and climate forcing of an ice flow model, J. Geophys. Res.-Ea. Surf., 120, 29–54,, 2015. a

Schlegel, N.-J., Wiese, D. N., Larour, E. Y., Watkins, M. M., Box, J. E., Fettweis, X., and van den Broeke, M. R.: Application of GRACE to the assessment of model-based estimates of monthly Greenland Ice Sheet mass balance (2003–2012), The Cryosphere, 10, 1965–1989,, 2016. a

Schlegel, N.-J., Seroussi, H., Schodlok, M. P., Larour, E. Y., Boening, C., Limonadi, D., Watkins, M. M., Morlighem, M., and van den Broeke, M. R.: Exploration of Antarctic Ice Sheet 100-year contribution to sea level rise and associated model uncertainties using the ISSM framework, The Cryosphere, 12, 3511–3534,, 2018. a, b, c, d, e

Seroussi, H., Nakayama, Y., Menemenlis, D., Larour, E., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199, 2017. a, b, c

Slangen, A. B. A., Katsman, C. A., van de Wal, R. S. W., Vermeersen, L. L. A., and Riva, R. E. M.: Towards Regional Projections of Twenty-First Century Sea-Level Change Based on IPCC SRES Scenarios, Clim. Dynam., 38, 1191–1209,, 2012. a, b

Slangen, A. B. A., Adloff, F., Jevrejeva, S., Leclercq, P. W., Marzeion, B., Wada, Y., and Winkelmann, R.: A Review of Recent Updates of Sea-Level Projections at Global and Regional Scales, Sur. Geophys., 38, 385–406,, 2017. a

Smith-Johnsen, S., Schlegel, N.-J., de Fleurian, B., and Nisancioglu, K. H.: Sensitivity of the Northeast Greenland Ice Stream to geothermal heat, J. Geophys. Res.-Earth Surf., 125, e2019JF005252,, 2020. a

Spada, G. and Stocchi, P.: SELEN: A Fortran 90 program for solving the “sea-level equation”, Comput. Geosci., 33, 538–562, 2007. a

Stepanov, V. N. and Hughes, C. W.: Parameterization of ocean self-attraction and loading in numerical models of the ocean circulation, J. Geophys. Res.-Oceans, 109, C03037,, 2004. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the experiment design, B. Am. Meteorol. Soc., 93, 485–498,, 2012. a

Thompson, P. R., Hamlington, B. D., Landerer, F. W., and Adhikari, S.: Are long tide gauge records in the wrong place to measure global mean sea level rise?, Geophys. Res. Lett., 43, 10403–10411, 2016. a

Vinogradova, N. T., Ponte, R. M., Quinn, K. J., Tamisiea, M. E., Campin, J.-M., and Davis, J. L.: Dynamic Adjustment of the Ocean Circulation to Self-Attraction and Loading Effects, J. Phys. Oceanogr., 45, 678–689,, 2015. a

Wessel, P. and Smith, W. H. F.: A Global, Self-Consistent, Hierarchical, High-Resolution Shoreline Database, J. Geophys. Res.-Sol. Ea., 101, 8741–8743,, 1996. a

Whitehouse, P. L., Bentley, M. J., Milne, G. A., King, M. A., and Thomas, I. D.: A new glacial isostatic adjustment model for Antarctica: calibrated and tested using observations of relative sea-level change and present-day uplift rates, Geophys. J. Int., 190, 1464–1482, 2012. a

Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland Drainage Systems, Tech. rep., GSFC Cryospheric Sciences Laboratory, available at: (last access: 29 September 2020), 2012. a, b, c, d

Short summary
ISSM-SLPS is a new projection system for future sea level that increases the resolution and accuracy of current projection systems and improves the way uncertainty is treated in such projections. This will pave the way for better inclusion of state-of-the-art results from existing intercomparison efforts carried out by the scientific community, such as GlacierMIP2 or ISMIP6, into sea-level projections.