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 5 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 10 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 15 systems. The approach can also account for numerous improvements in our understanding of sea-level processes. Copyright statement. The research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Copyright 2020. All rights reserved.


Introduction
Reliable projections of local sea-level change, together with robust uncertainties, are a key quantity for stakeholders to shape 20 adequate and cost-effective mitigation and adaptation measures to sea-level rise ). 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. (2014Kopp et al. ( , 2017 on assessment reports 25 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 most simple form, these process-based projections (we generally refer to these as KOPP14, in reference to Kopp et al. (2014)) can be expressed as: ( 1) and Jevrejeva (2016) and Kopp et al. (2017), the GRD effects of ocean dynamics (Richter et al., 2013) are explicitly taken into 50 account, with Kopp et al. (2017) computing these effects over the entire projection time series.
One of the key strengths of this approach is how simple and transparent it is, 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 55 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.
Yet recently, a growing body of research indicates that additional processes should be considered in this process-based approach. Indeed, inclusion of such processes is critical to improving the quantification of uncertainties in local sea-level 60 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 GRD,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 ( Caron et al., 2018) can also feed back into ice-mass-loss projections, thus considering these processes as independent ignores 85 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 explicitely described in terms of local mass change evolution (or thickness changes, in m/yr water equivalent). 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 90 large Modeling Intercomparison Projects (MIPs) such CMIP5 and CMIP6, as well as ISMIP6 or GlacierMIP2.
Similarly, additional strong positive feedbacks between ice sheet and ocean dynamics have been evidenced in work from among others, Goldberg et al. (2012Goldberg et al. ( , 2018Goldberg et al. ( , 2019 and Seroussi et al. (submitted). 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 95 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 Sheeet area (Seroussi et al., submitted). Significant work remains in calibrating 100 such melt rate parameterizations to correctly account for all afore-mentioned effects. While more work is required in terms of constraining such parameterizations, the impact 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 105 et al., 2012b, a;Schlegel et al., 2013Schlegel et al., , 2015Schlegel et al., , 2016Schlegel et al., , 2018. 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 modelling of GIA. Similarly, for ice sheet models, significant knowledge has been generated about how the mass balances of both 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 these advances need to be fully integrated into new probabilistic projections of 110 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 e.g. Fig. 2b-d), This applies to among 115 others SMB, basal friction, or ice and solid-Earth rheological properties.
For example, Eq. 1 relies on masses that are agregated at the basin/continental level. However, most ice sheet models compute high-resolution thickness change patterns that are not agregated. This agregation greatly reduces the complexity in representation of model physics and uncetainty propagated at the interface between ice-sheet models and PSLMs. A more comprehensive approach that reestablishes interfaces between forward models and PSLMs is therefore necessary, where model 120 outputs are not agregated or simplified. Model ice dynamics (ISSM-ICE) runs, sampled using the ISSM-DAKOTA uncertainty quantification framework. (Schlegel et al., 2018) 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 covariances between input processes, both temporally and spatially through the computation of high-resolution barystatic-GRD patterns. The latter feature builds 125 the basis for a geodetically compliant projection system where GRD patterns and their computation is done systematically and does not introduce biases in the projections.

Theory
Sterodynamic sea-level changes form a significant contributor to both global-mean sea-level rise and are responsible for large 130 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-135 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 (eg. 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 140 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.

145
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 Greenland and Antarctic Ice Sheets and terrestrial water storage (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 owing to the fact that the fingerprints do 150 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 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 155 1 cm 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.
In our approach ( Figure 3) ISSM Sea-Level Projection System (ISSM-SLPS) solves for RSL as follows: The first two terms on the right hand side, i.e. RSL ST R (t)+RSL DSL (θ, φ, t), together represent the sterodynamic sea level ( Barletta et al., 2018), acting essentially over time scales of 50-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 "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 equation (3). Note also that the 170 global-mean OBP is removed from the ocean models, since ocean dynamics don't add or remove any mass from/to the ocean.
In fact, our projections rely on CMIP5 and CMIP6 fields 'zos' (the sea-surface height change above geoid) and 'zostoga' 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; . Assuming that all of land ice/water mass change directly modulates the ocean mass, we define a global mass conserving loading 175 function, M global (θ, φ, t), that describes the change in mass per unit area on the solid Earth surface as follows: where the land loading function (with dimensions of mass per unit area) M land (θ, φ, t) is given by : 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) the respective cryospheric domains, and H T W S is the freshwater height change in the non-cryospheric land domain. Note that we invoke an ocean function O(θ, φ) in equation (3) to ensure mass conservation in the system. This function is equal to 1 over the oceans, 0 everywhere else.
The contemporary mass transport function M global (θ, φ, t) loads the underlying solid Earth that is self-gravitating, rotating, 185 and viscoelastically compressible. The induced spatial pattern of RSL GRD (θ, φ, 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

190
where G GRD (θ, φ, t) and B GRD (θ, φ, t) represent the change in geoid and bedrock elevation induced by the loading of the solid Earth (equation 3), respectively. Spatial invariant C(t) is invoked to ensure mass conservation in the Earth system, and it may be readily derived by inserting equation (5) into equation (3) and integrating it over the solid Earth surface.
Both G GRD (θ, φ, t) and B GRD (θ, φ, t) appearing in equation (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 global (θ, φ, t)

195
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 global (θ, φ, t). In order to solve for R GRD (θ, φ, t), we require an a priori knowledge of M global (θ, φ, t), which in 200 turn depends on R GRD (θ, φ, t) itself. The system of equations (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 farfield 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 205 presented in .

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 210 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 2D 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 215 constraints, we develop an approach based on meshing of a set of 2D continental areas of the Earth, projection of such 2D meshes onto the 3D Earth surface and then stitching of the resulting meshes into one seamless global 3D mesh.
A plot of the 2D regions is given in Fig. 4, which include South America, North America, Australia, Eurasia and the Pacific regions. At the North and South, we have regions defined for Antarctica, and Greenland. Greenland itself has been further refined into 18 regions drawn along the main ice divides of Greenland, following (Zwally et al., 2012, Fig 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 2D using BAMG, it is projected onto latitude and longitude, and concatenated together to create a 3D mesh. This is possible because each 2D mesh relies on the same set of boundaries as shown in Fig. 4. The resulting mesh is shown in Fig. 6, and comprises 38944 surface elements for 19,486 vertices. For comparison, an equi-angular

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 240 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.
The algorithm for sampling through SLPS is explained below, in the generic case where spatial covariances are available between variables. where t is the time variable (ranging from year 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., 265 2008;Larour et al., 2012b), alpha is the j th 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), 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' specify the vector of means and covariance matrix of alpha between each other and between partitions.

270
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 275 using one geographical partition (corresponding to Fig. 7 partition value of 2, in blue). We display the average thinning rate µ, µ + 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.

Modularity
The advantage of the partition approach as implemented in SLPS is that various approaches to probabilistic projections can be executed with 280 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 sheet (contrary to Fig. 12 where Greeland 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 (Schlegel et al.,285 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.

290
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 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, accepted), 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.

Results and Discussion
SLPS probabilistic projections were validated using model inputs from the Intergovernmental Panel on Climate Change World Fifth Assessment Report (IPCC AR5) (Church et al., 2013a). AR5 supplies several projection components in SLR equivalents: the 'expansion' term (ST R), 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 305 and the 'landwater' term for TWS contribution to SLR (which can be converted into an average change rate for HT W S ). For each of these terms, AR5 supplies the mean projection, and the 5-95% percentile 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  for HGLA, HAIS and HGIS. For HT W S , we assume a uniform spatial distribution over all the spatial partitions. ST R 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  Chapter 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 GMSL PDF distribution is shown in a) (in time) and b) (at a sub-set of time steps). The 5-95% confidence interval (likely range, following AR5 definition) is plotted in black in a), along with the temporal mean. Each time step is fully decorrelated from the previous time steps, this test being used to validate against existing an existing AR5 projection Fig. 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 320 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.  . DSL is fully deterministic, from the CMIP5 NorESM runs (Bentsen et al., 2013). RSLGIA was deterministically set to 0. Each time step was sampled for using 10,000 LHS samples. chains (MCMC) method (Metropolis and Ulam, 1949;Metropolis et al., 1953). They can be used directly in SLPS, either during a standard 325 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.  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 1 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 335 probabilistic RSL projections are significantly overestimating (by 20-30%) the width of the "likely" (5-95%) range in ice-melt contribution to RSL. This is understandable because of the fact that in a 1 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 Figure 12. Impact of sub-sampling 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 mean in PDF distributions for both scenarios are identical, however the tails are much larger for the 1 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.
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 at 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 345 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 point 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 350 that is consistent with the local behaviour of separate basins, glaciers, 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 355 is currently implemented in new instantiations of the KOPP14 projection system based on sampling of glacier projections across the 19 Randold 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 Latitues or North-Asia areas) and should be broken down. Our approach scales for any barystatic contributor, at any spatial scale (example, sub-basin, or at the glacier level) required by the structure of the error distribution of model inputs.
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 reestablishes temporal correlation in projections where they were previously lacking, and allows for better constraints on spatial and temporal covariances in the model inputs. In particular, it is capable of systematically computing geodetically compliant patterns of sea-level that are consistent with space and terrestial 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 which 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 370 realism of sea-level probabilistic projection frameworks.
Code availability. The ISSM code and its SLPS components are available at http://issm.jpl.nasa.gov. The instructions for the compilation of ISSM and SLPS modules are available at http://issm.jpl.nasa.gov/download. The public svn repository for the ISSM code can also be found directly at https://issm.ess.uci.edu/svn/issm/issm/trunk, 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.