The Framework for Assessing Changes To Sea-level (FACTS) v1.0: a platform for characterizing parametric and structural uncertainty in future global, relative, and extreme sea-level change

Future sea-level rise projections are characterized by both quantifiable uncertainty and unquantifiable structural uncertainty. Thorough scientific assessment of sea-level rise projections requires analysis of both dimensions of uncertainty. Probabilistic sea-level rise projections evaluate the quantifiable dimension of uncertainty; comparison of alternative probabilistic methods provides an indication of structural uncertainty. Here we describe the Framework for Assessing Changes To Sea-level (FACTS), a modular platform for characterizing different probability distributions for the drivers of sea-level change and their consequences for global mean, regional, and extreme sea-level change. We demonstrate its application by generating seven alternative probability distributions under multiple emissions scenarios for both future global mean sea-level change and future relative and extreme sea-level change at New York City. These distributions, closely aligned with those presented in the Intergovernmental Panel on Climate Change Sixth Assessment Report, emphasize the role of the Antarctic and Greenland ice sheets as drivers of structural uncertainty in sea-level change projections.

1 Introduction

Quantitative projections of future sea-level change have been of interest to both scientists and decision-makers since at least the 1980s (Garner et al., 2018; Hall et al., 2019; Horton et al., 2018; Kopp et al., 2019). To our knowledge, the first peer-reviewed scientific article projecting 21st century global mean sea-level rise appeared in Science in 1982 (Gornitz et al., 1982). The US Army Corps of Engineers and the Dutch Ministerie van Verkeer en Waterstaat first employed planning-oriented sea-level scenarios just a few years later (National Research Council, 1987; van der Kley, 1987; US Army Corps of Engineers, 1989). Thus, sea-level projections have always been one of the more practically relevant parts of scientific assessments of climate change, including all six of the Intergovernmental Panel on Climate Change (IPCC) Working Group 1 reports (Kopp et al., 2023a).

At the same time, scientific projections of sea-level rise have also long acknowledged the presence of factors – particularly associated with Antarctic Ice Sheet instability – that limit the ability to generate quantitative sea-level projections (e.g., Mercer, 1978; Gornitz et al., 1982). These limits give rise to what is sometimes called ambiguity or deep uncertainty – uncertainty that cannot be represented by singular probability distributions due to limited amount, reliability, and unanimity of information (ambiguity as defined by Ellsberg, 1961) or, similarly, to ignorance or disagreement among analysts (subtypes of deep uncertainty as defined by Lempert et al., 2003). The question of how to integrate such ambiguity into the assessment and communications of sea-level projections has long challenged the authors of scientific assessments (Oppenheimer et al., 2019b; Kopp et al., 2023a).

Until about 15 years ago, comprehensive, localized projections of relative sea-level (RSL) change were uncommon (e.g., Katsman et al., 2008; National Research Council, 2012).¹ Many users simply augmented global mean sea-level (GMSL) projections with estimates of vertical land motion (VLM) to project local RSL change. The IPCC Fourth Assessment Report (AR4) considered only deviations from GMSL driven by ocean dynamic sea-level change as represented in coupled atmosphere–ocean general circulation models (Meehl et al., 2007). Other researchers focused on the contemporary gravitational, rotational, and deformational (GRD) RSL changes caused by redistribution of mass within the cryosphere and hydrosphere (e.g., by melting land ice) (e.g., Mitrovica et al., 2009). These two threads began to come together in the literature leading up to the IPCC Fifth Assessment Report (AR5) (e.g., Kopp et al., 2010; Slangen et al., 2012). AR5 was the first IPCC report to consider both sterodynamic sea-level change and contemporary GRD, along with the effects of glacial isostatic adjustment (GIA), in its RSL projections (Church et al., 2013a).

The AR5 projections and numerous subsequent studies taking on the challenge of producing comprehensive, localized RSL projections (see Sect. 9.6.3.1 of  Fox-Kemper et al., 2021a, for an overview) are generally referred to as “probabilistic” projections in that, under different emissions scenarios, they estimate probability distributions for the change in each of the driving factors of GMSL and RSL change and their total. Producing such projections requires combining different lines of information: global climate models (GCMs) can simulate sterodynamic sea-level change but do not generally include coupled glaciers, ice sheets, or anthropogenic changes in land water storage. They also require using relatively simple representations of sea-level drivers; models of the complexity of GCMs do not lend themselves to the Monte Carlo sampling used to estimate sea-level distributions in probabilistic sea-level projections. Examples of open-source probabilistic sea-level projection frameworks include the ProjectSL/LocalizeSL framework (Kopp and Rasmussen, 2021), developed by Kopp et al. (2014, 2017), and BRICK (Wong et al., 2017). Additional studies present probabilistic RSL projection methodologies without associated open-source software releases (e.g., Slangen et al., 2014; Grinsted et al., 2015; Jackson and Jevrejeva, 2016; Jevrejeva et al., 2019; Le Cozannet et al., 2019; Palmer et al., 2020).

Probabilistic sea-level projection frameworks are limited in that they assume that future changes under a single emissions scenario can be represented by a single probability distribution. By definition, this assumption is not true for processes characterized by ambiguity (Kopp et al., 2023a; Hinkel et al., 2019). While some studies (e.g., Kopp et al., 2017; Jevrejeva et al., 2019) have worked around this problem to explore structural uncertainties by substituting different modeling approaches for different sea-level components, probabilistic projection frameworks have not generally been engineered to facilitate such explorations.

This paper describes the Framework for Assessing Changes To Sea-level (FACTS), a scalable, modular, open-source framework for global mean, local, and extreme sea-level projection that is designed to support the characterization of ambiguity in sea-level projections. FACTS is built using modern computational practices and in the spirit of open science (e.g., Wilkinson et al., 2016). It is designed so users can easily explore deep uncertainty by investigating the implications for GMSL, RSL, and extreme sea level (ESL) of different choices for different processes. Its modularity allows components to be represented by either simple or complex models. Because it is built upon the RADICAL-Cybertools computing stack (Merzky et al., 2021), different modules can in principle be dispatched for execution on resources appropriate to their computational complexity.

FACTS is, specifically, a tool for sea-level assessment. It is not intended as a substitute for detailed process-based analyses of individual sea-level contributions (for example, GCM studies of ocean dynamics or ice sheet modeling studies) or of integrated projections made with high-complexity Earth system models that are moving toward including coupled ice sheets (e.g., Muntjewerf et al., 2021; Smith et al., 2021). Such studies provide the scientific bases underlying FACTS modules. Rather, FACTS is intended to support scientists – like those participating in the IPCC and in numerous national and subnational assessment processes – who seek to develop projections that are internally consistent, represent the richness of approaches present in the scientific literature, and assess multiple types of uncertainty. Such assessment outputs, rather than individual projections in the primary scientific literature, are generally the primary way in which climate risk practitioners interact with estimates of future sea-level change (Kopp et al., 2023a).

Development versions of FACTS modules underlie the GMSL and RSL projections of the IPCC Sixth Assessment Report (AR6) (Fox-Kemper et al., 2021a; Slangen et al., 2023) and the 2022 US Government sea-level rise technical report (Sweet et al., 2022). For these implementations, several key steps were run offline, and modules were invoked outside the execution and data management framework provided by the FACTS Manager. FACTS 1.0 allows replication of the AR6 approach entirely within FACTS, starting from specification of emissions scenarios and ending with the production of multiple alternative probability distributions for GMSL, RSL, and ESL.

Box 1 Key sea-level terminology.

2 Model description

2.1 Overview

FACTS consists of the FACTS Manager, which oversees the execution of FACTS experiments, and an extendable suite of modules, which provide the scientific and analytical core that allows FACTS to simulate the different process contributing to GMSL, RSL, and ESL change. Modules represent independent processes (e.g., sterodynamic sea-level change or VLM) and can be run in parallel on high-performance computing (HPC) resources. Modules can also be run in sequence when their outputs depend upon inputs from other modules (e.g., the modules that compute total RSL change and ESL distribution shifts).

A FACTS experiment consists of a series of experiment steps (Fig. 1). Typical experiment steps include (1) a climate experiment step, which translates an inputted emissions scenario into projections of global mean surface air temperature (GSAT) and ocean heat content change; (2) a sea-level component experiment step, which simulates the different physical processes driving sea-level change; (3) an integration experiment step, which adds up the different components into projections of total GMSL and RSL change; and (4) an ESL experiment step, which uses tide gauge data and RSL projections to project the change in extreme sea-level occurrences over time.

Each experiment step runs one or more modules in parallel. Exchange of information between modules happens between experiment steps. This exchange is mediated by the file system, so experiment steps can be bypassed simply by providing appropriate input files (e.g., stored GSAT and ocean heat content trajectories) to the subsequent experiment step. Though the existing usage of FACTS contains only one sea-level component experiment step, and therefore treats the output of each module as independent conditional upon their common dependence on the climate simulated in the climate experiment step, the FACTS Manager allows experiment steps to be subdivided and could thus support between-module coupling.

The core concept of workflow provides FACTS with the flexibility required to explore structural uncertainty. A workflow consists of a set of sea-level component modules that are added together in the integration experiment step to produce a probabilistic estimate of their combined contribution to sea-level change. Workflows can be overlapping: for example, two workflows might use the same module for simulating sterodynamic sea-level change but use different modules for simulating ice sheet change. Modules run in the sea-level component experiment step are tagged as belonging to one or more workflows; those workflows are then aggregated at the integration experiment step. This structure allows a single sea-level components experiment step to include multiple modules representing alternative methods to simulate the same sea-level component and avoids redundant execution of modules employed in multiple workflows.

In practice, for a specific set of climate inputs (e.g., emissions scenario-forced GSAT projections), a single workflow produces a single (climate input-conditional) probabilistic projection of sea-level change. Multiple workflows can be compared to examine the structural uncertainty of GMSL, RSL, and ESL change to the choice of component methods (i.e., the ambiguity of projections) and combined (for example, in a p-box, as discussed in Sect. 4.1) to produce summary outputs that capture ambiguity (Kopp et al., 2023a).

2.2 FACTS Manager and RADICAL-Cybertools

Though most of the FACTS modules implemented to date can be run on a desktop computer and all can run on small-scale HPC clusters, FACTS is designed to allow modules of a broad range of computational demands, including those requiring supercomputer resources. This objective is achieved by using the RADICAL-Cybertools software stack in the FACTS Manager.

RADICAL-Cybertools are software systems designed to support the execution, across computing scales, of applications comprised of multiple tasks. A task can be any executable or Python function; tasks can have a short (𝒪 (s)) or long (𝒪 (h to d)) duration and can run on single or multiple cores, nodes, and threads, either locally or remotely.

The RADICAL Ensemble-Toolkit (hereafter, EnTK) (Balasubramanian et al., 2016, 2018) is the top-level system of the middleware stack used to implement FACTS. EnTK is an ensemble execution system, implemented as a Python library, that offers components to encode and execute ensemble applications on HPC systems. EnTK uses RADICAL-Pilot (RP) (Merzky et al., 2021) to decouple the description of ensemble applications from their execution, separating three concerns: (i) specification of tasks and resource requirements, (ii) resource selection and acquisition, and (iii) management of task execution. EnTK sits between the user and the HPC system(s), abstracting resource and execution management complexities from the user.

EnTK exposes an API with three user-facing constructs: pipeline, stage, and task. Those constructs allow the user to encode an ensemble application in terms of concurrency and sequentiality of tasks. Each pipeline is a sequence of stages, and each stage is a set of tasks. Consistent with their formal definition, EnTK executes the members of a set concurrently and the members of a sequence sequentially. For example, all the stages of each pipeline execute sequentially, and all the tasks of each stage execute concurrently. In this way, EnTK describes an ensemble application in terms of the concurrency and sequentiality of tasks, without requiring the explicit specification of tasks' data or control dependencies.

In the context of the FACTS Manager, each experiment step contains a set of pipelines that are run concurrently. Each pipeline is associated with one FACTS module, and each module runs a series of sequential, single-task stages described in its configuration file. Most typically, these stages consist of (1) a pre-processing stage with a task that prepares associated data, (2) a fitting stage with a task that calibrates the module based on the data prepared by the pre-processing stage, (3) a projection stage, and (4) a post-processing stage. In the existing sea-level component modules, the projection stage generates the projection of GMSL contributions, while the post-processing stage generates the projection of RSL contributions. For example, in a module computing Greenland Ice Sheet contributions, the projection stage might project the Greenland contribution to GMSL, while the post-processing stage might incorporate the contemporary GRD effects that modulate the Greenland contribution to RSL change at specific sites. Note that alternative specifications are possible, e.g., the totaling module runs in a single stage.

Figure 1 Schematic illustration of the FACTS experiment described in this paper. Large boxes represent the four experiment steps. Smaller boxes represent different modules run in each experiment step. Circles represent workflows (WFs) generated by combining sets of sea-level component modules in the integration experiment step (and carried forward to the extreme sea-level experiment step). Grey modules are applied to and included in all workflows, while colored modules are included in some but not all workflows in different combinations. See Table 2 for details of the modules making up each workflow.

2.3 Modules

FACTS 1.0 includes a library of different modules (Table 1) that both illustrate functionality and allow simulation of projection workflows analogous to those employed in the IPCC AR6 (Fox-Kemper et al., 2021a). Each of the included modules is described below. Configuration options such as the number of samples to run, the time points at which calculations are reported, and the reference period used for output can be globally specified but are implemented on a module-by-module basis.

2.3.1 Climate module

Climate simulation is provided by the fair/temperature module. This module wraps around the FaIR v1.6.4 climate model emulator (Smith et al., 2018; Millar et al., 2017) using the AR6 calibrated and constrained parameter set (Smith, 2021). Taking an emissions scenario as an input, this module samples uncertainty in key climate model parameters (e.g., equilibrium climate sensitivity and transient climate response) and generates probability distributions of GSAT and ocean heat content (using the two-layer temperature function of Geoffroy et al., 2013). The climate simulation experiment step can also be bypassed by providing to the modules run in the sea-level component experiment step an output file containing these probability distributions. For application in the AR6, for example, the climate simulation was run offline and passed as an input to modules depending on these inputs.

2.3.2 Sea-level component modules

The bulk of the modules distributed with FACTS simulate physical processes that contribute to GMSL and/or RSL change. Consistent with IPCC AR6, existing sea-level component modules output quantities that are relative to the 19-year average of GMSL and/or RSL centered in the year 2005.

Generic module

The simplest module in FACTS 1.0 is the generic direct sampling module (facts/directsample), which simply translates an ensemble of time series samples specified in a text file into FACTS.

IPCC AR6 offline land ice modules

For the implementation of FACTS used to develop the IPCC AR6 sea-level projections, several of the modules used to simulate ice sheet and glacier contributions (ipccar6/ismipemuicesheets, ipccar6/gmipemuglaciers, ipccar6/larmipAIS, deconto21/AIS, bamber19/icesheets) were based upon variants of facts/directsample, with the sample inputs being generated through offline simulation.

In the case of ipccar6/ismipemuicesheets and ipccar6/gmipemuglaciers, climate output generated offline using FaIR by the AR6 Working Group 1 Chapter 7 authors was run offline through the emulandice emulator of Edwards et al. (2021), the output of which was then transferred to the FACTS modules as static data. Similarly, in the case of ipccar6/larmipAIS, the Chapter 7 climate output was run through the LARMIP-2 emulator of Levermann et al. (2020), then transferred to the FACTS module. (Details of both the emulandice and LARMIP-2 emulators are described below.) For replicability reasons, the original AR6 direct-sample versions of the emulandice ISMIP6 and LARMIP modules (ipccar6/ismipemuicesheets and ipccar6/larmipAIS, respectively) are retained in FACTS 1.0, though their use is deprecated.

Greenland and Antarctic Ice Sheet modules

In FACTS 1.0, the larmip and emulandice modules bring the formerly offline-coupled emulandice and LARMIP-2 emulators into FACTS. These modules are both driven by sampled projections of GSAT.

The emulandice modules are structured as wrappers around separately developed, R-language Gaussian process emulators for ISMIP6 ice sheet simulations and GlacierMIP glacier simulations. They demonstrate the ability of FACTS to incorporate independently developed models (Nowicki et al., 2016, 2020; Hock et al., 2019; Marzeion et al., 2020; Edwards et al., 2021). ISMIP6 (Ice Sheet Model Intercomparison Project for CMIP6: Nowicki et al., 2016, 2020) generated around 600 simulations from 2015–2100 from 27 modeling groups under very high (RCP8.5, SSP5-8.5) and low (RCP2.6, SSP1-2.6) emissions scenarios, systematically varying a small number of ice sheet model parameters driving the response. These simulations were used in constructing Gaussian process emulators of the Greenland and Antarctic (west, east, and peninsula) contributions to sea level as a function of GSAT and of these parameters (Edwards et al., 2021). Note that emulandice emulates sea-level contributions in each year independently: the outputs are samples drawn from independent distributions for each year. This means it does not include temporal autocorrelation in uncertainty and therefore does not emulate the rates of change between years, although they can be approximated by smoothing the annual percentiles with a temporal filter (temporal correlation emerges from the underlying simulations). Because the ISMIP6 experiments end in 2100, and Gaussian process emulation should not be used for significant extrapolation (being nonparametric), the emulandice modules cannot generate projections beyond 2100.

For the Greenland Ice Sheet, FittedISMIP/GrIS provides a parametric emulator for 21 models participating in the ISMIP6 exercise. The parametric emulator is based on fitting each model's projected sea-level contributions under different scenarios as a cubic function of GSAT and quadratic function of time. Details are provided by Fox-Kemper et al. (2021b). In contrast to emulandice, FittedISMIP/GrIS can be used to estimate rates of change.

The larmip module is an adaptation of separately developed code (Levermann et al., 2020), modified to achieve substantial speed improvements. Within the Linear Antarctic Response Model Intercomparison Project (LARMIP-2), 16 state-of-the-art ice sheet models performed experiments in which they applied a constant additional basal ice shelf melt forcing of 8 m yr⁻¹ underneath each of five distinct regions of the Antarctic coast for 200 years. The time derivative of the ice loss response from these experiments yielded a linear response function for each of the regions in each of the models. To apply these linear response functions to generate new projections, GSAT projections are scaled and time-delayed in according with the response of the CMIP6 climate models' subsurface oceanic warming to surface warming. This subsurface warming signal is then scaled with the observed sensitivities of basal melting to warming outside of the Antarctic ice shelf cavities. The resulting basal melt forcing is convolved with the linear response function to project the dynamic response of the Antarctic Ice Sheet.

Because the LARMIP-2 experiment examined only the dynamic response of the Antarctic Ice Sheet, projecting the full Antarctic response requires incorporating a separate term representing surface mass balance changes. This is done within the larmip module using the same approach as applied by the IPCC Fifth Assessment Report (Church et al., 2013b) and in the ipccar5/icesheets modules, as described below. Whereas in AR6, LARMIP-2 projections (including surface mass balance) are extrapolated beyond 2100 assuming a fixed rate of ice sheet mass loss after 2100, here we allow the rate of loss to evolve following the linear response function formulation.

The ipccar5/icesheets module implements the Greenland and Antarctic Ice Sheet projection methods used in the IPCC Fifth Assessment Report (Church et al., 2013b). Greenland surface mass balance is projected using a cubic polynomial of GSAT (Fettweis et al., 2013). The polynomial is multiplied with a lognormally distributed factor representing methodological uncertainty. Another multiplier varying randomly between 1 and 1.15 is added to account for positive elevation feedback. For Antarctic surface mass balance, accumulation is projected to increase by 5.1±1.5 % per degree Celsius warming in Antarctica, with a 1.1±0.2 ratio of warming in Antarctica to GSAT increase. The uncertainties in both of these numbers are assumed to be normally distributed, and a negative rate term that scales with accumulation is added to account for the feedback between enhanced accumulation and dynamic ice discharge. The ice dynamic contributions of Greenland and Antarctica are parameterized by quadratic functions of time, starting at either the lower or upper end of the uncertainty range of observed rates of ice loss over 2005–2010 and reaching the minimum or maximum contributions, respectively, of the ice sheets in 2100 that the Fifth Assessment Report assessed based on the available literature at that time. Samples are drawn assuming a uniform probability density between these extreme quadratic functions (Church et al., 2013b).

FACTS 1.0 also includes direct sampling modules used to incorporate ice sheet projections that include, either by structured expert judgment (bamber19/icesheets) (Bamber et al., 2019) or physical modeling (deconto21/AIS) (DeConto et al., 2021), processes such as marine ice cliff instability that are not included in most ice sheet models but that might have the potential to substantially accelerate the ice sheet contribution to sea level. Bamber et al. (2019) used formal structured expert judgment with calibrated expert responses to probabilistically evaluate Antarctic and Greenland mass loss through 2300 under 2 and 5^∘ C GSAT stabilization scenarios. In the IPCC AR6, these two GSAT scenarios were mapped to SSP1-2.6 and SSP5-8.5 projections. DeConto et al. (2021) projected future Antarctic Ice Sheet changes under the Representative Concentration Pathway (RCP) scenarios using a model that incorporates hydrofracturing of ice shelves and the gravitational instability of marine ice cliffs without the protection of a buttressing ice shelf. In the IPCC AR6, the RCP scenario projections were employed in the context of the corresponding SSP projections (e.g., RCP2.6 projections from DeConto et al., 2021, applied to SSP1-2.6).

All the existing ice sheet modules include in their post-processing stage a regional scaling based on GRD fingerprints for West Antarctica, East Antarctica, and Greenland (e.g., Mitrovica et al., 2001; Gomez et al., 2010; Mitrovica et al., 2011). The fingerprints include both gravitational and rotational effects on sea-surface height, as well as deformational effects on seafloor height. They are implemented as static fingerprints that do not change over time; as in Kopp et al. (2014), mass change is assumed to be uniform across the respective regions. The fingerprints were pre-computed (outside the FACTS framework) by solving the sea-level equation with a pseudo-spectral approach up to spherical harmonic degree and order 512 (equivalent to a spatial resolution of about 0.4^∘). They assume a radially symmetric, elastic, and compressible Earth model based on the Preliminary Reference Earth Model (Dziewonski and Anderson, 1981).

Glacier modules

The emulandice/glacier module, like the emulandice/AIS and emulandice/GrIS module, is based on Gaussian process emulation of a multi-model intercomparison exercise, specifically the GlacierMIP2 ensemble (Marzeion et al., 2020), and is driven by inputted GSAT trajectories. The GlacierMIP2 project generated nearly 300 simulations of 2015–2100 glacier loss from 11 modeling groups under four RCP scenarios. These simulations were used in constructing Gaussian process emulators of the 19 glacier region contributions to sea level as a function of GSAT (Edwards et al., 2021). Because the GlacierMIP experiments end in 2100 (as for the ice sheets), the emulandice modules cannot generate projections beyond 2100.

The ipccar5/glaciers module is based on the glacier projection approach used in the IPCC Fifth Assessment Report (Church et al., 2013b), which models the global mean sea-level change due to the melt of glaciers as f×I(t)^p, where I(t) is the time integral of GSAT at time t, and f and p are parameters estimated from simulations of a set of four glacier models (Giesen and Oerlemans, 2013; Marzeion et al., 2012; Radić et al., 2014; Slangen and Van De Wal, 2011). The glacier models are equally weighted and systematic uncertainty in the glacier projections is accounted for by Monte Carlo sampling, assuming a normal distribution with a time-dependent model-specific standard deviation. For the glacier projections of the IPCC AR6, these parameters were also derived from the simulations of GlacierMIP and GlacierMIP2 and added as calibration options to the ipccar5/glaciers module (Hock et al., 2019; Marzeion et al., 2020; Fox-Kemper et al., 2021a) (in this paper, as in IPCC AR6, we focus on the GlacierMIP2 calibration of this module, denoted as GMIP2 in tables). As the IPCC AR5 model itself does not disaggregate the glacier contribution into separate regions, this disaggregation is based upon the time-varying proportion of the contributions of different glaciers in the median projection of Kopp et al. (2014).

As described in Kopp et al. (2014), the kopp14/glaciers module projects the contribution of 17 different glaciers and ice cap regions for different RCPs by employing a multivariate t distribution of ice mass change estimated from the model simulations of Marzeion et al. (2012) for different source regions.

As with the ice sheet modules, the glacier modules scale their output in the post-processing stage using offline-calculated fingerprints. As in Kopp et al. (2014), the lookup library includes separate GRD fingerprints for 17 different glacier regions, and thus the spatial pattern associated with glaciers as a whole can change over time in response to the spatial distribution of glacier mass loss.

Sterodynamic modules

Several modules are included to project sterodynamic sea-level change, i.e., the sum of global mean thermosteric sea-level rise and ocean dynamic sea-level change (Box 1). As described in Fox-Kemper et al. (2021a), the tlm/sterodynamics module does so by taking as input the emulated ocean heat content from the fair/temperature module and pre-processed gridded simulations of CMIP6 models (as noted above, fair/temperature is run using a two-layer model representation of the forcing–temperature coupling, from whence comes the abbreviation “tlm”). Global mean thermosteric sea-level rise is projected by sampling from a distribution of time-invariant global thermal expansion coefficients derived from CMIP6 simulations (Fox-Kemper et al., 2021a) and multiplying the emulated ocean heat content by those coefficients. The CMIP6 simulations that were used for the calibration of the expansion coefficients are shown in Table A3. The resulting global mean thermosteric sea-level rise is then combined with ocean dynamic sea-level change and the inverse barometer (IB) effect using the gridded output of CMIP6 models (see the right column of Table A3 for the models that were used in Fox-Kemper et al., 2021a) based on the time-varying correlation structure between global mean thermosteric sea-level rise and ocean dynamic sea-level change in the multi-model ensemble. The tlm/sterodynamics module expects the CMIP6 input to be pre-processed (e.g., dedrifted and regridded) a priori. The approach used in the provided data set is described in Fox-Kemper et al. (2021b), and further details are provided in Appendix A.

The sterodynamic component is also provided by the kopp14/sterodynamics module, which implements the methodology of Kopp et al. (2014). In this module, drift-corrected global mean thermosteric sea-level rise is characterized for specific Representative Concentration Pathway scenarios using a t distribution with the mean and covariance derived from the CMIP5 multi-model ensemble. As in tlm/sterodynamics, ocean dynamic sea-level change is then projected using the time-varying correlation structure between global mean thermosteric sea-level change and ocean dynamic sea-level change in the multi-model ensemble.

As described in Church et al. (2013b), the ipccar5/thermalexpansion module projects the distribution of global mean thermosteric sea-level rise. It is calibrated to the time-dependent mean and standard deviation of the global mean thermosteric sea-level rise simulated by a multi-model ensemble. Samples are drawn from the mean and standard deviation assuming a normal distribution. The same method was applied by several studies and reports published between the Fifth and Sixth Assessment Reports of the IPCC (Palmer et al., 2018, 2020; Hermans et al., 2021).

Land water storage modules

Two modules provide the land water storage component of sea-level change. As described in Kopp et al. (2014), the first, kopp14/landwaterstorage, estimates this component based on the relationship between changes in land water storage and global population change using United Nations population projections. Reservoir storage is assumed to follow a sigmoidal function of population change, calibrated based on Chao et al. (2008). The relationship between groundwater depletion and population change is based on linear fits to estimates of Wada et al. (2012) and Konikow (2011). The groundwater projection of Pokhrel et al. (2012), based upon a water resource assessment model, is included as an option for sensitivity analysis. Uncertainty in the projections is generated by sampling the parameters of the sigmoidal fit for reservoir storage and linear fit for groundwater depletion.

The second module, ssp/landwaterstorage, follows the methods of Kopp et al. (2014), except for three aspects: (1) instead of using scenario-independent global population projections, population projections of the different SSPs were used (Samir and Lutz, 2017); (2) the groundwater depletion component was multiplied by 0.8 to account for only 80 % of depleted groundwater reaching the ocean (Wada et al., 2016); and (3) the capability to add a temporally linear adjustment for projected reservoir storage based on planned dam construction was added (and applied in AR6 projections using the Hawley et al., 2020, projections for 2020–2040). The GRD fingerprint used is based on the groundwater source pattern of Wada et al. (2012), as described in Slangen et al. (2014).

Long-term vertical land motion and glacio-isostatic adjustment modules

Long-term VLM (as well as the sea-surface height contribution from GIA) is provided by the kopp14/verticalandmotion module. As described in Kopp et al. (2014), this module estimates a constant trend at each spatial location based upon a Gaussian process spatiotemporal analysis of tide gauges in the Permanent Service for Mean Sea Level database. This Gaussian process analysis modifies an estimated rate of long-term trend derived from a single GIA model (ICE5G-VM2-90), treating this GIA model as a prior mean rate estimate whose misfits are statistically corrected. Sensitivity tests show this approach exhibits little sensitivity to the choice of initial GIA model in the vicinity of tide gauge records; more substantial differences can occur in parts of the polar region that do not have good observational constraints (Fig. A1).

The spatiotemporal model assumes observed RSL can be described as the sum of a uniform (and independently estimated) global component, a regionally varying autocorrelated nonlinear component (with a decorrelation timescale of the order of 1–3 years), and a regionally varying constant trend. The spatial and temporal correlation scales of the regional components are separately tuned (via maximum likelihood optimization) along different coastal segments. The constant trend is assumed to equal the long-term contribution from VLM (including the VLM term arising from GIA), as well as from the sea-surface height trend arising from GIA, and is propagated into the projection. Uncertainty in the projection is generated based on the uncertainty in the estimate of the constant trend.

Because the statistical model is constructed to extract a century-scale, climate-uncorrelated trend, there should be minimal double-counting of the deformational effects associated with recent land ice mass loss and land water redistribution. This may be a concern along coastlines with only short tide gauge records, but the resulting bias remains small because future projected rates of land ice changes are substantially larger than the average rates over the last several decades. VLM associated with future land ice mass loss and land water redistribution is incorporated into the GRD projections of those components' respective modules.

An alternative direct-sampling-based VLM approach is demonstrated by the NZInsarGPS/verticallandmotion module, which reads and samples gridded land motion data described in an external file and extrapolates these rates linearly into the future. In Naish et al. (2023), this module applies a gridded data file describing rates of land motion inferred from interferometric synthetic aperture radar (InSAR) data.

2.3.3 Totaling module

The facts/total module handles the aggregation of sea-level component probability distributions into probability distributions for total GMSL and RSL change. This module takes as an input a configuration file pointing to the output files that constitute different workflows (see Sect. 2.4).

2.3.4 ESL module

The extremesealevel/pointsoverthreshold module, which is based on the methods of Oppenheimer et al. (2019a) and Frederikse et al. (2020), first derives declustered ESLs from tide gauge data from the GESLA2 database (Woodworth et al., 2016) using a peak-over-threshold method with a user-defined threshold percentile. After removal of the annual means, a generalized Pareto distribution is fitted to the declustered extremes using maximum likelihood estimation. The estimated parameters and their uncertainty are used to generate ESL return-period curves. Below the threshold of the generalized Pareto distribution, a Gumbel distribution with support between mean higher high water and the threshold is assumed and used to compute return periods, following Buchanan et al. (2016). In the projection stage, the module augments the return-period curves by projected RSL change to project how the expected frequency of ESL events of different magnitudes changes as the baseline height of the events is increased (Frederikse et al., 2020). Note that this approach assumes that the ESL distribution, relative to a changing mean sea level, is stationary; it does not account for factors such as changes to storm frequency, intensity, or tracks.

Table 1 Modules included in FACTS 1.0.

The ipccar6 modules are direct-sample modules that were used only in IPCC AR6 and have been deprecated in FACTS 1.0 in favor of the emulandice and larmip modules. The ipccar5 modules indicate the methods described in Church et al. (2013b), which in some cases and contexts were used by AR6, as described in Fox-Kemper et al. (2021a) and Table 2. The ipccar5/glaciers module includes, in addition to the original IPCC Fifth Assessment Report calibration, recalibrations to GlacierMIP and GlacierMIP2 (Hock et al., 2019; Marzeion et al., 2020). The GlacierMIP2 recalibration is used in IPCC AR6 and in this paper and is denoted by a parenthetical “(GMIP2)” in Tables 2 and 3.

2.4 Workflows

In this paper, we demonstrate the FACTS capabilities by implementing seven different workflows (i.e., sets of sea-level component modules) (Table 2). The workflows align with those implemented by IPCC AR6. As previously described in the description of the IPCC AR6 land ice modules, in FACTS 1.0, we replace the direct sampling of offline calculated values used in AR6 with GSAT-driven emulandice and larmip modules. The workflows share a common set of modules used for projecting VLM (kopp14/verticallandmotion), sterodynamic sea level (tlm/sterodynamics), and land water storage (ssp/landwaterstorage). They differ based on their handling of the cryospheric components (ice sheets and glaciers).

Workflows 1e and 2e employ Gaussian process emulation of ice sheet and glacier intercomparison exercise outputs for Greenland, glaciers, and, in the case of workflow 1e, Antarctica (i.e., emulandice in Table 2). However, the Gaussian process emulator of Edwards et al. (2021) models each time point independently and thus does not estimate rates. Because emulandice uses a nonparametric (Gaussian process) model, where no functional form is assumed, rather than a parametric model, in which dependencies are asserted, workflows using emulandice modules can only project up to the end of the original simulations (rather than extrapolate beyond them) and therefore end in 2100. Workflows 1f and 2f therefore substitute alternative parametric representations for GrIS and glaciers. Workflows 2e and 2f differ from workflows 1e and 1f by employing an alternative Antarctic Ice Sheet emulator, provided by the larmip module. These four workflows together form the basis of the medium-confidence projections presented by AR6 (Fox-Kemper et al., 2021a) (for example, in the unshaded columns of Table 9.9 of Fox-Kemper et al., 2021a). Workflows 3e, 3f, and 4, by contrast, take alternative approaches to ice sheet representation intended to capture processes not represented in most ice sheet models. Workflows 3e and 3f employ the deconto21 projections for Antarctica, while workflow 4 employs structured expert-judgment-based projections (bamber19) for both the Greenland and Antarctic ice sheets. These three workflows are combined with the medium-confidence workflows to form the basis of the broader AR6 low-confidence projections (for example, for SSP5-8.5, in the final column of Table 9.9 of Fox-Kemper et al., 2021a).

Table 2 Workflows used in this paper.

Workflows used in this paper match those of AR6 (Fox-Kemper et al., 2021a), except that in AR6, results from ISMIP/GlacierMIP emulation and LARMIP-2 were computed offline by those models' authors and then added into the projections as static data rather than online in a FACTS experiment as done by the emulandice and larmip modules.

Table 3 Component projections for 2100.

Median (17th–83rd percentile) projections produced by FACTS modules. All components except GSAT are in meters of GMSL contribution relative to a 1995–2014 baseline. Global mean surface air temperature (GSAT) is in degrees Celsius (^∘C) relative to a 1850–1900 baseline. For certain modules (marked with an asterisk), projections for Representative Concentration Pathways 2.6, 4.5, and 8.5 are shown in lieu of those for the SSP scenarios.

3 Results

All results presented are based on 2000 pseudo-random Monte Carlo samples. To illustrate the application of FACTS, we focus on GMSL projections and on RSL and ESL projections at a single site, New York City.

3.1 Temperature projections

FACTS experiments begin with the estimation of the GSAT response to emissions forcing, as projected by the FaIR climate emulator. By construction, these projections are generally consistent with those of AR6 (Lee et al., 2021), with median warming in 2100 above 1850–1900 of 1.6 ^∘C in SSP1-2.6, 2.6 ^∘C in SSP2-4.5, and 4.7 ^∘C in SSP5-8.5 (Table 3). Note that SSP1-2.6 is aligned with the Paris Agreement goal of limiting warming to well below 2 ^∘C, while SSP2-4.5 is closer to projected emissions under current policy. SSP5-8.5 emissions represent a high-end trajectory that would require a reversion to fossil-fuel-intensive development (Riahi et al., 2022).

Figure 2 GMSL contributions from the Greenland Ice Sheet (GrIS), the Antarctic Ice Sheet (AIS), glaciers, and thermal expansion (TE) for SSP1-2.6 and SSP5-8.5 based upon different FACTS modules. Curves show median projections. Thick and thin bars on the right show the 17th–83rd and 5th–95th percentile projections for 2100.

3.2 Global mean contributions from sea-level components

In the sea-level component experiment step, FACTS estimates the contributions to future GMSL and RSL rise from the cryosphere, land water storage, and sterodynamic sea-level change. Some sea-level component modules (for example, the sterodynamic, ice sheet, and glacier modules used in workflows 1e, 1f, 2e, and 2f) take the FaIR-projected warming as an input. Others rely upon pre-computed projections, in some cases indexed by an SSP or RCP emissions scenario (for example, the deconto21 and bamber19 ice sheet modules and the deprecated ipccar6 ice sheet and glacier modules) (Table 1).

Projected median and 17th–83rd percentile GMSL contributions are shown in Table 3 and Fig. 2. The cryosphere as a whole (including glaciers and polar ice sheets) dominates median projections for 2100 under all emissions scenarios, but the relative contribution of polar ice sheets in particular varies substantially across modules. This is particularly the case under very high emissions (SSP5-8.5, RCP8.5), where one module (deconto21) projects the Antarctic contribution to be the single largest term. For the polar ice sheet contributions, GMSL contributions projected by different modules are similar until 2040 but begin to diverge beyond 2050, and this divergence is larger for higher emission scenarios (Fig. 2). By contrast, both glacier modules (ipccar5 and emulandice) remain consistent throughout this century; this is to be expected, given that both are calibrated to the same underlying GlacierMIP ensemble of glacier model projections (Marzeion et al., 2020).

3.3 Total global mean sea-level change projections

Total GMSL projections (Table 4, Fig. 3) are generally in close agreement between workflows using the emulandice emulators of ISMIP6 and GlacierMIP projections (i.e., workflows 1e, 2e, and 3e) and the corresponding workflows that substitute parametric emulators (i.e., workflows 1f, 2f, and 3f) (Table 4). For example, under SSP2-4.5, total projections are 0.50 m (0.42–0.60 m, 17th–83rd percentile range) under workflow 1e and 0.49 m (0.40–0.59 m) under workflow 1f; differences are smaller for lower emissions scenarios and for other emulandice–parametric workflow pairs (i.e., 2e vs. 2f and 3e vs. 3f) and larger for higher emissions. For the remainder of this paper, we therefore focus primarily on workflows 1f, 2f, 3f, and 4. The text primarily discusses SSP5-8.5, for which different workflows show the greatest distinctions. Figures highlight the difference between SSP1-2.6 and SSP5-8.5, as these are the two scenarios that can be projected using all workflows.

Figure 3 Total GMSL projections under four different workflows under different SSP scenarios. Curves show median projections. Thick and thin bars on the right show the 17th–83rd and 5th–95th percentile projections for 2150. Note change of y-axis scale between workflows 1f and 2f (top two rows), representing medium-confidence processes, and workflows 3e and 4, which include low-confidence processes.

Substantial differences arise between workflows based particularly on the choice of Antarctic module. Under SSP5-8.5, median projections for 2100 differ by 0.14 m between workflow 1f (Antarctica calibrated as per IPCC AR5: 0.66 (0.55–0.78) m) and workflow 2f (Antarctica calibrated to LARMIP2: 0.80 (0.60–1.00) m), with the latter projections also exhibiting fatter tails. This reflects the differences seen at the component level (Table 3, Fig. 2). Larger differences are seen under higher emissions scenarios with the two workflows that AR6 employed to incorporate low-confidence processes. Both workflow 3f (Antarctic Ice Sheet modeling marine ice cliff instability: 0.97 (0.81–1.17) m) and workflow 4 (both Antarctica and Greenland based on structured expert judgment: 1.00 (0.69–1.64) m) have median projections for SSP5-8.5 exceeding those of the medium-confidence workflows by at least 0.17–0.20 m. The median projections for workflow 3e and 4 are closely aligned, but the structured expert-judgment-based projections (workflow 4) span a larger range, primarily reflecting greater Greenland Ice Sheet uncertainty than in workflow 3f.

Consistent with these observations, by the end of the century, total projection variance is generally dominated by polar ice sheet uncertainty, particularly under workflows 2f, 3f, and 4 (Fig. 4). In addition, workflows 1f and 3f reveal a positive interaction term: i.e., the variance of GMSL projections is greater than the sum of the variances of the individual components. This positive interaction term arises because global mean thermosteric sea-level rise, glacier loss, and (in the medium-confidence workflows) polar ice sheet loss share a common dependence on GSAT and are thus positively correlated.

Table 4 Total projections for 2100.

Median (17th–83rd percentile) projections are shown relative to a 1995–2014 baseline.

Table 5 Frequency amplification factors in the years 2050 and 2100 for the historic 1 % average annual probability (100-year return period) extreme sea-level event at New York City.

Median (17th–83rd percentile) projections are shown as a ratio of event probability in 2050 or 2100 to event probability in 1995–2014.

Figure 4 Variance decomposition of GMSL change through 2100 under SSP1-2.6 (left column) and SSP5-8.5 (right column) under workflows 1f, 2f, 3f, and 4 (top to bottom) in the style of Hawkins and Sutton (2009). Each colored wedge represents the variance across Monte Carlo samples for a particular component under the specified scenario and workflow normalized by the variance of projections for total sea-level change in the same scenario and workflow. The difference between the sum of component variances and the total variance (normalized to 1.0) represents the interaction among components.

Figure 5 Projected New York City RSL contributions from the Greenland Ice Sheet, the Antarctic Ice Sheet, glaciers, and sterodynamic sea level from different FACTS modules under SSP1-2.6 and SSP5-8.5. Curves show median projections. Thick and thin bars on the right show the 17th–83rd and 5th–95th percentile projections for 2100.

3.4 Relative and extreme sea-level projections for New York City

The differences between projected GMSL rise and projected RSL rise at New York City are consistent with past studies (e.g., Kopp et al., 2014) (Table 4; Figs. 5, 6, 7). The median contribution and variance arising from the distant Antarctic are increased due to GRD effects, which cause West Antarctic ice sheet loss to cause about 20 % greater sea-level rise at New York City than in the global mean, while the median contributions and variance arising from the Greenland Ice Sheet and global glaciers are reduced due to relative proximity. The median sterodynamic contribution and its variance are larger than global mean thermosteric sea-level rise due to the potential contribution from a slowdown of the Atlantic Meridional Overturning Circulation (Yin et al., 2009; Yin and Goddard, 2013). A long-term GIA trend, arising primarily from land subsidence, adds a steady 1.5±0.2 mm yr⁻¹ to RSL rise, shifting all projections upward but contributing little to variance.

As with GMSL projections, substantial differences arise between workflows based particularly on the choice of Antarctic module. Under SSP5-8.5, workflow 1f (0.90 (0.71–1.10) m) and workflow 2f (1.07 (0.86–1.34) m) differ by 0.17 m in the median. Further differences are seen with the two workflows that AR6 employed to incorporate low-confidence processes. Notably, because high-end GMSL projections in workflow 4 include a larger Greenland contribution than in other workflows, and because Greenland's effects on RSL rise at New York City are less than in GMSL rise, median workflow 4 RSL projections (1.22 (0.88–1.73) m) are lower than workflow 3e (1.27 (1.04–1.51) m), which relies more heavily on Antarctica to drive high-end GMSL rise. While the workflow 4 tail remains the fattest of all workflows, workflow 3e's tail is fattened substantially compared to GMSL because of the heightened response of New York City RSL to Antarctic mass loss.

Figure 6 Total RSL projections for New York City under four different workflows under different SSP scenarios. Curves show median projections. Thick and thin bars on the right show the 17th–83rd and 5th–95th percentile projections for 2150.

Figure 7 Variance decomposition of New York City RSL change through 2100 under SSP1-2.6 (left column) and SSP5-8.5 (right column) under workflows 1f, 2f, 3f, and 4 (top to bottom) in the style of Hawkins and Sutton (2009). Each colored wedge represents the variance across Monte Carlo samples for a particular component under the specified scenario and workflow, normalized by the variance of projections for total sea-level change in the same scenario and workflow. The difference between the sum of component variances and the total variance (normalized to 1.0) represents the interaction among components.

Differences in RSL projections translate into differences in ESL projections (Table 5, Fig. 8). For example, under workflow 1f, the historic 1 % average annual probability extreme sea level at New York City (estimated at 1.83 m above mean higher high water) is projected to occur 2.6 (1.8–4.2) times more often by 2050 and 6.5 (3.2–17.2) times more often by 2100 under SSP1-2.6 due to the effects of RSL rise, as well as 2.8 (1.9–4.7) times more often by 2050 and 22.1 (8.0–90.3) times more often by 2100 under SSP5-8.5. The 83rd percentile projected amplification factors are all <6.8 by 2100 under SSP1-2.6, but under SSP5-8.5 and workflow 4, the 83rd percentile SSP5-8.5 amplification factor exceeds 10 000 – meaning that, under the high end of this fat-tailed projection, the historic 100-year ESL event might occur over 100 times per year. (Note that ESL return periods do not translate directly into flooding or flood damage; see Rasmussen et al., 2022, for a critique of ESL amplification factors as a metric and Hermans et al., 2023, for presentation of a related approach.)

Figure 8 Extreme sea-level return-period curves at New York City under SSP1-2.6 and SSP5-8.5 under workflows 1f, 2f, 3f, and 4 in the historical period (black) as well as the years 2050 (blue) and 2100 (red). The solid line shows the median projection; shading shows 17th–83rd percentile estimates. The dashed vertical line indicates the extreme sea-level height associated with the historic 1 % average annual probability event.

4 Discussion

4.1 Applications to date

The modular approach adopted by FACTS intentionally lends itself to careful consideration of both parametric and structural uncertainty in sea-level projections. Indeed, FACTS modules have already been used to support several major assessments of sea-level change. As previously noted, the IPCC AR6 sea-level projections were developed using FACTS modules, and the example workflows described in this paper replicate the AR6 analysis within rounding errors (Table A1). Slightly larger discrepancies with total projections (Table A2) are attributable to the combination of rounding errors and differences in sampling. (Note that AR6 used 20 000 samples per workflow compared to the 2000 per workflow in the results shown here.)

AR6 followed the development of workflow probability distributions with a particular approach to combine alternative probability distributions based upon probability boxes, or p-boxes (Kriegler and Held, 2005; Le Cozannet et al., 2017). The p-boxes describing a set of probability distributions encompass the cumulative distribution functions of the underlying probabilities; for example, the outer 17th–83rd percentile range of a p-box spans from the lowest 17th percentile of all distributions considered to the highest 83rd-percentile. All the distributions considered by construction agree that there is at least a 66 % chance that the true value falls within this particular range. Fox-Kemper et al. (2021a) used outer 17th–83rd percentile p-box ranges to characterize its likely ranges, where likely in IPCC terminology means a 66 %–100 % chance (Mastrandrea et al., 2010). (This is a difference from the definition of likely range used in the rest of the IPCC AR6 Working Group 1 report, which specifically refers to the 17th–83rd percentile of a single estimated probability distribution.) Workflows employing ISMIP and GlacierMIP emulators (1e, 2e, and 3e) were preferred over those with simple parametric representations for land ice where possible, but workflows employing these simple representations (1f, 2f, and 3f) were used when required for rates which were not emulated. Workflows 1e–1f and 2e–2f were combined in a p-box to produce the AR6 medium-confidence projections, while workflows 3e–3f and 4 were added for low-confidence projections.

The US Interagency Task Force on Sea-level Rise Scenarios (Sweet et al., 2022) built upon the same FACTS output as AR6 but took a different approach to summarizing their results. Intending to produce a set of plausible global and regional sea-level scenarios to guide decision-making – rather than, as in AR6, to characterize the likelihood of different future outcomes – the task force filtered the samples of sea-level rise associated with workflows 1f, 2f, 3f, and 4 to identify five subsets consistent with a range of 21st century GMSL rise. This range was semi-independently defined, based in part on an interpretation of the range of values presented in the AR6, to span from as low as 0.3 m (roughly, a continuation of late 20th century GMSL range) to as high as 2.0 m, the latter informed by AR6's conclusion that low-likelihood, high-impact processes could elevate GMSL above the likely range by more than 1 m. The median of each subset forms the center of each set of GMSL and RSL scenarios, while the 17th and 83rd percentiles of each subset provide within-scenario high- and low-sensitivity cases.

Both IPCC and the US Interagency Task Force invested significant effort in communicating these projections. For example, the NASA Sea Level Change Team, in partnership with these two groups, developed interactive projection viewers (at https://sealevel.nasa.gov/ipcc, last access: 14 December 2023 and https://sealevel.nasa.gov/task-force-scenario-tool, last access: 14 December 2023) to allow practitioners to explore the projections for sites around the world and the US, respectively.

FACTS has also been used to develop national RSL projections for New Zealand (Levy et al., 2020; Naish et al., 2023). In these studies, existing workflows (either based on Kopp et al., 2014, or matching those employed in AR6) were amended by replacing the existing Kopp et al. (2014)-based projections of GIA and VLM with gridded estimates based on interferometric synthetic aperture radar (InSAR) data, calibrated with ground-based Global Navigation Satellite System (GNSS) measurements. This substitution reflects a need common to many national and subnational sea-level assessments, which seek consistency with broader assessments while substituting in information that can be assessed in greater detail at a local scale.

4.2 Directions for improvement

From a scientific perspective, a number of different directions promise improvement in FACTS projections.

At present, many but not all the modules within FACTS accept climate information as an input. In particular, the ice sheet modules used to project deeply uncertain ice sheet processes (the bamber19/icesheets and deconto21/AIS modules) rely upon direct sampling of output generated by individual studies (Bamber et al., 2019; DeConto et al., 2021). This means they can be applied only to a limited set of climate scenarios. For example, Bamber et al. (2019) produced projections for 2 and 5 ^∘C GSAT stabilization scenarios. With some caveats, these are used by AR6 to inform the projections for SSP1-2.6 and SSP5-8.5 (though SSP1-2.6 most likely stabilizes below 2 ^∘C and SSP5-8.5 continues after 2100 to warm well above 5 ^∘C). Generalizing emulation approaches to encompass these alternative sources of information would allow projections of low-confidence processes to be generated for arbitrary climate scenarios; doing this cautiously might require either advances in the primary literature or a great deal of humility and uncertainty regarding the assumptions used for scenario interpolation. (It is difficult, for example, to infer the warming level associated with critical thresholds in ice sheet behavior from only two climate scenarios.)

The existing VLM modules assume a constant-rate trend into the future. While perhaps the best assumption that can be undertaken at a global scale, more refined approaches might be possible at a local scale. For example, in many regions, VLM is driven in part by highly localized subsidence associated with anthropogenic interventions, such as fluid withdrawal and/or surface loading (Shirzaei et al., 2021). Indeed, in the regions of the world experiencing the fastest rates of RSL rise, these currently tend to be the largest drivers. In such cases, assumption of a constant-rate trend in the future may not be the most suitable assumption. A module capable of representing alternative scenarios of such factors and their evolution over time could be helpful in assessments in such regions (e.g., Minderhoud et al., 2020). Alternatively, such scenarios might be independently added to RSL projections that have the VLM component removed.

The current VLM modules also do not explicitly address uncertainty in GIA (e.g., Melini and Spada, 2019). In the kopp14/verticallandmotion module, GIA uncertainty does not make a substantial contribution in locations where tide gauges are available to constrain long-term changes; however, this uncertainty can be significant at sites distant from tide gauges, particularly in polar regions where the GIA contribution is largest (Fig. A1). New approaches to fusing model projections with geological, tide gauge, and satellite observations could better characterize this uncertainty (e.g., Caron et al., 2018).

The existing ESL module treats the shape and scale of ESL return-period curves as stationary, with the distribution only shifted vertically by the increment of RSL change. In fact, ESL return periods will change due to processes such as shifts in tropical cyclone intensity and tracks (Fox-Kemper et al., 2021a), and an enhanced ESL module could incorporate parametric representations of such changes. A further refinement could capture the relationship between these changes and GSAT, allowing the analysis to reflect correlations between RSL change and storm changes (Lockwood et al., 2022). Alternative ESL modules might also use different data sources, such as regional or global hydrodynamic models (e.g., Dullaart et al., 2020).

Contemporary GRD processes in current FACTS modules are currently based on a library of scaling factors (sometimes called “fingerprints”) applied to ice sheet, regional glacier, and land water storage projections in each module's post-processing stage. Such a library approach is most appropriate for glaciers, as the glacier regions are geographically small enough that the shifts in the locus of mass loss within a region will not substantially modify that region's fingerprint. For the larger ice sheets, however, the locus of ice mass change can significantly affect the contemporary GRD spatial pattern (Larour et al., 2017; Mitrovica et al., 2018; Cederberg et al., 2023). This variability could be incorporated into FACTS through more spatially resolved ice sheet emulation, as well as potentially through a new integration module that includes an online GRD solver (e.g., Larour et al., 2020).

The existing FACTS modules start projections in the 21st century. This choice is, in part, a limitation of the underlying studies on which these modules are built. While CMIP6 historical climate simulations start in 1850, neither the GlacierMIP nor ISMIP6 model intercomparison exercises include historical simulations (Nowicki et al., 2016, 2020; Hock et al., 2019; Marzeion et al., 2020). Implementing work-arounds to these absences – and incorporating new historical simulations as they become available – would allow FACTS projections to start in the 19th or 20th century and thus enable model–data comparison for changes in both sea level and individual components. This, in turn, could allow the probability distributions generated by FACTS to be updated in a Bayesian sense based on historical, current, and (as time proceeds) future observations. While current observations are unlikely to significantly reduce the deep uncertainty associated with late 21st century high-end projections (Kopp et al., 2017), they could have a substantial impact on nearer-term projections. A Bayesian approach could also be coupled to economic models (e.g., Depsky et al., 2023) to assess the value of information associated with additional observational constraints and process model enhancements.

Relevant to such a model–data comparison is that the existing FACTS sea-level component modules are focused on projecting changes in tidal-datum-epoch (i.e., 19-year) average mean sea level, not higher-temporal-frequency (e.g., interannual) variability. Alternative ocean dynamics modules could introduce this higher-frequency variability, or, alternatively, the autocorrelation structure of such variability could be incorporated into model–data comparisons.

Some modeling approaches may require more communication between modules. At the moment, all sea-level components are computed independently, conditional upon a common input of projected GSAT and/or ocean heat content in the climate step. As a consequence, uncertainty within individual workflows is probably underestimated (e.g., Le Bars, 2018; van de Wal et al., 2022). While correlations between GSAT and global ocean heat content change will tend to lead to some correlation between projected sea-level components (e.g., Palmer et al., 2020), correlations associated with regional or systematic changes are not represented. This is not in the case in Earth system models, where (for example) meltwater input affects sterodynamic sea level (e.g., Lambert et al., 2021). Representations of such interactions could be incorporated into FACTS by subdividing the sea-level component experiment steps, either recursively refining projections with one-way coupling (e.g., modifying an initial dynamic sea-level projection for meltwater input) or proceeding in incremental time steps with two-way coupling.

More broadly, to date, FACTS has been developed by a small team, with a primary objective being to support specific assessment processes, particularly that of the IPCC AR6. A critical objective moving forward is to transform FACTS into a larger-scale community project, with modules developed autonomously by different research and assessment teams. The structure of FACTS – which enables modules to serve as wrappers around independently developed code – is intended to facilitate such efforts.

5 Conclusions

Sea-level rise is a major driver of climate risk to coastal communities and ecosystems around the world. Appropriately managing this risk requires planners to be cognizant of both quantifiable and structural uncertainty in projections of future sea-level change, and synthesizing this information is an important task of scientific assessment processes. FACTS provides a flexible, modular, and open-source platform that allows comparable probabilistic outputs to be generated in parallel through multiple modeling approaches. Its flexibility allows it to be customized based on the needs of specific assessment processes (e.g., substituting alternative approaches to VLM or higher-resolution sterodynamic sea level), while its parallel workflow structure supports the characterization of deep uncertainty.

Appendix A: The tlm/sterodynamics methodology

The ocean dynamic sea-level projection method used by the tlm/sterodynamics module is a modification of that described in Kopp et al. (2014). Whereas in Kopp et al. (2014) global mean thermosteric sea-level rise projections are derived directly from a GCM ensemble, in tlm/sterodynamics they are generated from the two-layer model, as described in Fox-Kemper et al. (2021b).

As in Kopp et al. (2014), ocean dynamic sea level is assumed to have a degree of correlation with global mean thermosteric sea-level rise, with the correlation assessed on a grid cell basis. In the case of tlm/sterodynamics, the correlation is calculated based on the CMIP6 ensemble for a particular (specified) SSP scenario. Given a sample of 19-year-average global mean thermosteric sea-level rise y at a particular point in time t, 19-year-average ocean dynamic sea level z is taken as distributed following a t distribution with a conditional mean of

$$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{‾}}{z_t}\left(r\right)+{\mathit{\sigma }}_t\left(r\right)k_t\left(r\right){\displaystyle \frac{y_t-\stackrel{\mathrm{‾}}{y_t}}{s_t}}\end{array}$$

and a conditional standard deviation proportional to

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_t\left(r\right)\mathrm{1}-k_t(r{)}^{\mathrm{2}},\end{array}$$

where z_t(r) is the multi-model mean ocean dynamic sea level at time t and location r, σ_t(r) is the multi-model standard deviation, k_t(r) is the correlation between global mean thermosteric sea-level rise and z_t(r), $\overline{y_t}$ is the multi-model mean of global mean thermosteric sea-level rise, and s_t is the standard deviation across models of global mean thermosteric sea-level rise. The standard deviation is inflated relative to that of the ensemble to account for the expert judgment that the 5th–95th percentile of the ensemble may have as much as a 33 % chance of being exceeded on either end (i.e., the 5th–95th percentile range is treated as a likely range). Though the parameters of this regression model are refit for each time point, correlation across time is preserved (perhaps excessively) in sampling by drawing (via Latin hypercube sampling) a single quantile of the variance characterized by the conditional standard deviation to use at all time points for a given time series sample. In sampling the t distribution, the number of degrees of freedom is taken as the number of GCMs providing ocean dynamic sea-level projections for a particular grid cell in the scenario used for calibration.

In some ways, the approach is similar to that of a linear-regression-based scaling of ocean dynamic sea level on global mean thermosteric sea-level rise, as in Palmer et al. (2020). The commonality is the assumption that the distribution of ocean dynamic sea level at a given point may be constrained by information about global mean thermosteric sea-level rise (“may” is an operative word here – it is also possible for the scaling factor or correlation coefficient to be zero).

One important difference is that this approach is recalibrated for each time step, whereas the Palmer et al. (2020) approach finds a single regression coefficient for a given GCM across time. A second is that the uncertainty not captured in the characterized correlation is sampled, whereas in Palmer et al. (2020), all variance is assumed to be captured by the spread of regression coefficients across GCMs. The approach used here is more focused on the distributional characteristics across GCMs, as opposed to representing each individual GCM by a regression coefficient. As a consequence of these differences, the Kopp et al. (2014) approach loses a degree of traceability to individual GCMs, instead being focused on preserving the distributional properties assessed based on the ensemble.

Note that where global mean thermosteric sea-level rise and ocean dynamic sea level are uncorrelated, this approach simply returns the multi-model mean and scaled standard deviation for the scenario.

Table A1 GMSL component projections for 2100 including AR6 projections.

Median (17th–83rd percentile) projections are shown relative to a 1995–2014 baseline. All are in meters. For certain modules, projections for Representative Concentration Pathways 2.6, 4.5, and 8.5 are shown in lieu of those for the SSP scenarios. AR6 results are taken from Fox-Kemper et al. (2021a) Tables 9.2, 9.3, 9.4, 9.8, and 9.9, as indicated by numbers in parentheses after the label.

Table A2 Total projections for 2100 compared to AR6.

Median (17th–83rd percentile) projections are shown in meters relative to a 1995–2014 baseline. AR6 values are taken from Fox-Kemper et al. (2021a) Table 9.9, except for low-confidence SSP1-2.6 values, which are taken from Garner et al. (2021). Table 9.9 results are based on workflows 1e and 2e (medium-confidence projections) as well as workflows 1e, 2e, 3e, and 4 (low-confidence projections).

Table A3 CMIP6 models used for calibrating the thermal expansion coefficients of Fox-Kemper et al. (2021a) (TE, left column) and for projecting ocean dynamic sea-level change and the IB effect (zos+psl, right column) in the tlm/sterodynamics module .

Figure A1 Comparison of mean rate estimates from the kopp14/verticallandmotion module, which uses the ICE5G ice history and VM2-90 viscosity profile (Peltier, 2004), with rate estimates derived using the same methodology but with a prior based on the ICE-6G-C ice history and VM5a viscosity profile (Stuhne and Peltier, 2015). (a) Differences in absolute rates (mm yr⁻¹) at tide gauges and grid cells. Pale areas have an absolute mean rate difference of <0.25 mm yr⁻¹. (b) For tide gauge locations, the rate estimate is derived using the alternative prior with the kopp14/verticallandmotion estimate. The green line is the 1:1 line. Uncertainties shown are ±1σ.