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

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 Jevrejeva, 2016; 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 Boston, 2016; 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

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\text{RSL}}_{\text{total}}\left(\mathit{\theta },\mathit{\varphi },t\right)=& \sum _{i=\mathrm{1}}^n{F}_{\text{GRD},i}(\mathit{\theta },\mathit{\varphi })\cdot B_i\left(t\right)\\ & +{\text{RSL}}_{\text{sterodynamic}}\left(\mathit{\theta },\mathit{\varphi },t\right)\\ & +{\text{RSL}}_{\text{GIA}}\left(\mathit{\theta },\mathit{\varphi },t\right).\end{array}\end{array}$$

where RSL${}_{\text{total}}\left(\mathit{\theta },\mathit{\varphi },t\right)$ 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 B_i(t) by the associated barystatic-GRD fingerprint (abbreviated by “fingerprint” from here on), or $F_{\text{GRD},i}(\mathit{\theta },\mathit{\varphi })$, 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 RSL${}_{\text{sterodynamic}}\left(\mathit{\theta },\mathit{\varphi },t\right)$, 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 RSL${}_{\text{GIA}}\left(\mathit{\theta },\mathit{\varphi },t\right)$ 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 B_i 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 Pollard, 2016). 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 B_i(t) with a fingerprint $F_{\text{GRD},i}(\mathit{\theta },\mathit{\varphi })$, 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 1 Normalized 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⁻¹). 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 2 Contribution 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 Spada, 2019) or ice-sheet mass balance (Larour et al., 2012b, a; Schlegel and Larour, 2019; 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 (Ray, 1998; Stepanov and Hughes, 2004; 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 3 Diagram 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:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{RSL}(\mathit{\theta },\mathit{\varphi },t)& ={\text{RSL}}_{\text{STR}}\left(t\right)+{\text{RSL}}_{\text{DSL}}(\mathit{\theta },\mathit{\varphi },t)\\ & +{\text{RSL}}_{\text{GIA}}(\mathit{\theta },\mathit{\varphi },t)+{\text{RSL}}_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t).\end{array}\end{array}$$

The first two terms on the right-hand side, i.e., RSL${}_{\text{STR}}\left(t\right)+{\text{RSL}}_{\text{DSL}}(\mathit{\theta },\mathit{\varphi },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, RSL_GIA, is given, for example, by Caron et al. (2018). The last term, RSL_GRD, 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., Ray, 1998). 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 (Greatbatch, 1994).

We compute RSL_GRD 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, $M_{\text{global}}(\mathit{\theta },\mathit{\varphi },t)$, that describes the change in mass per unit area on the solid-Earth surface as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{M}_{\text{global}}(\mathit{\theta },\mathit{\varphi },t)& =M_{\text{land}}(\mathit{\theta },\mathit{\varphi },t)[\mathrm{1}-\mathcal{O}(\mathit{\theta },\mathit{\varphi }\left)\right]\\ & +{\mathit{\rho }}_{\text{o}}\left[H_{\text{OBP}}\right(\mathit{\theta },\mathit{\varphi },t)\\ & +{\text{RSL}}_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)\left]\mathcal{O}\right(\mathit{\theta },\mathit{\varphi }),\end{array}\end{array}$$

where the land loading function (with dimensions of mass per unit area) $M_{\text{land}}(\mathit{\theta },\mathit{\varphi },t)$ is given by Adhikari et al. (2016):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{M}_{\text{land}}(\mathit{\theta },\mathit{\varphi },t)& ={\mathit{\rho }}_{\text{i}}\left[H_{\text{AIS}}\right(\mathit{\theta },\mathit{\varphi },t)+H_{\text{GIS}}(\mathit{\theta },\mathit{\varphi },t)\\ & +H_{\text{GLA}}(\mathit{\theta },\mathit{\varphi },t)]+{\mathit{\rho }}_{\text{w}}{H}_{\text{TWS}}(\mathit{\theta },\mathit{\varphi },t).\end{array}\end{array}$$

Here, ρ_i is the ice density, ρ_w is the freshwater density, ρ_o is the mean density of ocean water, and H_OBP is the (ocean) water equivalent height of the ocean-bottom-pressure change. Similarly, H_AIS, H_GIS and H_GLA are the ice height change in the respective cryospheric domains, and H_TWS 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 $M_{\text{global}}(\mathit{\theta },\mathit{\varphi },t)$ loads the underlying solid Earth that is self-gravitating, rotating and viscoelastically compressible. The induced spatial pattern of RSL${}_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)$ is dictated by the perturbation in Earth's gravitational and rotational potentials and associated viscoelastic deformation of the solid Earth (Farrell and Clark, 1976; Milne and Mitrovica, 1998). 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

$$\begin{array}{}\text{(5)}& {\displaystyle }{\text{RSL}}_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)=C\left(t\right)+G_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)-B_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t),\end{array}$$

where $G_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)$ and $B_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },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 $G_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)$ and $B_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },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 $M_{\text{global}}(\mathit{\theta },\mathit{\varphi },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 $M_{\text{global}}(\mathit{\theta },\mathit{\varphi },t)$. In order to solve for $R_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },t)$, we require an a priori knowledge of $M_{\text{global}}(\mathit{\theta },\mathit{\varphi },t)$, which in turn depends on $R_{\text{GRD}}(\mathit{\theta },\mathit{\varphi },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 Clark, 1976; Mitrovica and Peltier, 1991; Mitrovica and Milne, 2003; Spada and Stocchi, 2007; Adhikari et al., 2019), version 1.0 of the SESAW algorithm where RSL_GRD is solved for is presented in Adhikari et al. (2016).

Figure 4 A 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 (Hecht, 2006a, 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 5 A 2-D adaptive mesh of South America using GRACE observations of ice and hydrological mass change (in cm yr⁻¹) 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⁻¹ 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 Ivins, 2016) 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 Smith, 1996) 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 $\mathrm{1}\times \mathrm{1}{}^{\circ }$ grid would require 64 800 vertices, which is 3 times as many as for a coarse grid resolution.

Figure 6 A 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⁻¹) from 2003 to 2016 (Adhikari and Ivins, 2016) 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 7 Partition 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 RSL_GIA 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 8 Random sampling of thinning rates across Greenland. (a) GRACE-generated thinning rate pattern at 2005 (in m yr⁻¹). (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 9 ISSM-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: H_GIS, H_AIS, thermal expansion of the ocean (STR), H_TWS and glacier contributions H_GLA (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 10 AR5 calibrated projection of RSL for nine cities around the world from 2007 to 2100. Sampling was carried out for H_AIS, H_GIS, H_GLA, STR and H_TWS using mean and standard deviations from AR5 (Church et al., 2013a). The patterns for ice thickness are from GRACE 2003–2016 trends (Adhikari and Ivins, 2016). DSL is fully deterministic, from the CMIP5 Norwegian Earth System Model (NorESM) runs (Bentsen et al., 2013). RSL_GIA 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 H_GLA), “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 H_AIS and H_GIS, and the “landwater” term for TWS contribution to SLR (which can be converted into an average change rate for H_TWS). 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 H_GLA, H_AIS and H_GIS. For H_TWS, 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 Ulam, 1949; 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 11 AR5 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 12 Impact 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 H_GIS, H_AIS and H_GLA 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.