As the mass losses in Greenland and Antarctic ice sheets have accelerated in recent years (Otosaka et al., 2023), it becomes more important to accurately model the ice sheets’ behavior and their interaction with climate change. Although several numerical models have been proposed based on the physical understanding of ice flow dynamics (e.g., Stokes equations), those numerical models require intensive computation to solve complex partial differential equations (PDEs). In order to emulate such computationally expensive

1Department of Computer Science and Engineering, Lehigh University, Bethlehem, PA, USA 2Department of Civil and Environmental Engineering, Lehigh University, Bethlehem, PA, USA. Correspondence to: Maryam Rahnemoonfar <maryam@lehigh.edu>.

ICML 2024 Workshop on Machine Learning for Earth System Modeling, Copyright 2024 by the author(s).

numerical models, computationally cheap machine learning models have been developed leveraging the parallel processing capability of graphic processing units (GPUs).

Most of them relied on convolutional neural network (CNN) architectures (Jouvet et al., 2022; Jouvet & Cordonnier, 2023) that are only compatible with regular Euclidean or grid-like structures, such as images. However, given that finite-element numerical models are implemented on unstructured meshes rather than regular grids, CNN cannot fully represent unstructured meshes of numerical ice sheet models. Therefore, instead of traditional CNN architecture, we choose graph neural network (GNN) as the backbone deep learning architecture to emulate finite-element ice sheet models. Unlike CNN, GNN can operate on any irregular non-Euclidean graph structures (i.e., any data structures with nodes and edges) by updating node features iteratively through message-passing processes between neighboring nodes (Zhang et al., 2019). In particular, considering the dynamical behavior of ice sheets, we adopt a special graph convolutional network (GCN) architecture designed to preserve equivariance to rotation and translations of dynamic systems: so-called equivariant graph convolutional network (EGCN) (Satorras et al., 2022).

In this study, we aim to develop GNN emulators, particularly EGCN, for the Ice-sheet Sea-level System Model (ISSM) (Larour et al., 2012). The computational efficiency of the ISSM numerical model results from adaptive mesh refinement (AMR), which allocates computational resources depending on the expected precision of ice velocity for individual finite elements (Larour et al., 2012; dos Santos et al., 2021). We take the Helheim Glacier, Greenland (Fig. 1a), and Pine Island Glacier (PIG), Antarctica (Fig. 1b), as our testing sites because they are the representative ice sheets that have experienced rapid acceleration and mass loss. Given that mass loss of Helheim Glacier is primarily driven by calving (Choi et al., 2018; Cheng et al., 2022) and PIG is driven by basal melting (Joughin et al., 2021; Jacobs et al., 2011), we predict how the ice sheet dynamics would change by calving and melting parameters.

We train GNN models using the simulation data acquired from the ISSM numerical model and assess their fidelity and computational efficiency for modeling ice dynamics in the Helheim Glacier and PIG. The main contributions of this research are the following:

• We adopt EGCN as a statistical emulator to reproduce the ice dynamics simulated from ISSM. This study is the first application of EGCN on ice sheet modeling.

• We conduct extensive experiments for the Helheim Glacier and PIG to evaluate the potential and superiority of EGCN architecture over traditional GCN and CNN architectures in predicting ice dynamics.

Figure 1. (a) Location of the Helheim Glacier in Greenland; (b) Location of the Pine Island Glacier (PIG) in Antarctica.

The unstructured meshes of ISSM can be regarded as graph structures where the connectivity between nodes is expressed via adjacency matrices. Additionally, since the ice movements change the spatial domains of ice sheets, modeling of ice sheet dynamics can be regarded as a dynamic graph problem. Hence, based on the graph representation of the ISSM meshes, we adopt and test the EGCN, which is specialized for the modeling of dynamically changing graph structures (Satorras et al., 2022). We also test two baseline deep learning models: normal GCN and CNN.

2.1. Equivariant Graph Convolutional Network

EGCN is designed to conserve equivariance to rotations, translations, reflections, and permutations in a graph structure (Satorras et al., 2022): this EGCN structure has shown greater generalizability to any graph structures of dynamics systems. Let the lth graph convolutional layer receives a set of node features has the input and produces a new set of node features, h, for the next l + 1th layer. N is the number of nodes; and is the number of features in each node at lth layer and l + 1 layer, respectively. Then, an equivariant graph convolutional layer updates node features by using the following equations:

where is the edge attributes, are the coordinate embeddings for node i and j, respectively, and C is a constant for normalization computed as 1/|N(i)|. N(i) is the set of neighbors of node i. For the edge attributes, we extract five attributes from the connecting nodes: distance, surface slope, base slope, acceleration of x-component velocity, and acceleration of y-component velocity. , and are the edge, position, and node operations, respectively, which are approximated by single-layer MLPs with 128 hidden features. Herein, we regard the x- and y-component ice velocities as the displacement causing the coordinate changes of graphs (i.e., and ), and ice thickness is represented by hidden features (i.e., ). Thus, the ice thickness is equivariant to the displacement caused by ice flow. Our EGCN consists of one input layer, five equivariant graph convolutional layers, and one output layer (Fig. 2).

Figure 2. Schematic illustration of GCN architectures.

2.2. Graph Convolutional Network (GCN)

We compare the EGCN with a normal GCN architecture proposed by Kipf & Welling (2017). Our GCN consists of one input layer, five graph convolutional hidden layers, and one output layer (Fig. 2). For each graph convolutional layer, the number of features is set to 128. The weights of graph convolutional layers are updated via the layer-wise propagation rule as follows:

where is the product of the square root of node degrees (i.e., ), Wis a layer-specific trainable weight matrix (W), and is an activation function; we use the Leaky ReLU activation function of 0.01 negative slope in this study.

2.3. Convolutional Neural Network (CNN)

As a baseline model to be compared with GNN emulators, we train and test a fully convolutional network (FCN), which has a similar architecture to Jouvet et al. (2022). In the FCN architecture, the hidden graph convolutional layers are replaced with convolutional layers; it consists of five hidden convolutional layers. Each convolution has a 33 kernel size and 128 features. The leaky ReLU activation function of 0.01 negative slope is applied after each convolutional layer. Since the FCN takes regular grids as the input and output, we interpolate all the irregular mesh construction of the ISSM into a 1 km regular grid using bilinear interpolation.

In this study, we apply the EGCN, GCN, and FCN emulators to predict the ice sheet dynamics in the Helheim Glacier varying by calving threshold parameter () of the von Mises (VM) calving law (Morlighem et al., 2016).

3.1. Preparation of training data

To generate training datasets for the GNN emulators, we run transient simulations of ice dynamics in the Helheim Glacier between 2007 and 2020. We use the Shelfy-Stream Approximation (SSA) (MacAyeal, 1989) to explain ice flow in the Helheim Glacier (Cheng et al., 2022; Choi et al., 2018). For the initial condition of the model, we use the surface velocities from satellite observations (Mouginot et al., 2017a; 2019), bed topography from BedMachine Greenland v6 (Morlighem et al., 2017), surface mass balance (SMB) from the Regional Atmosphere Model (Tedesco & Fettweis, 2020), and the ocean thermal forcing from the reanalysis data (Wood et al., 2021). To examine the sensitivity of ice dynamics to the calving parameter , we implement the transient solutions for 7 different values (i.e., 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.0 MPa) based on the values proposed by Choi et al. (2018).

We convert the finite-element simulation results into graph structures by extracting adjacent matrices that represent the connectivity between nodes: for a triangular mesh composed of three nodes, these nodes are connected to each other. Our GNN takes 10 features of graph nodes as inputs, including , time, SMB, initial x- and y-component ice velocities, initial surface elevation, bed elevation, initial ice thickness, and initial ice mask. Then, the output layer predicts 4 features of the graph nodes: x-component ice velocity, y-component ice velocity, ice thickness, and ice mask. All those input and output features are normalized

Table 1. Test accuracy of ice velocity and thickness for deep learning emulators for the Helheim Glacier. All matrices are averaged for 2 testing values: 75 and 95 MPa.

into [-1, 1] for stable learning using the nominal maximum and minimum values of those variables.

A total of 1,827 graphs are divided into the train, validation, and test datasets based on the values to assess if our emulator can be generalized for out-of-sample The data with of 0.70, 0.80, 0.85, 0.90, and 1.0 MPa are used for training and validation: we randomly split them into 70 % and 30 % for training and testing, respectively. The rest of the data with 0.75 and 0.95 MPa are test datasets. Consequently, the number of training, validation, and test datasets is 913, 392, and 522, respectively. The model is optimized by Adam stochastic gradient descent algorithm with the mean square error (MSE) loss function, 400 epochs, and 0.001 learning rate.

3.2. Results

The overall test accuracy of GCN, EGCN, and FCN is shown in Table 1. All deep-learning emulators show outstanding performances in ice velocity and thickness prediction with R-values greater than 0.990 in most cases. EGCN shows the best accuracy among them, with around 10 m/year of lower ice velocity RMSE and 20 m of lower ice thickness RMSE than normal GCN. Whereas GCN uses the spatial distance of neighboring nodes in determining their relative weights in the propagation process (Eq. 5), EGCN uses the message passing from all nodes to preserve the equivariance of the entire graph. The equivariance architecture to the graph rotation and translation could guarantee more generalizability to any graph conditions, leading to the improvement of model fidelity. In particular, EGCN can potentially represent the ice thickness changes induced by ice velocity because the spatial embeddings are used to predict hidden feature embeddings (Eq. 1-4). Additionally, we highlight that GNN emulators outperform traditional FCN. Since FCN uses the interpolated 1-km regular grid for all points without any adjustment to ice velocity, FCN inhibits delineating boundary conditions in a fine resolution, leading to significant errors. Since ISSM and GNN emulators use irregular meshes, which assign a finer resolution for fast ice and a coarser resolution for static ice, they can describe the ice boundaries more accurately with fine resolutions.

We applied the same GCN, EGCN, and FCN architectures to the PIG in Antarctica. Although the same input and output features are used for the PIG emulator, only the calving parameter is replaced with basal melting rate because the ice dynamics in PIG are likely driven by melt-driven thinning near the grounding line rather than calving.

4.1. Preparation of training data

Similar to the Helheim experiment, the training datasets are also collected from the ISSM transient simulation for 20 years using the SSA. We implement the ISSM simulations for different annual basal melting rate scenarios ranging from 0 to 70 m/year for every 2 m/year. Consequently, we collect 8,640 graphs from the experiments with 36 different melting rates. We divide the 8,640 graphs into train, validation, and test datasets based on the melting rate values: melting rates of 10, 30, 50, and 70 m/year are used as validation datasets, melting rates of 0, 20, 40, and 60 m/year as test datasets and the rest are used as training datasets. The total graphs for training, validation, and testing are 6,720, 960, and 960, respectively. We use the following observation data as the initial ice conditions: ice velocity from the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) (Mouginot et al., 2017b), 1 km Antarctic digital elevation model (DEM) (Bamber et al., 2009), bedrock topography data of the Amundsen Sea continental shelf (Nitsche et al., 2007), Antarctic surface temperature data (Comiso, 2000), and Antarctic surface mass balance data (Vaughan et al., 1999).

4.2. Results

Table 2 summarizes the accuracy of ice velocity and thickness prediction from deep learning emulators. All deep learning emulators show significant performance in both ice velocity and thickness predictions with R greater than 0.997. GNN emulators, including EGCN and GCN, show better accuracy than FCN in ice velocity prediction, and EGCN shows the best accuracy. EGCN shows a remarkably low ice thickness RMSE compared to the other emulators, which agrees with the results from the Helheim Glacier. Therefore, the equivariant architecture of EGCN helps improve the predictability of ice thickness changes caused by ice dynamics in the PIG as well.

We record and compare the time to generate the final transient simulations for the Helheim Glacier and PIG experiments (Table 3). The computation time of ISSM is the total elapsed time spent on a single node of the Texas Advanced Computing Center (TACC) Frontera supercomputing clus-

Table 2. Test accuracy of ice velocity and thickness for deep learning emulators for the PIG. All matrices are averaged for 4 testing melting rates: 0, 20, 40, and 60 m/year.

Table 3. Total computational time (in seconds) for producing ice sheet 20-year transient simulations in the PIG for 36 melting rates and training time to train deep learning emulators.

ter, which is equipped with 56 cores of Intel 8280 Cascade Lake CPUs (192 GB memory). The computation times of deep learning emulators are the total elapsed time on a CPU (Intel(R) Core(TM) i7-11700F; 32 GB memory) and GPU (NVIDIA GeForce RTX 3070; 24GB memory). In the Helheim experiment, GCN shows the most dramatic speed-up by around 560 times faster computation time than ISSM, followed by EGCN (260 times speed-up) and FCN (250 times speed-up). In the PIG experiment, GCN and EGCN show 50 times and 44 times faster computation time than ISSM, respectively. FCN takes the longest time to reproduce the 20-year ice sheet dynamics, with only 21 times speed-up compared to the ISSM simulation.

This study introduces an equivariant graph convolutional network (EGCN) as a deep learning emulator to reproduce ice dynamics modeled by the Ice-sheet and Sea-level System Model (ISSM) and compares its performance with traditional graph convolutional network (GCN) and convolutional neural network (CNN). When the EGCN is trained with transient simulations in the Helheim Glacier, Greenland, and Pine Island Glacier (PIG), Antarctica, EGCN shows better accuracy than the others in predicting ice velocity and thickness by preserving the equivariance of graph structures. By reducing computational time 50-500 times compared to the CPU-based numerical models, the EGCN and GCN emulators will be promising tools to investigate the impacts of parameterization on future ice behavior, which will contribute to the improvement of the prediction accuracy of ice sheet mass loss and sea level rise.

Bamber, J. L., Gomez-Dans, J. L., and Griggs, J. A. A new 1 km digital elevation model of the antarctic derived from combined satellite radar and laser data – part 1: Data and methods. The Cryosphere, 3(1):101–111, 2009. doi: 10.5194/tc-3-101-2009. URL https://tc. copernicus.org/articles/3/101/2009/.

Cheng, G., Morlighem, M., Mouginot, J., and Cheng, D. Helheim glacier’s terminus position controls its seasonal and inter-annual ice flow variability. Geophysical Research Letters, 49(5):e2021GL097085, 2022.

Choi, Y., Morlighem, M., Wood, M., and Bondzio, J. H. Comparison of four calving laws to model greenland outlet glaciers. The Cryosphere, 12(12):3735–3746, 2018. doi: 10.5194/tc-12-3735-2018. URL https://tc. copernicus.org/articles/12/3735/2018/.

Comiso, J. C. Variability and trends in antarctic surface tem- peratures from in situ and satellite infrared measurements. Journal of Climate, 13(10):1674 – 1696, 2000. doi: 10.1175/1520-0442(2000)0131674:VATIAS2.0.CO;2. URL https://journals.ametsoc.org/view/ journals/clim/13/10/1520-0442_2000_ 013_1674_vatias_2.0.co_2.xml.

dos Santos, T. D., Morlighem, M., and Seroussi, H. As- sessment of numerical schemes for transient, finite-element ice flow models using issm v4.18. Geoscientific Model Development, 14(5):2545–2573, 2021. doi: 10.5194/gmd-14-2545-2021. URL https://gmd. copernicus.org/articles/14/2545/2021/.

Jacobs, S. S., Jenkins, A., Giulivi, C. F., and Dutrieux, P. Stronger ocean circulation and increased melting under pine island glacier ice shelf. Nature Geoscience, 4(8): 519–523, 2011.

Joughin, I., Shapero, D., Smith, B., Dutrieux, P., and Barham, M. Ice-shelf retreat drives recent pine island glacier speedup. Science Advances, 7(24):eabg3080, 2021. doi: 10.1126/sciadv. abg3080. URL https://www.science.org/ doi/abs/10.1126/sciadv.abg3080.

Jouvet, G. and Cordonnier, G. Ice-flow model emulator based on physics-informed deep learning. Journal of Glaciology, pp. 1–15, 2023. doi: 10.1017/jog.2023.73.

Jouvet, G., Cordonnier, G., Kim, B., L¨uthi, M., Vieli, A., and Aschwanden, A. Deep learning speeds up ice flow modelling by several orders of magnitude. Journal of Glaciology, 68(270):651–664, 2022. doi: 10.1017/jog. 2021.120.

Kipf, T. N. and Welling, M. Semi-supervised classification with graph convolutional networks, 2017.

Larour, E., Seroussi, H., Morlighem, M., and Rig- not, E. Continental scale, high order, high spatial resolution, ice sheet modeling using the ice sheet system model (issm). Journal of Geophysical Research: Earth Surface, 117(F1), 2012. doi: https://doi.org/10.1029/2011JF002140. URL https://agupubs.onlinelibrary.wiley. com/doi/abs/10.1029/2011JF002140.

MacAyeal, D. R. Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica. Journal of Geophysical Research: Solid Earth, 94(B4):4071–4087, 1989. doi: https://doi.org/10.1029/JB094iB04p04071. URL https://agupubs.onlinelibrary.wiley. com/doi/abs/10.1029/JB094iB04p04071.

Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S. Modeling of store gletscher’s calving dynamics, west greenland, in response to ocean thermal forcing. Geophysical Research Letters, 43(6):2659–2666, 2016. doi: https://doi.org/10.1002/2016GL067695. URL https://agupubs.onlinelibrary.wiley. com/doi/abs/10.1002/2016GL067695.

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauch´e, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., No¨el, B. P. Y., O’Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B. Bedmachine v3: Complete bed topography and ocean bathymetry mapping of greenland from multibeam echo sounding combined with mass conservation. Geophysical Research Letters, 44(21):11,051–11,061, 2017. doi: https://doi.org/10.1002/2017GL074954. URL https://agupubs.onlinelibrary.wiley. com/doi/abs/10.1002/2017GL074954.

Mouginot, J., Rignot, E., Scheuchl, B., and Millan, R. Comprehensive annual ice sheet velocity mapping using landsat-8, sentinel-1, and radarsat-2 data. Remote Sensing, 9(4):364, 2017a.

Mouginot, J., Scheuchl, B., and Rignot, E. Measures annual antarctic ice velocity maps, version 1, 2017b. URL https://nsidc.org/data/NSIDC-0720/ versions/1.

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., No¨el, B., Scheuchl, B., and

Wood, M. Forty-six years of greenland ice sheet mass balance from 1972 to 2018. Proceedings of the National Academy of Sciences, 116(19):9239–9244, 2019. doi: 10. 1073/pnas.1904242116. URL https://www.pnas. org/doi/abs/10.1073/pnas.1904242116.

Nitsche, F. O., Jacobs, S. S., Larter, R. D., and Gohl, K. Bathymetry of the amundsen sea continental shelf: Implications for geology, oceanography, and glaciology. Geochemistry, Geophysics, Geosystems, 8(10), 2007. doi: https://doi.org/10.1029/2007GC001694. URL https://agupubs.onlinelibrary.wiley. com/doi/abs/10.1029/2007GC001694.

Otosaka, I. N., Shepherd, A., Ivins, E. R., Schlegel, N.-J., Amory, C., van den Broeke, M. R., Horwath, M., Joughin, I., King, M. D., Krinner, G., Nowicki, S., Payne, A. J., Rignot, E., Scambos, T., Simon, K. M., Smith, B. E., Sørensen, L. S., Velicogna, I., Whitehouse, P. L., A, G., Agosta, C., Ahlstrøm, A. P., Blazquez, A., Colgan, W., Engdahl, M. E., Fettweis, X., Forsberg, R., Gall´ee, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B. C., Harig, C., Helm, V., Khan, S. A., Kittel, C., Konrad, H., Langen, P. L., Lecavalier, B. S., Liang, C.-C., Loomis, B. D., McMillan, M., Melini, D., Mernild, S. H., Mottram, R., Mouginot, J., Nilsson, J., No¨el, B., Pattle, M. E., Peltier, W. R., Pie, N., Roca, M., Sasgen, I., Save, H. V., Seo, K.-W., Scheuchl, B., Schrama, E. J. O., Schr¨oder, L., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T. C., Vishwakarma, B. D., van Wessem, J. M., Wiese, D., van der Wal, W., and Wouters, B. Mass balance of the greenland and antarctic ice sheets from 1992 to 2020. Earth System Science Data, 15(4):1597–1616, 2023. doi: 10.5194/essd-15-1597-2023. URL https://essd. copernicus.org/articles/15/1597/2023/.

Satorras, V. G., Hoogeboom, E., and Welling, M. E(n) equivariant graph neural networks, 2022.

Tedesco, M. and Fettweis, X. Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the greenland ice sheet. The Cryosphere, 14(4):1209–1223, 2020.

Vaughan, D. G., Bamber, J. L., Giovinetto, M., Rus- sell, J., and Cooper, A. P. R. Reassessment of net surface mass balance in antarctica. Journal of Climate, 12(4):933 – 946, 1999. doi: 10.1175/ 1520-0442(1999)0120933:RONSMB2.0.CO;2. URL https://journals.ametsoc.org/view/ journals/clim/12/4/1520-0442_1999_ 012_0933_ronsmb_2.0.co_2.xml.

Wood, M., Rignot, E., Fenty, I., An, L., Bjørk, A., Broeke, M. v. d., Cai, C., Kane, E., Menemenlis, D., Millan, R., Morlighem, M., Mouginot, J., No¨el, B., Scheuchl, B.,

Velicogna, I., Willis, J. K., and Zhang, H. Ocean forcing drives glacier retreat in Greenland. Science Advances, 7(1):eaba7282, 2021. ISSN 2375-2548. doi: 10.1126/ sciadv.aba7282.

Zhang, S., Tong, H., Xu, J., and Maciejewski, R. Graph convolutional networks: a comprehensive review. Computational Social Networks, 6(1):1–23, 2019.