Uncertainty-enabled machine learning for emulation of regional sea-level change caused by the Antarctic Ice Sheet

Abstract

Projecting sea-level change in various climate-change scenarios typically involves running forward simulations of the Earth’s gravitational, rotational and deformational (GRD) response to ice mass change, which requires high computational cost and time. Here we build neural-network emulators of sea-level change at 27 coastal locations, due to the GRD effects associated with future Antarctic Ice Sheet mass change over the 21st century. The emulators are based on datasets produced using a numerical solver for the static sea-level equation and published ISMIP6-2100 ice-sheet model simulations referenced in the IPCC AR6 report. We show that the neural-network emulators have an accuracy that is competitive with baseline machine learning emulators. In order to quantify uncertainty, we derive well-calibrated prediction intervals for simulated sea-level change via a linear regression postprocessing technique that uses (nonlinear) machine learning model outputs, a technique that has previously been applied to numerical climate models. We also demonstrate substantial gains in computational efficiency: a feedforward neural-network emulator exhibits on the order of 100 times speedup in comparison to the numerical sea-level equation solver that is used for training.

Keywords Regional sea level change  $\cdot$ Emulator  $\cdot$ Machine learning

1 Introduction

Sea-level change is a major impact of climate change that threatens coastal communities and infrastructure (Fox-Kemper et al.,, 2021; Hauer et al.,, 2020). Global mean sea level (GMSL) rose by 0.20 m over the last century, and the recent Intergovernmental Panel on Climate Change report (IPCC AR6; Fox-Kemper et al., (2021)) indicates that it is "virtually certain" that GMSL will continue to rise through 2100, regardless of anthropogenic greenhouse gas emissions scenario, due to the likelihood of all major processes contributing to GMSL continuing to operate. Sea-level change varies spatially, and relative sea level (RSL, i.e., ocean depth) can differ by up to a factor of two relative to GMSL across regions of the globe (Fox-Kemper et al.,, 2021; Kopp et al.,, 2015). Variations in RSL are caused by steric effects (changes in ocean density and salinity) and dynamic sea level changes, as well as gravitational, rotational, and deformational (GRD) effects (Slangen et al.,, 2014; Kopp et al.,, 2015; Gregory et al.,, 2019; Hamlington et al.,, 2020; Fox-Kemper et al.,, 2021). Unlike the other two factors, GRD effects are induced by changes in surface mass loading such as ice sheets, glaciers, and terrestrial water storage (reservoirs, lakes, groundwater). That is, mass redistribution on the Earth’s surface perturbs the gravitational field and rotation potential of the Earth. In addition, the solid Earth deforms over time in response to the changing surface loads. These complex processes yield spatially and temporally varying sea-level change referred to as "sea-level fingerprints" unique to the geometry of ice mass change (Farrell and Clark,, 1976; Mitrovica et al.,, 2011). Future sea-level change and uncertainties in projections will be dominated by changes in the Antarctic Ice Sheet (AIS) (Kopp et al.,, 2017; Bakker et al.,, 2017; Edwards et al.,, 2021; Fox-Kemper et al.,, 2021), making the sea-level change due to GRD effects an important factor to accurately predict.

Projections of future AIS contribution to sea level are improving through coordinated community model intercomparison efforts (Nowicki et al.,, 2020; Seroussi et al.,, 2020; Edwards et al.,, 2021), as well as large ensembles of AIS simulations from single models that quantify uncertainties (Bulthuis et al.,, 2019; Berdahl et al.,, 2023; Coulon et al.,, 2024). These ice-sheet projections have been used in probabilistic frameworks for analyzing regional sea level (Jevrejeva et al.,, 2018; Kopp et al.,, 2023). However, these frameworks have so far assumed static sea-level fingerprints calculated by only accounting for the total ice thickness changes for each ice sheet and glacier region (Kopp et al.,, 2023). One could improve the projections by considering temporally evolving sea-level fingerprint, as recent studies that predict sea-level changes due to GRD effects associated with projected land-ice mass change for 2015-2100 CE have shown that the variable geometry of AIS mass changes has a significant impact on future sea-level change (Roffman et al.,, 2023; Cederberg et al.,, 2023). Because probabilistic sea-level change frameworks require evaluation of tens of thousands of samples (Wong et al.,, 2017; Kopp et al.,, 2023), they typically rely on fast-to-evaluate emulators or parameterizations of complex models and processes. Tools for rapidly calculating sea-level fingerprints due to GRD effects for different patterns of land-ice mass change would facilitate probabilistic sea-level projections.

Because of the complex relationships between mass loading changes and sea-level fingerprints (Mitrovica et al.,, 2011; Roffman et al.,, 2023; Cederberg et al.,, 2023) and costs of calculating time-evolving sea-level fingerprints using sea-level models, machine learning (ML) provides a promising avenue for creating computationally inexpensive emulators of the sea-level equation. Recent work has demonstrated the potential of a graph-based spherical convolutional neural network to emulate a sea-level model applied to ice-age glacial isostatic adjustment over the last deglaciation and achieved a 100–1000$\times$ speedup (Lin et al.,, 2023). Another study used an artificial neural network to develop an emulator for the effect of lateral variations in Earth rheological structure that are often left out of sea-level models (Love et al.,, 2023). ML has also been applied to the ice-sheet model projections that form the inputs to sea-level models. Van Katwyk et al., (2023) created an emulator for the contribution from the AIS to GMSL by training a neural network on an ensemble of ice-sheet model experiments. Rohmer et al., (2022) used a similar ensemble, in this case for the Greenland Ice Sheet, in conjunction with an ML-based attribution approach to quantify the influence of different ice-sheet modeling decisions to projections of future ice-sheet contribution to GMSL. ML methods have also been utilized to predict complex regional sea level signals from observed historical sea-level change datasets at month to year timescales (e.g. Tur et al.,, 2021; Ayinde et al.,, 2023; Nieves et al.,, 2021; Guillou and Chapalain,, 2021; Altunkaynak and Kartal,, 2021; Song et al.,, 2021) and applied to various individual components of the ice-sheet models (e.g. Brinkerhoff et al.,, 2021; Jouvet et al.,, 2022; Jouvet,, 2023; He et al.,, 2023; Riel et al.,, 2021; Riel and Minchew,, 2023; Rosier et al.,, 2023; Verjans and Robel,, 2024). However, the use of ML to accelerate sea-level fingerprint calculations for projections of future regional sea level has yet to be attempted.

Here, we develop a set of neural-network ML emulators (sometimes referred to as surrogate models) of the sea-level equation for predicting regional sea-level changes caused by changes in the mass of AIS and compare their performances to baseline Gaussian process and random forest based emulators. The emulators are trained on a set of sea-level model simulations forced with a recent ensemble of ice-sheet model projections of AIS through 2100 (ISMIP6, Serrousi et al., 2020). Rather than predicting the entire global field of sea-level change, we train the emulators to only consider the locations of 27 major coastal cities. In this paper, we begin by comparing four different ML methods for building emulators and evaluate their strengths and weaknesses for this application. We then develop an approach for quantifying uncertainty in the ML emulators using a simple linear regression technique to estimate prediction intervals of each emulator for each city, a method that has been used for calibrating numerical climate models (Glahn and Lowry,, 1972; Glahn et al.,, 2009; Lovegrove and Siegert,, 2023). Finally, we provide an illustration of applying the emulators to generate probabilistic projections of regional sea-level change caused by AIS mass changes by generating a large ensemble of artificial ice-sheet projections.

2 Data description

Training the ML emulators required datasets of ice-sheet mass change projections and the corresponding regional sea-level fingerprints, which were divided into training, validation and test groups. The following sections describe the generation of these datasets in detail.

2.1 Ice-sheet mass change dataset

The input dataset for our workflow was an ensemble of future AIS mass change calculated by ice-sheet models. For this, we adopted the results of the AIS experiments from the Ice Sheet Model Intercomparison Project for Coupled Model Intercomparison Project Phase 6 (ISMIP6) (Seroussi et al.,, 2020). The ISMIP6-Antarctica protocol specified standard experiment configurations for ice-sheet models to generate projections of AIS from 2015 through 2100 (Nowicki et al.,, 2020). The ice-sheet model simulations were forced by atmosphere and ocean output from six different Coupled Model Intercomparison Project Phase 5 (CMIP5) climate models for the high greenhouse gas emissions RCP8.5 scenario and two climate models for the low emissions RCP2.6 scenario. In addition to testing different climate forcings, different experiments tested different representations of ice-shelf basal melt and the potential for removal of ice-shelf regions through hydrofracture. In total, there were 21 experiments defined; here we select eight experiments for which standardized protocols were prescribed. The experiments were executed by 13 modeling groups around the world, and the results were analyzed on a common reference grid. Each modeling group also performed a control run through 2100 using steady constant climate forcing to quantify model drift. The ensemble represents a wide variety of ice-sheet model complexity, grid resolution, and model characteristics. Notably, none of the ice-sheet models in the ISMIP6 ensemble employed any of GRD effects.

To prepare these ice-sheet projection data as input to a sea-level model, we performed a number of processing steps: masking out floating ice and removing model drift from projections on the 8 km-resolution standard polar stereographic grid, and interpolating the data onto a global Gaussian grid of spherical Earth. Since changes in floating ice do not contribute to sea level (except by a negligible amount due to density difference between ice and ocean water), ice thickness input to the sea-level model should only comprise grounded ice. We calculated changes in the ISMIP6 spatial field “land_ice_thickness” from its initial value. We restricted the calculation to the region of the domain where ice is grounded, which we calculate based on the flotation criterion, assuming an ice density of 910 kg m^-3 and an ocean density of 1028 kg m^-3. Second, because there is a large drift in the control run of many simulations (Aschwanden et al.,, 2021), we removed the thickness change in the control run for each model for each projection, as was done in the ISMIP6 analysis (Seroussi et al.,, 2020). This step is particularly important for modeling regional sea-level changes, as some portion of ice thickness change from raw projections will be related to the bias inherent to each ice-sheet model (i.e., model drift), which we do not want to feed into our sea-level fingerprint calculation. While ISMIP6 output is posted annually, we subsampled it to 5 year intervals because annual changes are relatively small, leading to 17 time levels of grounded ice thickness change. Finally, we remapped the ice thickness change fields from the 8 km polar stereographic projection to the $512\times 1024$ latitude/longitude Gaussian grid used by the sea-level model using a conservative remapping method to ensure total AIS mass change is unaltered by the remapping.

A total of 194 AIS projections were generated for ISMIP6 across the different experiments and models used, and here we restrict ourselves to 71 projections derived from the eight experiments with standardized protocol (experiments 5-10 and 11-12), eliminating a small number of projections due to processing errors related to issues with data format conventions. Accounting for the 17 time levels retained for each projection, we had a total of 1207 fields of grounded ice thickness change in our input dataset used for training, validation, and testing of ML emulators.

2.2 Sea-level dataset: Calculation of sea-level fingerprints

We generated sea-level fingerprints associated with different patterns of AIS mass loss using the sea-level model described by Han et al., (2022), an updated version of sea-level models described in Kendall et al., (2005) and Gomez et al., (2010). The model computes barystatic sea-level changes due to GRD effects associated with AIS mass changes by solving the static sea-level equation (Farrell and Clark,, 1976; Mitrovica and Milne,, 2003) using the pseudo-spectral algorithm described in detail in Kendall et al., (2005), with truncation up to spherical harmonic degree and order 512. The sea-level framework accounts for shoreline migration, including inundation of presently glaciated marine basins as the extent of the grounded ice sheet retreats, as well as perturbations in the Earth’s rotation axis (Mitrovica and Milne,, 2003; Mitrovica et al.,, 2005).

In addition to the projected patterns of AIS ice mass loss, the simulation requires a model for the rheological structure of the solid Earth. We adopt a radially symmetric Maxwell viscoelastic model for the Earth with elastic and density properties derived from the seismic PREM model (Dziewonski and Anderson,, 1981). The depth varying viscosity of the Earth is given by the viscosity profile VM7, which is based on recent inferences from wide range of global geophysical observables (Roy and Peltier,, 2017, 2018). We note that these values represent a global average; while appropriate for East Antarctica and our 27 far field coastal locations, this Earth structure does not accurately capture the Earth structure in West Antarctica where upper mantle viscosity is thought to be 2-3 orders of magnitude lower than the global average (Barletta et al.,, 2018; Nield et al.,, 2014). We also note that while it is well-established that the Earth structure in Antarctica shows strong lateral heterogeneity (Lloyd et al.,, 2023) likely having an effect on projections of local sea-level change, the sensitivity of far-field sea level to the 1D versus 3D Earth structure only becomes large enough to consider in a long-term (multi-centennial or longer) scale projection (Pan et al.,, 2021; Powell et al.,, 2022), which is beyond the timescale we consider in this work. While we include viscoelastic deformation of the solid Earth in the model, with the chosen viscosity profile over centennial timescales, the viscous component is significantly smaller than the elastic component, yielding results similar to assuming an elastic Earth.

Since the sea-level model runs on a global grid but the original ISMIP6-AIS datasets were provided on a local grid, we adopt topography from the ETOPO2 dataset (NOAA National Centers for Environmental Information,, 2022) for the region outside of the ice-sheet model domain. With the pre-processed ice-sheet thickness change time-series and the Earth structure described above, we perform a suite of sea-level simulations with 5-year time intervals and generate sea-level change output datasets for emulator training, validation, and testing. Examples of the maps of sea-level change associated with AIS thickness changes are illustrated in Figure 1.

We initially attempted to emulate the entire $512\times 1,024=524,288$ dimension output space of the sea-level model, with the best-performing model using an autoencoder (Kramer,, 1991) for dimension reduction on inputs and outputs and a feed-forward neural network linking between reduced-dimension spaces. Our initial results varied considerably in quality – while some of the entire fields in the test set could be emulated reliably, about 10% of the global test fields were poorly emulated (normalized root mean squared error exceeding 0.11). Moreover, the latitude-longitude coordinates of the sea-level model spatial output added complexity in accurately representing a spherical surface in two dimensions (Lin et al.,, 2023). To create a more tractable but useful demonstration case, we chose to only evaluate RSL change at a subset of coastal locations. We chose the 27 long-term tide gauge sites identified by Meyssignac et al., (2017) (Fig. 1c), which come from the Permanent Service for Mean Sea Level (Holgate et al.,, 2013).

Figure 1: Illustration of RSL fingerprints for given AIS thickness changes for two samples. a) Modeled AIS thickness change between 2015 and 2100 for sample from training dataset with third largest total AIS mass change. b) Calculated RSL change fingerprint for AIS mass loss in a. c) Modeled AIS thickness change between 2015 and 2100 for sample from training dataset with third smallest total AIS mass change. d) Calculated RSL change fingerprint for AIS mass loss in c. Note different colorbar ranges for each plot and that values exist beyond the ranges shown. The 27 long-term tide gauge sites selected by Meyssignac et al., (2017) that we used as emulation targets are shown with black dots in b and d.

2.3 Partitioning of training, validation, and test datasets

The entire dataset of ISMIP6 AIS mass change fields, along with their corresponding numerically solved sea-level change outputs, was split into training, validation, and test datasets, as is commonly done in ML (Hastie et al., (2009)), with 80% used for training, 10% used for validation, and 10% used for a hold-out test set. As mentioned above, this constitutes 1,207 pairs of AIS mass change and global relative sea-level change fields across the different ice-sheet models, experiments, and time levels.

In order to ensure representation of each year across the three data splits, the input-output pairs from the 71 models per each year were randomly selected to be in the training, validation, and test sets. That is, for a given year, 57 were randomly selected to be in the training set, 7 were randomly selected to be in the validation set, and 7 were randomly selected to be in the test set. The final split of data consisted of 969 in training, 119 in validation, and 119 in testing, where the inputs of each split were standardized to have 0 mean and unit variance. Before applying a variety of ML models (described in Section 3) to this dataset, we preprocessed the sea-level dataset by removing those elements that were zeros throughout the entire dataset (i.e., cells without grounded ice at any time, almost two-thirds of the grid cells, Fig. 1a,b). After preprocessing, the length of the flattened input vector was 205,928, and the output vector was the 27 long-term tide gauge locations.

3 ML models

Our overarching objective is to construct a machine learning (ML) model that takes as input a vector of ice thickness changes across Antarctica derived from ISMIP-6 simulations and outputs a prediction of the simulation of resultant RSL change at 27 long-term tide gauge locations (Fig. 1c,d). That is, in our RSL change emulator, the input is the flattened vector of ice-thickness changes, and the output is a prediction of the resultant RSL change for the 27 long-term tide gauge locations. Ideally, the ML model should accurately predict simulated RSL change values, be computationally more affordable than running the full numerical simulation, and also enable statistically valid uncertainty quantification of predictions. In what follows, we describe a variety of ML implementations that are constructed and compared: a feedforward neural network (NN), a conditional variational autoencoder (CVAE), a Gaussian process (GP) acting on a reduced-dimensional space, and a random forest (RF) acting on a reduced-dimensional space. Our main finding is that the NN, CVAE, and GP models yield similar levels of accuracy for this application, but the NN and CVAE yield the best uncertainty quantification in terms of prediction intervals using a simple linear regression postprocessing technique; moreover, prediction intervals are of a useful length. In what follows, we describe these ML approaches and our implementations in more detail.

3.1 Feedforward neural network (NN)

A feedforward neural network (NN) is a fundamental component of deep neural networks. The primary objective of an NN is to learn (or approximate) a mapping $\bm{y}=\bm{f}(\bm{x};\bm{\theta})$ from input $\bm{x}$ to output $\bm{y}$, accomplished through the composition of functions governed by parameters $\bm{\theta}$ that consist of weights and bias terms. The approximation task is equivalent to determining the best estimate $\hat{\bm{\theta}}$ that minimizes a loss function (e.g., mean-squared error) using a training set, which is typically updated iteratively with backpropagation (Rumelhart et al.,, 1986).

It is common to use well-known optimization routines such as stochastic gradient descent (SGD) (Bottou,, 2010) and ADAM (Kingma and Ba,, 2017) for fitting neural networks. A feedforward neural network (NN) consists of an input layer, $L$ hidden layers, and an output layer, where $L$ denotes the number of hidden layers. $\ell^{th}$ hidden layer is equipped with $d_{\ell}$ neurons, which serve as the dimensions of hidden states.

The layers’ units (neurons) are interconnected through a linear transformation and an activation function, denoted as $g(\bm{x}^{\top}\bm{w}+\bm{b}$), where $\bm{w}$, $\bm{b}$, and $g(\cdot)$ correspond to the weight, bias, and an activation function, respectively. A sequence of these operations across each layer can yield a high-quality approximation to complex functions. The potential for such approximation is enhanced by the activation function, such as the rectified linear unit (RELU), which introduces non-linearity. In addition to its rich approximating capabilities, user-friendly software, such as PyTorch by Paszke et al., (2019), has contributed to the popularity of NNs. For a more detailed introduction to NNs, refer to Goodfellow et al., (2016).

In our implementation, we utilize a neural network architecture comprising $L=4$ hidden layers. We use the RELU activation function between layers, while a linear activation function is applied between the last hidden layer and the output layer. The number of neurons within each layer is $(d_{1},d_{2},d_{3},d_{4})=(3200,1600,800,400)$. Using root mean squared error (RMSE) as the loss function, we train the model with the SGD optimizer with a learning rate of 0.001. Throughout the training process, we monitor the convergence of the loss on the validation set. To prevent overfitting, we also introduce dropout (Srivastava et al.,, 2014) with a probability 0.3 between layers. We implement the NN using Pytorch version 2.0.1 with Python version 3.9.12.

3.2 Conditional variational autoencoder (CVAE)

The CVAE is a widely-used modern neural network architecture that couples an encoder to a decoder, with the intermediate layer between them representing a latent, lower-dimensional space. The CVAE is appealing and can be effective in constructing generative models (e.g., a molecular generative model (Lim et al.,, 2018) and zero-shot learning generative model (Mishra et al.,, 2018)). We consider a CVAE in addition to the feed-forward NN in order to attempt a generative approach for the sea-level change emulator, since by sampling from the appropriate latent space, one can generate samples of predictions for a given input as opposed to a single prediction; indeed this approach serves as a precursor for generative artificial intelligence (AI) approaches such as diffusion models (Ho et al.,, 2020).

The origin of the CVAE can be traced back to the autoencoder (AE) (Kramer,, 1991) and the variational autoencoder (VAE) (Kingma and Welling,, 2013). In general, an AE is trained to reproduce the original input, thereby learning a reduced-dimension latent representation that captures much of the information of the higher dimensional input data. Although the architecture of the VAE is similar to the AE, the VAE assumes that the reduced-dimensional latent space follows a certain distribution, often a multivariate normal with 0 mean and identity covariance matrix. The key advantage of the VAE over the AE is that it can serve as a generative model by drawing a sample from the distribution of the latent, intermediate layer.

While both the CVAE and the VAE are similar, the CVAE possesses an additional feature: the ability to instruct the model to predict the output based on conditional information. In essence, a prediction for a new input is determined by concatenating a random draw from the latent space with conditional information, which in our setting is the vector of ice-thickness changes. This concatenated vector, which combines conditional information with latent information, is passed through the decoder to make a prediction. For a visual representation, see Figure 2; more details are described in the tutorial of Doersch, (2016).

Figure 2: CVAE example with an intermediate layer that follows a multivariate normal distribution. Note that conditional information is provided in the decoder.

Our implementation of the CVAE for the sea-level change emulator takes the entire vector of ice-thickness changes as conditional information, concatenates it with a latent draw, and passes the concatenated vector through the decoder to make a prediction. Specifically, the encoder comprises $L=4$ hidden layers with neurons ($d_{1}$, $d_{2}$, $d_{3}$, $d_{4}$) = (3200, 1600, 800, 400), and it outputs 200 pairs of parameters ($\mu_{i}$, $log\sigma^{2}_{i}$), where $i=1,...,200$. These parameters correspond to the parameters of the posterior intermediate layer distribution $P(\bm{Z}|\bm{X})$, which is assumed to follow a multivariate normal distribution with a mean vector $\bm{\mu}=(\mu_{1},...,\mu_{200})$ and a diagonal covariance matrix $diag(\sigma^{2}_{1},...,\sigma^{2}_{200})$. Then, the decoder takes as input the concatenated vector of the ice-sheet thickness vector (i.e., the conditional information) and the transformed samples from the prior distribution $P(\bm{Z})$ for the intermediate layer. Assuming $P(\bm{Z})$ follows a multivariate normal distribution with 0 mean and identity covariance matrix, the transformation is achieved through re-parametrization $Z_{i}|\bm{X}=\mu_{i}+\sigma_{i}Z_{i}$. Similar to the encoder, the decoder also consists of $L=4$ hidden layers with the same number of neurons as the encoder. It outputs the predicted sea-level change at 27 long-term tide gauge locations. Because the CVAE is predicting sea-level change values as opposed to reconstructing its own input, this is a modification of the usual AE approach.

The loss function used to train is a sum of two terms: the first is a KL-divergence between the variational approximation of posterior latent distribution $P(\bm{Z}|\bm{X})$ and prior distribution $P(Z)$, which can be thought of as a regularization penalty. The second is a mean squared error term for predictions from the decoder. In our implementation, we draw 500 samples from the latent distribution $P(\bm{Z})$ during the prediction generation process, and the final predicted value is obtained by calculating the mean of these 500 predictions. The CVAE is implemented using PyTorch version 2.0.1 with Python version 3.9.17.

3.3 Baseline comparisons using Gaussian process and random forest

We use baseline comparison emulators to assess the performance of the NN and CVAE emulators. These baseline comparisons are a stationary Gaussian process with a squared-exponential covariance function (GP) and a random forest (RF). Gaussian processes have been a common choice to build emulators for physical applications (Higdon et al.,, 2008; Gu et al.,, 2018) and hence serve as a natural comparison to NN and CVAE emulators. They are primarily used because of closed-form multivariate normal distributions conditioning on training data; see Gramacy, (2020) for exact formulations.

For the baseline GP emulator, we use a GP that maps from a reduced dimensional input space to the sea-level change at each of the 27 long-term tide gauge locations. The dimension-reduction technique we adopt in this work is discussed subsequently. The GP uses a squared-exponential kernel (Rasmussen,, 2003; Gramacy,, 2020) whose parameters, including a nugget term, are fit with a statistical technique known as profile likelihood maximization, the mathematical steps of which are described in Chapter 5 of Gramacy, (2020). Our implementation of profile likelihood was coded in the R programming language. While there exist some recent variants of GPs such as deep Gaussian processes (e.g., Sauer et al., (2023)) that are capable of handling nonstationarity, we use a stationary GP with a squared-exponential covariance function only to serve as a baseline comparison.

RFs have been widely adopted as off-the-shelf ML models, making an RF model an ideal additional baseline comparison (e.g., see Hooten et al., (2011) and Gopalan and Wikle, (2022) for emulators of physical processes based on RFs). A detailed introduction to RF models is provided in Hastie et al., (2009) and James et al., (2013). For baseline comparisons, we use the “RandomForestRegressor" function in a package “scikit-learn" (Pedregosa et al.,, 2011) to train an RF in Python 3.9.17. We use 100 decision trees (default) and mean squared error (MSE) as the loss function for this baseline RF. As is the case for the GP baseline, an RF is trained for each of the 27 long-term tide gauge locations using the reduced dimensional input.

While the NN and CVAE models take as input an entire vector of ice thickness changes across Antarctica and perform dimension reduction within their respective architectures, the GP and RF models use a reduced-dimension space as input. Dimension reduction was performed by projecting input vectors onto a reduced-rank principal-component basis; 250 principal components were used, which represent almost 98.7% of the total variation of the training set inputs. This technique is standard for emulating high-dimensional computer simulator output (e.g., Higdon et al., (2008), Hooten et al., (2011), Gopalan and Wikle, (2022)), but in this work we use dimension reduction on the space of inputs for the GP and RF emulators, which is likely required as the dimension of inputs grows substantially – the input vector length is on the order of 100,000 whereas there are on the order of 1000 data vectors. For instance, computing the distance matrix needed for a GP model is an $O(n^{2}m)$ operation when there are $n$ length $m$ vectors, which for $m>>n$ can be more substantial than the $O(n^{3})$ inverse and determinant calculations needed for GPs.

4 Uncertainty quantification

We construct prediction intervals as a means of quantifying uncertainty in emulated sea-level change values at the 27 long-term tide gauge locations. For the NN, CVAE, and RF emulators, we use a simple linear regression postprocessing technique, discussed in what follows, in order to generate prediction intervals. For GPs, the distribution of predictions given training points is multivariate normal with exact mean and covariance (Rasmussen, (2003); Gramacy, (2020)). Prediction intervals can thus be derived for the GP emulator by using the appropriate quantiles of the resultant Gaussian distribution. Before explaining the linear regression technique, we note that we found the empirical prediction intervals solely constructed with Monte Carlo dropout for both the NN and CVAE were far too narrow; a similar phenomenon was exhibited with generated samples from the CVAE.

As reviewed in Lovegrove and Siegert, (2023), linear regression methods have been used to postprocess outputs from numerical climate models by using the numerical model output as an input variable (i.e., covariate) for either a simple or multiple linear regression, the latter case when there are multiple covariates. In the weather forecasting literature, this has been referred to as model output statistics (Glahn and Lowry,, 1972; Glahn et al.,, 2009). Lovegrove and Siegert, (2023) show with a global temperature forecasting example that fitting simple linear regressions yields prediction intervals for temperature that have close to nominal coverage probabilities; i.e., a (1-$\alpha$)% prediction interval as derived from the simple linear regression covers the true value with probability near (1-$\alpha$)%. In addition, such an approach has been used for calibrating probabilities generated from ML models for categories via logistic regression, referred to as Platt scaling (Platt et al.,, 1999).

Our approach follows this idea: for a given location, we regress the simulated sea-level change values onto the ML predictions for all of the simulation data in the validation set using a simple linear regression. This is also similar to the split-conformal method of Lei et al., (2018) because we use a dataset split from the original training data to fit the linear regression models, as to avoid overfitting by training both the ML model and the simple linear model on the same data. However, we are able to obtain conditional coverage by fitting different linear regressions for each location. Moreover, additional covariates could be included into the linear regression to improve conditional calibration. Additionally, a per-location bias is accounted for by the intercept term of the simple linear regression models.

We use the standard prediction interval from the fitted simple linear regression model:

$$\widetilde{y}_{ij}^{test} \pm t_{1-\alpha/2,n-2}\times$$

$$\sqrt{S_{i}^{2,val}\times\bigg{(}1+\frac{1}{n}+\frac{(\hat{y}_{ij}^{test}-\bar{\hat{y}}_{i}^{val})^{2}}{\sum_{j}(\hat{y}_{ij}^{val}-\bar{\hat{y}}_{i}^{val})^{2}}\bigg{)}},$$ (1)

where $\hat{y}_{ij}^{test}$ denotes the prediction from the ML model for the $i^{th}$ location and $j^{th}$ test sample, $i=1,..,27$, $j=1,...,119$. Additionally, $\bar{\hat{y}}_{i}^{val}$ and $\widetilde{y}_{ij}^{test}$ indicate the location-wise mean of the prediction from the ML model $\hat{y}_{ij}^{val}$ and the calibrated point estimate derived through a simple linear regression:

$$\widetilde{y}_{ij}^{test}=\hat{\beta}_{i,0}^{val}+\hat{\beta}_{i,1}^{val}\times\hat{y}_{ij}^{test},$$ (2)

where $\hat{\beta}_{i,0}^{val}$ and $\hat{\beta}_{i,1}^{val}$ are derived by fitting a simple linear regression on the validation set. Also, a variance estimate for the $i^{th}$ location, $S^{2,val}_{i}$, is obtained by

$$S^{2,val}_{i}=\frac{1}{n-2}\sum_{j}(y_{ij}^{val}-\widetilde{y}_{ij}^{val})^{2},$$ (3)

where $y_{ij}^{val}$ represents the (simulated) sea-level change for the $i^{th}$ location and $j^{th}$ validation sample, while $\widetilde{y}_{ij}^{val}$ denotes the calibrated point estimate obtained through linear regression. Lastly, $t_{1-\alpha/2,n-2}$ denotes the 1-$\alpha$ quantile of a $t$ distribution with $n-2$ degree of freedom, where $n=119$ represents the sample size (i.e., the validation set size). We emphasize that we utilize the validation set not only to monitor potential overfitting but also to derive estimates for a simple linear regression. In other words, $\hat{y}_{ij}^{val}$, $\bar{\hat{y}}_{i}^{val}$, $S_{i}^{2,val}$, $\hat{\beta}_{i,0}^{val}$, and $\hat{\beta}_{i,1}^{val}$ are estimated based on the validation set, and these estimates are then employed to calibrate the ML predictions.

So long as a linear model is an appropriate fit between observations (i.e., the sea-level change simulations) and ML predictions, prediction intervals should cover with nominal probabilities. In Section 5 that follows, we demonstrate empirical coverage rates on a hold-out test set of simulations. While we do not explore alternate residual distributions or spatial correlation in the residuals, we believe these items would be worth exploring in future works.

5 Emulator results

Figure 3: R-squared comparison between calibrated ML model predictions and actual simulated sea-level change values. R-squared values indicate comparable accuracy for the NN, CVAE, and GP models, though the RF model performs worst.

Figure 4: Comparison of empirical prediction interval coverage rates, suggesting that the CVAE and NN have coverage closest to the nominal rate. Note that since there are on the order of 100 examples in the test hold-out set, it is difficult to distinguish between coverage rates of 1%.

Our main results consist of a comparison of accuracy, empirical prediction interval coverage rate, and prediction interval lengths computed on the test hold-out set. We calculate Pearson’s correlation coefficient of determination, $R^{2}$, between calibrated ML model predictions and the actual simulated sea-level change values as a measure of accuracy. The boxplots of Figure 3, 4, and 5 depict the range of values for each of the models across the 27 distinct global locations.

From Figure 3, it is evident that the CVAE, GP, and NN models exhibit similar prediction performance, while the RF model shows inferior performance. The average $R^{2}$ values for each model are 0.972 (NN), 0.969 (CVAE), 0.971 (GP), and 0.844 (RF). This suggests that the CVAE and NN can serve as alternatives to GP-based emulators, which have been widely used.

Figure 5: Average length comparison of prediction intervals in meters – the Gaussian process intervals are shortest but are overconfident as based on coverage rate. The RF tends to be longest.

Figure 4 and 5 show the range of (1-$\alpha$)% empirical coverage rates and the length of prediction intervals for the significance level $\alpha$=0.01, 0.05, 0.1, and 0.15. From Figure 4, one can see that the empirical coverage rates for the GP model underestimate the nominal rate more than other ML methods after calibration. It is noteworthy that the RF still exhibits some deviation in empirical coverage from the nominal rates, even after calibration. This is due to the less accurate prediction performance of the RF model. While the aim of calibration is to enhance the accuracy of prediction intervals and point estimates, its effectiveness can be compromised if the underlying model lacks prediction performance.

Figure 5 indicates that the lengths of the GP prediction intervals are generally the shortest but compensate with undercoverage; i.e., a smaller empirical coverage probability than is nominally stated for the four levels shown. On the other hand, due to the overall inaccuracy, the RF tends to have the longest intervals as calculated with simple linear regression, emphasizing the importance of the underlying model.

The prediction intervals from the NN and CVAE are overall quite similar in terms of length, coverage, and accuracy. It appears that the CVAE is essentially behaving almost identically to the NN; we believe this is because the conditional ice-thickness input is providing much of the information in comparison to the latent variable for the prediction, and the decoder architecture of the CVAE is the same as the NN. Overall, the NN and CVAE models seem to perform best taking into account all three characteristics being compared.

As an illustration of the simple linear regression approach for determining prediction intervals, Figure 6 shows a scatterplot depicting simulated sea-level change on the $y$-axis and NN predictions on the $x$-axis for the location Balboa. It is noticeable that the $y=x$ line (solid red) closely aligns with the simple regression line (solid blue), indicating a high degree of accuracy in the predictions ($R^{2}\approx 0.974$). Furthermore, one can see that a simple linear regression technique can yield well-calibrated 95% prediction intervals (solid green).

Figure 6: A scatter plot featuring predictions on the X-axis and simulated sea-level change on the Y-axis. Solid red line is the $y=x$ line while solid blue is the simple regression line. Solid green lines show prediction intervals for significance level $\alpha=0.05$.

Figure 7: Runtime Comparison. For comparison the full numerical model requires about 60 seconds per solution.

In addition, we considered a variation where the inputs are mean ice-thickness values by distinct geographic sectors in Antarctica where the sectors are illustrated in Figure 9 of Appendix A; using sector-based inputs has a precedent in Van Katwyk et al., (2023). However, for this particular emulator, we found that the sectors performed markedly worse for all models using the sector mean values as input. These comparisons can be seen in the Supplementary Materials.

We additionally assess the runtime for generating predictions with each of the ML models by performing a repeated single prediction task 1,000 times. Note we consider that a single prediction is complete when the model provides the predicted values for all 27 locations. For the CVAE model, the runtime is measured based on a single sample drawn from the assumed latent distribution. The runtime for the GP model is evaluated after its kernel parameters are estimated. Figure 7 indicates that all models can generate predictions within one second. The NN, CVAE, and RF models exhibit comparable and more efficient performance than the GP model. Nonetheless, it is noteworthy that when 500 samples from the latent distribution are utilized for the CVAE, as per our analysis, the runtime extends to approximately 125 seconds. We also observe that the additional calibration step does not contribute to a significant increase in runtime ($\approx$ 0.0006 seconds). Since a single run of the sea-level solver can take about 60 seconds, the demonstrated speed up is greater than two orders of magnitude, particularly for the NN emulator. Training each emulator takes on the order of a few hours on top of the overheard of  21 hours (17 minutes per a single run, times 71 runs) for running the sea-level fingerprint calculations in serial.

6 Illustration of emulator use for a Monte Carlo ensemble

We demonstrate a potential use of the sea-level change emulator for generating distributions of sea-level change outputs for locations across the globe, in order to show how such an emulator could be used to produce RSL distributions from future Antarctic ice-sheet simulations. The most straightforward Monte Carlo approach to do so is to sample many parameters from their respective distributions, run an ice-sheet simulator under the sampled parameters, and subsequently run the resultant ice-sheet states through the sea-level change emulator in order to generate a distribution of sea-level change values. For our purposes, running an ice-sheet simulator iteratively is prohibitively expensive and so we instead construct a Monte Carlo ensemble of ice-sheet states using the ISMIP6 ensemble of Antarctica simulations.

Our goal in the demonstration is to use the neural-network emulator to produce a distribution of sea-level change values for the 27 long-term tide gauge sites. For this illustration, we focus on the year 2100 and high-emissions scenarios, and cull 48 instances from the ISMIP-6 training data matching these criteria. We performed a principal-components analysis on these 48 instances, retaining the top 38 which explained 99% of the variation in the year 2100, high-emissions simulations. We then bootstrapped the principal-component weights 100 times in order to randomly generate a Monte Carlo ensemble of Antarctic ice-sheet states representative of the year 2100 and high emissions. Alternative dimension reduction techniques could be explored, though a principal-components approach was used as an illustrative method to generate a Monte Carlo ensemble. These 100 randomly generated ice-sheet states were then fed into the neural-network emulator. We illustrate emulated densities of sea-level change using the Monte Carlo simulations in Figure 8. We emphasize that these distributions are not calibrated RSL projections and are meant as a demonstration of the application of the NN emulator.

The three sites chosen were selected as they represent among the widest variation in RSL of the cities sampled (Figure 1). The emulated RSL derived from the synthetic ice-sheet states span a similar range of sea-level change as the GMSL estimates from the ISMIP6 ensemble from which they were derived; Seroussi et al., (2020) report mean sea-level contribution of about 0.025 m with a standard deviation of 0.049 m for the high emissions scenarios. Seroussi et al., (2020) also demonstrate a longer tail of positive GMSL change in the ensemble, which our distributions also show (Figure 8). We consider this consistent given we did not use all ISMIP6 samples, our synthetic ice-sheet states are not guaranteed to have the same statistical properties as the ISMIP6 ensemble, we are considering RSL rather than GMSL, and these results come from the NN emulator.

Figure 8: Densities of emulated RSL change for Dunedin, Montevideo, and Midway during year 2100 with high-emissions scenario. These densities are an example of the application of the NN emulator and do not represent calibrated projections of RSL.

While there are clear differences in the densities of RSL projections for the three locations, they share many similarities, despite being from different parts of the Earth. The mode of each density close to zero is a result of our synthetic ice-sheet states generating many samples with small net mass loss from AIS. However, the tails of the distributions exhibit significant differences. Positive RSL values have lower magnitude for the two Southern Hemisphere locations (Dunedin and Montevideo), consistent with the fingerprint of decreased RSL rise closer to AIS (Figure. 1). This has the effect of compressing the range of the density functions. The normalized RSL changes at southern hemisphere locations is more sensitive to the spatial pattern of AIS mass change (Cederberg et al.,, 2023; Roffman et al.,, 2023), but due to the smaller magnitude of RSL change there, also have smaller standard deviations (Cederberg et al.,, 2023). While more detailed evaluation of these synthetic distributions is not appropriate, this example shows the potential utility of the NN emulator for producing RSL estimates rapidly from very large numbers of ice-sheet simulations.

7 Discussion

7.1 Comparison to related ML sea-level applications

The most similar ML application to sea-level change is the “Graph Neural Network Based EmulatOR for Glacial Isostatic Adjustment” by Lin et al., (2023). This work used a similar sea-level model (equivalently called a glacial isostatic adjustment model) to generate a training dataset of RSL change from a range of ice-sheet mass change trajectories, which was used to construct a CNN. A key difference is that their training dataset used paleo ice-sheet changes over tens of thousands of years, targeting a different application. They also utilized a graph-based spherical CNN, which allows convolution operations on a spherical manifold appropriate for representing the entire Earth. That is a promising technique for whole-Earth RSL emulation, and their emulator presumably could be trained on future AIS projections (and other sea-level contributors). For quantifying uncertainty, Lin et al., (2023) utilized an ensemble learning approach, training 30 versions of their emulator. However, their 99% prediction intervals cover the true values reportedly 94.5% of the time, suggesting that actual coverage rates are less than nominal coverage rates. Additionally, such an ensemble approach has the increased computational burden of having to refit the neural network multiple times.

While this study focuses on 27 specific locations, it is of interest to build a RSL emulator that produces spatial output for the entire globe, where the dimension of the output $\bm{y}$ significantly increases thereby introducing more complexity. Our preliminary work on this global spatial emulation problem has yielded varying quality of emulated outputs. Better performance could be achieved by exploring modern NN architectures, including convolutional NNs such as U-Net (Ronneberger et al.,, 2015) and the transformer (Vaswani et al.,, 2017), which could prove beneficial by effectively capturing spatial dependencies but come at the potential cost of generating a prohibitively large training data set of simulations. An additional challenge of emulating the entire globe is ensuring spherical consistency of the emulated RSL field. While the native output of the sea-level model is a structured latitude-longitude grid, treating that as two-dimensional space would ignore the fact that the longitude points are periodic and lines of latitude converge at the poles. Using a spherical representation of the output space (Lin et al.,, 2023) or emulating the spherical harmonic coefficients directly are avenues for addressing this issue.

Beyond the small dimension of our output space, there are ways in which the training data could be improved. The sample of the ISMIP6 ensemble we used here was relatively small and exhibited a fairly limited range of potential AIS states. Recent single-model ice-sheet model ensembles have produced hundreds (e.g. DeConto et al.,, 2021; Berdahl et al.,, 2023; Coulon et al.,, 2024; Aschwanden and Brinkerhoff,, 2022) or thousands (e.g. Bulthuis et al.,, 2019) of model realizations, and by exploring parameter uncertainty produce a large spread of potential ice-sheet states. Utilizing more expansive training data would improve both the accuracy of the emulator as well as its applicability, but with the trade-off of increased cost of generating the RSL training data with the sea-level model. A more selective approach to choosing training data could also make the process more efficient — because many of the ensemble members had mass change close to zero, the emulators may have been focusing on behavior of little interest.

Other areas for improvement include extending the RSL emulation to additional source of sea-level change (Greenland Ice Sheet, mountain glacier and icecaps, terrestrial water storage) and considering more complex Earth structure. As mentioned above, lateral variations in mantle viscosity and lithospheric thickness (Barletta et al.,, 2018; Nield et al.,, 2014) have significant impacts on RSL derived from AIS (Pan et al.,, 2021; Powell et al.,, 2022). Models accounting for 3-D viscoelastic Earth structure are significantly more computationally expensive to run (e.g., factor of $\mathcal{O}(10^{4})$, Love et al.,, 2023), which would greatly increase the value of developing fast emulators. Love et al., (2023) demonstrated potential for ML to capture the effects of lateral variations in Earth structure at a fraction of the computational cost.

7.2 Recommendations

Our results, and the others aforementioned, have shown that NNs are viable options for sea-level change emulators. Software such as PyTorch (Paszke et al.,, 2019), amongst others, has enabled the straightforward use of a variety of neural network architectures and different types of loss functions, creating many possibilities for the scientist.

Despite their substantial predictive capabilities, designing the architectures of both the NN and CVAE, including parameters like the number of layers ($L$), the quantity of units ($d_{\ell}$), the selection of activation functions, and so on, may not be straightforward and can be application specific. The process of fitting multiple models with varying architectures to pinpoint an optimal design can be somewhat ad-hoc and time consuming. Alternatively, Neural Architecture Search (NAS) has proven to be effective in automating the architecture search (Elsken et al.,, 2019). Therefore, identifying appropriate model structure for sea-level change emulators through neural architecture search (NAS) could be helpful in constructing ML-based emulators.

The simple linear regression method of performing uncertainty quantification on top of a nonlinear ML model shares elements from the weather forecasting literature (Glahn and Lowry,, 1972; Glahn et al.,, 2009; Lovegrove and Siegert,, 2023) in addition to split conformal inference (Lei et al.,, 2018). Some advantages of the approach are that that it is relatively inexpensive to implement after the initial (nonlinear) ML model has been fit and is agnostic to the particular ML model and loss function used; it is therefore recommended as an additional approach to UQ for ML predictions.

8 Conclusion

We demonstrate the efficacy of ML models, primarily the NN and CVAE, in emulating sea-level change. By employing a postprocessing strategy using simple linear regression, both the CVAE and NN exhibit accurate prediction performance and high-quality prediction intervals. The CVAE and NN models show prediction performance comparable to a baseline GP model in terms of accuracy, while the baseline GP struggles to adequately quantify uncertainty in terms of prediction interval coverage rates.

We have also highlighted the use of a simple linear regression method for deriving statistically valid uncertainty quantification using (nonlinear) ML predictions. This method of calibrating uncertainties is rooted in the weather forecasting literature (Glahn and Lowry,, 1972; Glahn et al.,, 2009; Lovegrove and Siegert,, 2023) and also shares commonalities with split-conformal inference Lei et al., (2018). The method is analogous to Platt scaling (Platt et al.,, 1999), used for calibrating probabilities of discrete, as opposed to continuous, events. Our UQ approach in this work is straightforward to implement, computationally inexpensive, and has yielded well-calibrated prediction intervals or simulated sea-level change that are of a useful length.

Finally, we have demonstrated an application of the NN emulator to a large sample of synthetic ice-sheet states and showed its promise for producing probabilistic projections of RSL at different locations. Further work is necessary to develop a RSL emulator applicable to the whole Earth and that accounts for important lateral variations in Earth structure. However, ML shows promise as a tool to accelerate workflows for producing sea-level projections that are fast, reliable, and with quantifiable uncertainty.

Appendix A Appendix: Using sectors in Antarctica

We examine the prediction performance of ML models using the mean ice-thickness of geographic sectors across Antarctica as inputs, similar to the work of Van Katwyk et al., (2023). The model structure is the same as described in Section 3 except with using the mean thickness value for each of the 18 regions as inputs (Figure 9).

Figure 10 indicates a substantial decrease in prediction accuracy when using sectorized regions. In addition, we also observe that the empirical coverage rates for all models deviate from the nominal rates. For example, at a significance level of $\alpha=0.05$, the empirical coverage rates are 0.913 (NN), 0.913 (CVAE), 0.924 (GP), and 0.918 (RF), respectively, and the average prediction intervals are noticeably wider (e.g., 0.170 (NN), 0.169 (CVAE), 0.181 (GP), and 0.169 (RF) at $\alpha=0.05$, respectively). This suggests that more advanced dimension-reduction techniques should be employed to preserve the information in the input data.

Figure 9: A map of sectors in Antarctica from the ISMIP6 basins. Numbers represent the basin number.

Figure 10: R-squared comparison between ML model predictions and actual simulated sea-level change values using sections in Antarctica. Values are markedly lower compared to when using the full-dimensional input vector.