Calibrating calving parameterizations using graph neural network emulators: application to Helheim Glacier, East Greenland

Calving is responsible for the retreat, acceleration, and thinning of numerous tidewater glaciers in Greenland. An accurate representation of this process in ice sheet numerical models is critical to better predict the future response of the ice sheet to climate change. While traditional numerical models have been used to simulate ice dynamics and calving under specific parameterized conditions, the computational demand of these models makes it difficult to efficiently fine-tune these parameterizations, adding to the overall uncertainty in future sea level rise. In this study, we adopt three standard graph neural network (GNN) architectures, including graph convolutional network, graph attention network, and equivariant graph convolutional network (EGCN), to develop surrogate models for finite-element simulations from the Ice-sheet and Sea-level System Model. GNNs are particularly well-suited for this problem as they naturally capture the representation of unstructured meshes used by finite-element models. When these GNNs are trained with numerical simulations of Helheim Glacier, Greenland, for different calving stress thresholds, they successfully reproduce the observed evolution of ice velocity, ice thickness, and ice front migration between 2007 and 2020. Moreover, these emulators exhibit uncertainties of less than 10 %–20 % when extrapolating to out-of-sample calving parameterization cases. Among the three GNN architectures, EGCN outperforms the others by preserving the equivariance of graph structures. By leveraging the GPU-based GNN emulators, which are 30–34 times faster than traditional numerical simulations, we determine the temporal variations of the optimal calving threshold that minimizes the misfit between modeled and observed ice fronts. This fine-tuned calving parameterization, enabled by GNN emulators, can enhance the reliability of numerical models in capturing glacier mass loss driven by calving.

1 Introduction

Over the past three decades, the Greenland ice sheet has experienced an average annual loss of 170 Gt of ice, resulting in a global mean sea level rise exceeding 15 mm (Otosaka et al., 2023). This trend of mass loss has intensified in recent years. Between 1990 and 2000, the annual mass loss hovered around 40 Gt a⁻¹, but in the 2010s, it surged to approximately 280 Gt a⁻¹ (Otosaka et al., 2023; Mouginot et al., 2019). This escalating mass loss can be attributed to two primary processes: (1) the change in surface mass balance driven by enhanced surface melt and (2) calving and submarine melting of marine-terminating glaciers, commonly referred to as frontal ablation (King et al., 2020; Choi et al., 2021). In specific regions and seasons, ice discharge may be responsible for more than 50 % of the total ice sheet mass loss (Mouginot et al., 2019; King et al., 2020; Choi et al., 2021; Aschwanden et al., 2019).

Enhanced calving significantly impacts ice dynamics, often resulting in ice flow acceleration and subsequent thinning (Bondzio et al., 2017; Cheng et al., 2022; Lippert et al., 2024). Given that calving is sensitive to climate conditions (Greene et al., 2024; Wood et al., 2021), it is important to understand how future calving rates would impact ice sheet mass balance and sea-level rise (Choi et al., 2021). Numerous studies have utilized ice sheet numerical models to explore the applicability of various calving laws to Greenland and Antarctic ice sheets and identified the set of model parameters that minimize the misfit between model and observations (Wilner et al., 2023; Choi et al., 2018). According to Choi et al. (2018), the von Mises calving law (Morlighem et al., 2016, VM;) best replicates observed terminus positions of nine outlet glaciers in Greenland compared to other existing laws.

Although numerical models can provide reliable solutions for ice flow when ice extent is kept constant, capturing the precise impacts of spatiotemporally varying calving rates on terminus migration and ice flow in numerical models remains challenging. Additionally, assimilating remote sensing observations into numerical models, required to infer certain model parameters, is both complex to implement and computationally intensive (Choi et al., 2023). Furthermore, the integration of calving as a boundary condition in numerical models introduces significant complexity since calving directly alters the ice geometry and model domain during simulations. Consequently, identifying the optimal calving parameterizations consistent with observations is difficult and time-consuming, thereby limiting our ability to project the future mass balance of ice sheets under various parameter settings (Choi et al., 2018; Edwards et al., 2021; Morlighem et al., 2020).

To address the computational demands of numerical models, various statistical approaches have emerged relying on faster machine learning models in lieu of ice sheet numerical models. While traditional numerical models often necessitate using high-performance computing clusters to solve partial differential equations (PDEs) on central processing units (CPUs), machine learning emulators offer the advantage of operating on lighter computational resources, leveraging the parallel processing capabilities of graphic processing units (GPUs). For instance, Downs et al. (2023) used a Gaussian process (GP) emulator to infer the sensitivity of time-independent model parameters to the frontal ablation of Helheim Glacier in Southeast Greenland. Although they were able to identify the best set of calving threshold parameters in the VM calving law, their GP approach did not account for spatial relationships or interactions between neighboring nodes. Additionally, their emulator focused on matching the observed and modeled terminus positions along a central flow line rather than the entire glacier system.

Given that ice dynamics and calving are affected by the glacier geometry and underlying bed topography, it is important for the emulators to learn the spatial context across the entire glacier domain. To account for spatial relationships between nodes for emulating ice sheet dynamics, graph neural networks (GNNs) have emerged as an effective neural network architecture for handling irregular non-Euclidean data structures such as molecular structures, point clouds, social networks, and natural language (Zhang et al., 2019). GNNs are adaptable to any type of data structure organized as graphs, comprising nodes (i.e., data points) and edges (i.e., the connections between nodes), which make them suitable to manifest meshes of numerical simulations (Pfaff et al., 2021). GNNs make predictions by utilizing pairwise message passing between nodes, wherein information exchange occurs, updating individual node features through interactions with connected nodes. Inspired by the resemblance of the mesh structure in finite-element analysis to a graph structure, numerous studies have investigated the training of emulators of finite-element numerical simulations using GNNs (Fu et al., 2023; Shivaditya et al., 2022; Black and Najafi, 2022; Perera et al., 2022; Salehi and Giannacopoulos, 2022; Maurizi et al., 2022; Jiang and Chen, 2023).

In this study, we investigate whether GNNs can be used as the backbone architecture for statistical mapping between input and output variables of finite-element ice sheet modeling while embedding the spatial connections between nodes. GNN emulators take direct advantage of unstructured meshes of the Ice-sheet and Sea-level System Model (ISSM), allowing flexible spatial resolution and efficient allocation of computational resources, and can speed up model parameter search. Here, as an illustration, we use GNN emulators to determine how calving parameterization affects ice dynamics while maintaining the accuracy and precision of ISSM simulations. Our focus is on Helheim Glacier, chosen as the target site for training and evaluating our GNN emulators (Fig. 1a). Helheim Glacier is one of the largest outlet glaciers in Greenland, and it has been shown that its ice velocity is closely correlated to the position of its terminus (Cheng et al., 2022). Despite considerable advancements in observation and modeling techniques over the past decades, the mechanisms governing Helheim Glacier's ice front position remain elusive (Bevan et al., 2015; Miles et al., 2016; Cheng et al., 2022). A key contributing factor to this knowledge gap is the computational intensity of numerical models incorporating calving, impeding the fine-tuning of calving parameters to accurately reflect observations. As the main parameter that determines the terminus positions, the temporal variations of calving parameters and their impacts on calving should be explored (Downs et al., 2023). In this study, we train GNN models using simulation data derived from numerical models and evaluate their fidelity and computational efficiency in modeling the dynamics and calving front migration of Helheim Glacier. We assess the potential of GNN architectures as statistical emulators for numerical finite-element ice sheet models to represent spatial features of ice sheet dynamics and calving across the entire glacier domain.

Figure 1 (a) Location of Helheim Glacier, Greenland; (b) ice velocity; (c) bed elevation; and (d) ice thickness near the calving front of Helheim Glacier in 2007. The background image is the MODIS true color image from 20 June 2023.

The remainder of the paper is structured as follows. Section 2 provides a comprehensive review of relevant literature about calving parameterization in ice sheet models, other machine-learning emulators employed in ice sheet modeling, and the use of GNNs as emulators for finite-element models. In Sect. 3, we describe the training data collected from numerical simulations, while Sect. 4 outlines the specifics of GNN architectures. Section 5 presents the accuracy and computational efficiency of these models in replicating the dynamics of Helheim Glacier. In addition, we demonstrate the utility of GNN emulators in optimizing calving parameterization to match simulation results with observations. Finally, we discuss the scientific implications of calving calibration and the limitations of GNN emulators in Sect. 6.

2 Background

2.1 von Mises calving law

In the VM calving law, the calving rate, c, is assumed to be proportional to both the tensile stress and the magnitude of ice velocity, formulated as follows (Morlighem et al., 2016):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }c=\left|\right|\mathit{v}\left|\right|{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}{{\mathit{\sigma }}_{\text{max}}}},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is a scalar quantity representing the effective tensile stress of the ice, σ_max is a stress threshold that needs to be calibrated on a glacier-by-glacier basis, and v is the ice flow velocity at the ice front. The migration rate of the ice front is then determined from the calving rate c using the following equation (Morlighem et al., 2016):

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathit{v}}_{\mathrm{f}}=\mathit{v}-(c+\dot{M})\mathit{n},\end{array}$$

where v_f is the ice front migration rate, $\dot{M}$ is the melting rate on the calving front, and n is a unit normal vector pointing outward from the ice domain. In the ISSM, the numerical ice sheet model we use in this study, moving boundaries of ice sheets are represented implicitly using the level set method (Osher and Sethian, 1988; Bondzio et al., 2016; Morlighem et al., 2016; Cheng et al., 2024). The level set method defines a scalar field φ(x,t): φ is negative in the location of x that is ice-covered at time t, positive where there is no ice, and the zero contour of φ defines the ice boundary. With a given initial condition φ(x,t₀), the level set φ is updated by solving an advection equation using the migration rate v_f that

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\phi }}{\partial t}}+{\mathit{v}}_{\mathrm{f}}\cdot \mathrm{\nabla }\mathit{\phi }=\mathrm{0}.\end{array}$$

In VM, σ_max is the only parameter that needs to be calibrated (Morlighem et al., 2016). A lower value of σ_max correlates to weaker ice and, consequently, a larger calving rate. According to Morlighem et al. (2016), σ_max is consistent with the range of ice tensile strength, which is typically around 1 MPa but may be as low as 0.7 MPa and not exceed 3 MPa. In most of the previous studies, this σ_max value is simply assumed to remain constant over time within a predefined range. However, environmental conditions around the ice front (e.g., temperature, mélange, sea ice, ocean wave, tide) can vary over time, potentially influencing this calving parameter in response and leading to substantial changes in calving patterns (de Juan et al., 2010; Todd and Christoffersen, 2014; Kondo and Sugiyama, 2023; Xie et al., 2019; Wehrlé et al., 2023). Therefore, fine calibration of σ_max is essential for accurately reproducing observed glacier evolution.

2.2 Machine learning emulator for ice sheet modeling

Machine learning techniques have been extensively employed as statistical emulators for numerical ice sheet models and further for calibration of model parameters. Tarasov et al. (2012) applied Markov chain Monte Carlo methods to artificial neural networks (ANNs) for the Bayesian calibration of model parameters. This approach allows for obtaining posterior distributions for different model parameters given observational datasets. Similarly, Brinkerhoff et al. (2021) used an ANN surrogate model to infer the posterior distributions for parameters that govern the behavior of the sliding law and the hydrological model. Downs et al. (2023) employed the GP for the calving parameterization associated with the VM calving law in Helheim Glacier. Nevertheless, the ANNs of these previous studies are simple feedforward networks without spatial embeddings between nodes. Given that ice dynamics are determined by the topographical features, it is essential for the surrogate models to fully capture the interactions between neighboring elements to accurately represent realistic ice sheet behavior.

In replicating spatial characteristics of glaciers, convolutional neural networks (CNNs) have been predominantly adopted as the main architectures based on their advantages in capturing the spatial variations in topographical features of glaciers (Jouvet et al., 2022; Jouvet and Cordonnier, 2023; Jouvet, 2023; Verjans and Robel, 2024). For instance, the CNN of the Instructed Glacier Model (IGM) (Jouvet et al., 2022) reproduced the ice dynamics from the Parallel Ice Sheet Model (Winkelmann et al., 2011, PISM) and CfsFlow models (Jouvet et al., 2008). Jouvet (2023) extended this CNN emulator to address inversion problems, inferring optimal ice thickness distribution, ice flow velocity, and ice surface elevation to match both a Stokes model and observational data. Another CNN emulator introduced by Jouvet and Cordonnier (2023) employed a physics-informed loss function to minimize the energy associated with ice-flow equations during training. However, CNNs cannot take full advantage of finite-element ice sheet models on their native grid because they rely on regular grids (Zhang et al., 2019). Given the GNN's capacity to account for dynamic interactions between nodes and edges (Satorras et al., 2022), it can be a promising tool for predicting the dynamic behavior of ice. Indeed, Koo and Rahnemoonfar (2025) demonstrated the fidelity and computational efficiency of a GNN architecture in emulating a finite-element ice sheet numerical model. This capability of GNN in preserving the mesh structure from the numerical models can have significant advantages in accurately describing ice front migration.

2.3 Graph neural networks for finite-element analysis

GNNs have been broadly used to emulate finite-element numerical models due to the similarities between computational meshes and graphs. Perera et al. (2022) developed a GNN framework to simulate fracture and stress evolution in brittle materials, training their model with data generated from a finite-element method fracture solver. Similarly, Shivaditya et al. (2022) proposed a GNN surrogate model for finite-element simulations of metal forging processes, demonstrating superior performance compared to other machine learning models and achieving a 10-fold reduction in processing time. Salehi and Giannacopoulos (2022) also developed PhysGNN, a GNN framework tailored for simulating soft tissue deformation, and Maurizi et al. (2022) utilized GNNs to predict stress, strain, and deformation across various material systems, including fiber and stratified composites, and lattice metamaterials. Fu et al. (2023) proposed a boundary-oriented graph embedding (BOGE) approach within the GNN framework for solving finite-element cantilever beam problems, incorporating both boundary elements and local neighbor elements. Jiang and Chen (2023) introduced a novel graph attribute representation for triangular meshes in finite-element von Mises stress problems, effectively capturing geometry and boundary conditions to mitigate over-smoothing issues associated with deep GNNs. Black and Najafi (2022) introduced a multi-fidelity GNN for the cantilever beam problem, leveraging low-fidelity projections to inform high-fidelity modeling across arbitrary subdomains of subgraphs. However, despite the prevalence of finite-element analysis in ice sheet modeling (Larour et al., 2012; dos Santos et al., 2021; Gagliardini et al., 2013), to the best of our knowledge, GNNs have yet to be adopted exclusively for modeling calving dynamics and calibrating parameters involved.

3 Data

3.1 Ice sheet numerical simulation

To generate training datasets for the GNN emulators, we conduct transient simulations of ice dynamics and calving of Helheim Glacier between 2007 and 2020, using ISSM (Larour et al., 2012). The shelfy-stream approximation (SSA; MacAyeal, 1989) is used for describing ice flow. The SSA, which assumes depth-independent horizontal velocity and negligible vertical shear stresses, is appropriate for fast-flowing glaciers controlled by basal sliding such as the Helheim Glacier (Cheng et al., 2022; Choi et al., 2018).

The model setup is identical to the one described in Cheng et al. (2022). A two-dimensional unstructured mesh is constructed with a spatial resolution ranging from 100 m in the fast-flowing ice front to 1500 m in the inland domain, ultimately comprising 46 434 elements and 23 466 vertices (nodes). The transient simulations run forward in time with a time step of 1.825 d (0.005 years), and we output the state of the model every 10 time steps (∼18 d). Consequently, each transient simulation generates a total of 261 outputs between 2007 and 2020. Basal friction is calibrated and held constant using the surface velocities from satellite interferometry (Mouginot et al., 2017, 2019) (Fig. 1b). Bed topography and the initial ice thickness are from BedMachine Greenland v6 (Morlighem et al., 2017) (Fig. 1c and d). Surface mass balance (SMB) is from the Regional Atmosphere Model (Tedesco and Fettweis, 2020), and the ocean thermal forcing is from Wood et al. (2021). The melting rate at the calving front (i.e., $\dot{M}$ in Eq. 2) is parameterized based on Rignot et al. (2016). To examine the sensitivity of ice dynamics to σ_max of the VM calving law, we run transient solutions for 9 different σ_max values (i.e., 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.05, and 1.10 MPa) based on the values proposed by Choi et al. (2018).

The ISSM simulations provide the solutions of ice velocity, ice thickness, and a mask of ice-covered region every 10 time steps. We convert the triangular mesh from ISSM into a data structure that aligns with the input and output requirements of deep learning architectures. Specifically, for GNN architectures, we convert the meshes into graph nodes and edges by extracting adjacency matrices that represent the connectivity between nodes. In a triangular mesh, each element consists of three nodes that are interconnected by edges (Fig. 2). Using the nodes and elements of the mesh ensures that the resolution of this graph matches exactly with the finite-element mesh used in ISSM simulations.

Figure 2 Schematic workflow of using GNN emulators to predict calving front.

3.2 Observations

To determine the best σ_max that aligns with real observations, we compare the numerical simulation results with remote-sensing derived terminus positions and ice velocities for the same periods as the numerical simulations (2007–2020). First, we use surface ice velocities with a spatial resolution of 150 m (Mouginot et al., 2017, 2019), which have a precision better than 20 m a⁻¹ (Mouginot et al., 2017). Second, we use a time series of calving front positions of Helheim Glacier from Greene et al. (2024). During the targeted time frame from 2007 to 2020, we use monthly averaged ice front positions, yielding a total of 156 distinct ice front positions for analysis.

4 Method

Ice sheet modeling can be regarded as a node-regression problem within graph structures, where the output features of individual nodes are derived from the input features of nodes. The unstructured meshes of ISSM can be represented as graph structures, with node connectivity expressed via adjacency matrices. Based on the graph structures of the ISSM meshes, we develop three GNN architectures: graph convolutional network (GCN), graph attention network (GAT), and equivariant graph convolutional network (EGCN). Typical GNN architectures update graph nodes iteratively through message-passing processes between neighboring nodes, and the way to achieve this message passing determines the specific type of GNN architecture. For the undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with N nodes u_i∈𝒱, edges $(u_i,u_j)\in \mathcal{E}$, and an adjacency matrix $A\in {\mathbb{R}}^{N\times N}$, lth GNN layer receives a set of node features ${\mathbf{\text{h}}}^{\left(l\right)}=\mathit{\{}{h}_{\mathrm{1}}^{\left(l\right)},h_{\mathrm{2}}^{\left(l\right)},\mathrm{\dots },h_N^{\left(l\right)}\mathit{\}},h_i^{\left(l\right)}\in {\mathbb{R}}^{F_l},$ as the input and produces a new set of node features, ${\mathbf{\text{h}}}^{(l+\mathrm{1})}=\mathit{\{}{h}_{\mathrm{1}}^{(l+\mathrm{1})},h_{\mathrm{2}}^{(l+\mathrm{1})},\mathrm{\dots },h_N^{(l+\mathrm{1})}\mathit{\}},h_i^{(l+\mathrm{1})}\in {\mathbb{R}}^{F_{l+\mathrm{1}}},$ for the next l+1th layer. F_l and F_l+1 are the number of features in each node at lth layer and l+1 layer, respectively. The GCN, GAT, and EGCN operate on graph structures but use different message-passing approaches in updating h^(l+1) from h^(l). By comparing three representative GNN architectures, we evaluate what approach is more effective in replicating ice sheet dynamics and calving from the ISSM simulations.

4.1 Graph convolutional network

First, we employ a GCN proposed by Kipf and Welling (2017). We design a GCN with one input layer, five graph convolutional hidden layers, and one output layer (Fig. 2). The number of hidden layers is determined after conducting trial and error experiments with several options: we tested 1, 2, 5, and 10 hidden layers, and 5 hidden layers showed the best accuracy. The graph convolutional hidden layers are inspired by the localized first-order approximation of spectral graph convolutions on graph-structured data (Kipf and Welling, 2017). For each graph convolutional layer, the number of features is set to 128. Similar to the hidden layers, the number of hidden features was determined from trial and error with 4 options: 32, 64, 128, and 256. The weights of graph convolutional layers are updated via the layer-wise propagation rule as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{h}_i^{(l+\mathrm{1})}=\mathit{\sigma }\left(\sum _{j\in \mathcal{N}\left(i\right)}{\displaystyle \frac{\mathrm{1}}{c_{ij}}}{\mathbf{\text{W}}}^{\left(l\right)}{h}_j^{\left(l\right)}\right),\end{array}$$

where 𝒩(i) is the set of neighbors of node i, c_ij is an appropriately chosen normalization constant for the edge (u_i,u_j) defined as the product of the node degrees (i.e., $c_{ij}=\sqrt{\left|\mathcal{N}\right(j\left)\right|}\sqrt{\left|\mathcal{N}\right(i\left)\right|}$), and ${\mathbf{\text{W}}}^{\left(l\right)}\in {\mathbb{R}}^{F_{l+\mathrm{1}}\times F_l}$. W^(l) is a layer-specific trainable weight matrix (${\mathbf{\text{W}}}^{\left(l\right)}\in {\mathbb{R}}^{F_{l+\mathrm{1}}\times F_l}$), and σ(⋅) is an activation function; we use the Leaky ReLU activation function with a 0.01 negative slope in this study.

4.2 Graph attention network

Since the original GCN filters merely depend on graph structures and node connectivity (Eq. 4), a model trained with a certain graph structure can have limitations in general applicability to different graph structures. The GAT is proposed to address such shortcomings by adding masked self-attention layers (Veličković et al., 2018). This architecture assigns different weights to different nodes in a neighborhood by inserting a self-attention mechanism in a hidden layer, which can allow better generalizability for different ice conditions or topographies. The propagation process in graph attention layers can be expressed by the following equation (Veličković et al., 2018):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{h}_i^{(l+\mathrm{1})}=\mathit{\sigma }\left(\sum _{j\in \mathcal{N}\left(i\right)}{\mathit{\alpha }}_{ij}^{\left(l\right)}{\mathbf{\text{W}}}^{\left(l\right)}{h}_j^{\left(l\right)}\right),\end{array}$$

where α^(l) is the attention score between node i and node j defined as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{ij}^{\left(l\right)}={\text{softmax}}_j\left(e_{ij}^{\left(l\right)}\right)={\displaystyle \frac{\text{exp}\left(e_{ij}^{\left(l\right)}\right)}{{\sum }_{k\in \mathcal{N}\left(i\right)}\text{exp}\left(e_{ik}^{\left(l\right)}\right)}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }{e}_{ij}^{\left(l\right)}=a\left({\mathbf{\text{W}}}^{\left(l\right)}{h}_i^{\left(l\right)},{\mathbf{\text{W}}}^{\left(l\right)}{h}_j^{\left(l\right)}\right),\end{array}$$

where $a:{\mathbb{R}}^{F_{l+\mathrm{1}}}\times {\mathbb{R}}^{F_{l+\mathrm{1}}}\to \mathbb{R}$ is a self-attention mechanism to compute attention coefficient $e_{ij}^{\left(l\right)}$. This attention mechanism a is a single-layer feedforward neural network parameterized by a weight vector $\mathit{a}\in {\mathbb{R}}^{\mathrm{2}{F}_{l+\mathrm{1}}}$, normalized by LeakyReLU function afterward. The graph structure is applied to this attention mechanism by computing $e_{ij}^{\left(l\right)}$ for only nodes j∈𝒩(i), where 𝒩(i) is set of neighbors of node i. We execute three independent attention mechanisms of Eq. (5) and average them for the final graph attention layer (Veličković et al., 2018). Similar to the GCN, the GAT consists of one input layer, five graph attention hidden layers with 128 features, and one output layer (Fig. 2).

4.3 Equivariant graph convolutional network

Another graph neural network we adopt is EGCN, which is designed to conserve equivariance to rotations, translations, reflections, and permutations in a graph structure (Satorras et al., 2022). Since our emulator is intended to predict ice front migration, we anticipate that the preservation of equivariance to rotations and translations on spatial coordinates via the EGCN structure guarantees sufficient generalizability to various graph structures of dynamics systems. An equivariant graph convolutional layer can be expressed by the following equations:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{m}_{ij}={\mathit{\varphi }}_e\left(h_i^{\left(l\right)},h_j^{\left(l\right)},{\left|\left|x_i^{\left(l\right)}-x_j^{\left(l\right)}\right|\right|}^{\mathrm{2}},a_{ij}\right)\text{(8)}& {\displaystyle }& {\displaystyle }{v}_i^{(l+\mathrm{1})}={\mathit{\varphi }}_v\left(h_i^{\left(l\right)}\right)v_i^{\mathrm{0}}+C\sum _{j\ne i}\left(x_i^{\left(l\right)}-x_j^{\left(l\right)}\right){\mathit{\varphi }}_x\left(m_{ij}\right)\text{(9)}& {\displaystyle }& {\displaystyle }{x}_i^{(l+\mathrm{1})}=x_i^{\left(l\right)}+v_i^{(l+\mathrm{1})}\text{(10)}& {\displaystyle }& {\displaystyle }{m}_i=x_i^{\left(l\right)}+C\sum _{j\ne i}{m}_{ij}\text{(11)}& {\displaystyle }& {\displaystyle }{h}_i^{(l+\mathrm{1})}={\mathit{\varphi }}_h\left(h_i^{\left(l\right)},m_i\right),\end{array}$$

where a_ij is the edge attributes, x_i and x_j are the 2D coordinate embeddings for node i and j, respectively, and C is a constant for normalization computed as $\mathrm{1}/\left|\mathcal{N}\right(i\left)\right|$. For the edge attributes, we use three attributes that can be extracted from the connecting nodes: distance, surface slope, and base slope. ϕ_e, ϕ_x, ϕ_v, and ϕ_h are the edge, position, velocity, and node operations, respectively, which are approximated by single-layer neural networks with 128 hidden features. Herein, the x and y components of ice velocities are represented as velocity embeddings (v_i), which induce the displacement of coordinate embeddings (x_i) (Eq. 9), and ice thickness and stress represented as hidden embeddings ($h_i^{\left(l\right)}$). More details about the EGCN architecture are described in Satorras et al. (2022).

4.4 Model training

From the ice sheet simulation results, we collect 2349 graphs (261 sets of results per transient simulations × 9σ_max values). Each graph consists of 23 466 nodes and 139 798 edges. Individual nodes have numerical input and output features from the ISSM simulations and are connected to each other via edges, which are determined by the node connections of the finite elements. The input features include 7 variables that determine the ice sheet dynamic conditions: bed elevation, basal friction, ice thickness at time t−1, ice mask at time t−1, surface mass balance at time t−1, and x and y components of ice velocity at time t−1. The output includes 4 variables that represent ice dynamics at time t: the two components of ice velocity, ice thickness, and tensile stress $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ (Fig. 2). Our emulators are trained to establish a statistical mapping between input and output features: how the input variables at time t−1 change the output variables at time t in numerical modeling. Previous machine learning studies attempted similar types of mapping via various deep learning architecture (Brinkerhoff et al., 2021; Downs et al., 2023; Tarasov et al., 2012; Jouvet and Cordonnier, 2023; Jouvet et al., 2022; Jouvet, 2023). All input and output features are normalized to the range [−1, 1] for stable learning using the nominal maximum and minimum values of these variables. After the GNNs predict ice velocity, ice thickness, and tensile stress $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, the ice mask at t is computed from these outputs using a given σ_max using Eqs. (1)–(3) (Fig. 2).

All 2349 graph structures are divided into training and test datasets based on the σ_max values to assess whether our emulators can be generalized for different σ_max values. We conduct four sets of sensitivity tests to evaluate whether the emulator can be reliably generalized to unseen σ_max values (Table 1). In set A, the training dataset consists of samples with σ_max values of 0.75, 0.90, and 1.05 MPa, while the remaining σ_max values are used for testing. This set evaluates whether the model can represent the full range of σ_max values when trained with evenly sampled σ_max values. In set B, the training dataset includes low σ_max values (0.70, 0.75, and 0.80 MPa) to test whether the model can generalize to high-value samples when trained only with low-value samples. In set C, we use medium σ_max values (0.85, 0.90, and 0.95 MPa) as the training dataset to check whether the model can represent high- or low-value samples when trained only with medium-value samples. Finally, set D contains only high σ_max values (1.00, 1.05, and 1.10 MPa) to check the model performance when extrapolating to lower values. For each sensitivity test, the training and test datasets contain 783 and 1566 graph structures, respectively. The model is optimized using the Adam stochastic gradient descent algorithm with the mean square error (MSE) loss function, over 500 epochs, and a learning rate of 0.001.

Table 1 Division of training/test σ_max in four different sets.

4.5 Model evaluation

We evaluate the ability of our emulators to reproduce ice velocity, ice thickness, and calving front migration by comparing them to ISSM simulation results. For this evaluation, we calculate two metrics: (i) root mean square error (RMSE) and (ii) binary calving front error (BE); RMSE is used to evaluate ice velocity and ice thickness predictions, and BE is used to evaluate the calving front delineation. These metrics are calculated using the following equations:

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }\text{RMSE}(\widehat{y},y)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}\text{(13)}& {\displaystyle }& {\displaystyle }\text{BE}(\widehat{y},y)=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\mathbf{I}({\widehat{y}}_i=y_i),\end{array}$$

where $\widehat{y}$ denotes predicted values, y denotes true values, N is the number of data points, and $\mathbf{I}({\widehat{y}}_i=y_i)$ denotes the indication function that returns 1 if ${\widehat{y}}_i=y_i$ and returns 0 if ${\widehat{y}}_i\ne y_i$. RMSE represents the difference between the prediction and reference; a lower RMSE corresponds to better fidelity. In calculating BE, we convert the ice mask into binary values, representing ice nodes as 1 and non-ice nodes as 0. Therefore, BE is a metric indicating the proportion of incorrectly predicted nodes relative to the total number of nodes; a lower BE corresponds to better fidelity.

5 Results

5.1 Model fidelity

Our GNN models are trained with four different training/test sets, as outlined in Table 1. The training and test errors of ice velocity, ice thickness, and calving front position for the GNN emulators compared to ISSM simulation results are shown in Fig. 3. When the GNN emulators are trained with even-value samples (set A), the training and test errors are almost the same for ice velocity, ice thickness, and calving front. Indeed, t test results confirm no significant difference between training and test errors for all models (Table A1). This result indicates that the emulators trained with even-value samples can effectively represent the full range of σ_max values.

Figure 3 (a) RMSE of ice velocity, (b) RMSE of ice thickness, and (c) binary calving front error (BE) for four different sets. For each set, nine different σ_max datasets are divided into training (three σ_max values) and test (six σ_max values) datasets as described in Table 1. The circle markers indicate the median errors from training datasets, and the crossing markers indicate the median errors from test datasets. The shaded areas indicate the interquartile ranges (i.e., 25 % quantile to 75 % quantile) of errors.

In the case of low-value samples (set B), the test error is occasionally even lower than the training error, particularly in ice velocity prediction, which suggests that the emulators trained only with low-value samples can represent high σ_max value cases with <10 % of uncertainty in most cases (Table A1). Although the test errors are slightly greater than training errors in sets C and D, their differences are primarily within a 5 % uncertainty range for ice velocity and ice thickness and within a 20 % uncertainty range for the calving front (Table A1). The results of sets C and D demonstrate that the GNN emulators trained with medium or high σ_max values can generally have <10 %–20 % of uncertainty in representing the out-of-sample σ_max values. Overall, the sensitivity tests with four different sets show that GNN emulators trained with selected σ_max values can be generalizable for all other σ_max values with <10 %–20 % of uncertainties even if those σ_max values are not included in the training dataset. While the GP-based emulator by Downs et al. (2023) was limited to interpolating within the range of training samples, our GNN emulators can extrapolate to unseen σ_max values. This capability in extrapolation is particularly valuable for optimizing calving parameterization beyond the scope of training σ_max values.

Among the three GNN architectures, the EGCN exhibits notably lower RMSE values: ∼50–60 m a⁻¹ of ice velocity RMSE and ∼30 m of ice thickness RMSE (Fig. 3). That is, the EGCN reduces ice velocity RMSE by 40–150 m a⁻¹ and ice thickness RMSE by 15–20 m compared to the other GNN architectures (see Table A2 for statistical significance). On the other hand, the EGCN and GAT yield similar ice front predictions with BE ranging from 0.60 % to 0.75 %. The relatively lower errors achieved by the EGCN support that its equivariance architecture is beneficial in predicting ice dynamics. As mentioned in Sect. 4, the GCN, GAT, and EGCN have distinct characteristics in message-passing approaches: the GCN simply uses the adjacency status between neighboring nodes to determine the weights during the propagation process (Eq. 4); the GAT uses additional self-attention mechanisms to evaluate the relative importance of neighboring nodes (Eq. 6); the EGCN uses the message passing from all nodes to preserve the equivariance of the entire graph (Eqs. 7–9). We conjecture that the equivariance architecture of the EGCN contributes significantly to the improvement of model fidelity.

Figures 4 and 5 show the maps of ice velocity and thickness on two σ_max test cases of 0.70 and 1.10 MPa predicted by the GNN emulators trained with the training set A. Overall, all GNN emulators replicate the spatial and temporal patterns of ice velocity and thickness successfully. Additionally, they reproduce the trend of ice front retreat at lower σ_max and advance at higher σ_max. In terms of ice velocity (Fig. 4), the EGCN shows the lowest errors of approximately <100 m a⁻¹ across the glacier domain, followed by GAT and GCN. In terms of ice thickness (Fig. 5), the EGCN generally shows an average RMSE <40 m across the glacier domain, but the GCN and GAT have higher RMSEs of around >50 m in some areas.

Figure 4 Maps of ISSM-simulated ice velocity and difference in ice velocity between GNN emulators and ISSM simulations (GCN, GAT, and EGCN from top to bottom) for 0.70 and 1.10 MPa test σ_max values of set A.

Figure 5 Maps of ISSM-simulated ice thickness and difference in ice thickness between GNN emulators and ISSM simulations (GCN, GAT, and EGCN from top to bottom) for 0.70 and 1.10 MPa test σ_max values of set A.

5.2 Computational time

The most significant advantage of GNN emulators is their ability to reduce this computation time by leveraging GPUs while maintaining finite element structures, which is not completely available via other non-graph neural network architectures (Koo and Rahnemoonfar, 2025). When we run the ISSM numerical model, it takes 948 s (∼15.8 min) on average to complete one 13 year transient simulation for a given σ_max value. This computation time of ISSM reflects the total elapsed time on a single node of the Texas Advanced Computing Center (TACC) Frontera supercomputing cluster, which is equipped with 56 cores of Intel 8280 Cascade Lake CPUs (192 GB memory). To evaluate the computational efficiency of the GNN emulators, we also record the time required to generate a single transient simulation for a given σ_max value with our GNN emulators (Table 2). The GNN computations are performed on a CPU (Intel(R) Core(TM) i7-11700F; 32 GB memory) and a GPU (NVIDIA GeForce RTX 3070; 24 GB memory) within the same desktop (Lenovo Legion T5 26IOB6). We observe a dramatic speed-up with GPU-based deep learning emulators, achieving computation times 30–34 times faster than ISSM simulations. The GCN shows the highest speed-up, around 34 times faster than ISSM. Among the three GNN architectures, the GAT shows the longest computational time due to the complexity of its attention mechanisms. The EGCN also requires a longer computational time than the GCN because it operates on all graph nodes to preserve the equivariance in the entire graph, whereas the GCN operates only on the adjacent nodes. It is also noteworthy that using GPUs reduces the computation time of deep learning emulators by up to ∼4 times compared to using CPUs.

Table 2 Computation time for running one 13-year transient simulation from ISSM and GNN emulators and upstreaming computational time for training deep learning models. The shortest elapsed time and training time are highlighted in bold.

This experiment is also promising because, while ISSM computation is performed on a supercomputer, GNNs can be run on a personal desktop. By leveraging GPU-based GNN emulators, computationally intensive simulations of ice dynamics and calving can be quickly reproduced on personal desktops without the need for high-performance computing systems. We expect that these types of emulators will facilitate the efficient parameterization of different ice conditions, ice properties, and external climate forcings.

However, although GNNs can successfully replicate the finite-element structure of ISSM simulations and reduce the computational time by leveraging GPUs, it is important to consider the upstream costs for training deep learning emulators to assess their whole-process efficiencies. Table 2 presents the number of learnable parameters and model training time for the GCN, GAT, and EGCN. This training time is recorded from the 500-epoch training of each model on a multiple-GPU system equipped with 8 NVIDIA RTX A5000 GPUs. The GCN takes the least training time, followed by EGCN and GAT. Herein, these training times are based on the training datasets of set A; however, the training times of four different sets (Table 1) are similar to each other because they have the same number of training samples. It should also be noted that collecting various ISSM numerical simulations with various σ_max values is required to train GNN emulators: in our case, approximately 8500 s (∼2 h 22 min) are required to collect all 9 simulation results with different σ_max values. Nevertheless, given that the emulators show sufficient performances only with 3 training σ_max values, fewer simulations would be sufficient for training. Furthermore, we emphasize that these GNN emulators can be efficiently applied to this target glacier with a short implementation time once they are trained, which can accelerate the search for optimal parameterization settings that align the numerical model with observations.

5.3 Comparison with real observations

We use the EGCN emulator, trained with set A training samples, to determine the appropriate σ_max values that accurately capture the migration of Helheim's calving front, as the EGCN shows the best performances in predicting ice velocity, ice thickness, and calving front positions (Fig. 3). The satellite-derived ice velocity and ice front observations at each time step are used as input for the EGCN emulator to predict the ice front at the next time step. Then, we identify the optimal σ_max value for each time step that minimizes BE of ice front location between prediction and observation (Eq. 13). Given that previous studies (Morlighem et al., 2016; Wilner et al., 2023; Choi et al., 2018) suggest that σ_max should fall between 0.7 to 1.1 MPa, we vary σ_max incrementally by 0.01 MPa within this range and determine the σ_max value with the lowest BE.

Figure 6 shows the temporal variation of optimal σ_max derived from the EGCN spanning from 2007 to 2020. Interestingly, there is a significant temporal correlation between the optimal σ_max and the terminus positions with coefficient >0.5 (p value <0.01). This result indicates that σ_max is an important parameter that determines the calving front migration and should be tuned carefully in calving models. However, we should also note that the temporal variations of optimal σ_max and observed calving front do not agree well during the summers of 2017 and 2019. Given that these two periods coincide with prolonged plume–polynya events lasting more than 20 d in the Helheim Glacier (Melton et al., 2022), we conjecture that buoyant subglacial discharge plumes may have influenced the calving mechanisms during these periods, introducing additional uncertainties to the VM calving law (Everett et al., 2021).

Figure 6 Temporal variations of optimal σ_max derived by EGCN (red dashed line) and the temporal movement of terminus of flow line 1 (a) and flow line 2 (b) in Fig. 1b.

Year-to-year migrations of ice front from observations and numerical models with different σ_max are shown in Fig. 7. It should be noted that a temporarily constant σ_max does not represent the observed calving front migration well. If the σ_max value is set too high (>0.8 MPa), the calving front does not shift or slightly advances after 2017; on the contrary, if σ_max is set too low (≤0.8 MPa), the calving front retreats too much such that the modeling results no longer match the observations. On the other hand, the temporally varying optimal σ_max reproduces the long-term trend of calving front movement along the flow line 1 better than constant σ_max settings, although the short-term variations do not perfectly match (Fig. 8a). However, we should also note that the calving front migration of flow line 2 is not as well predicted as flow line 1 (Fig. 8b). Although some small advance or retreat along flow line 2 is predicted by numerical simulations with the optimal σ_max settings, the calving front remains almost the same for 13 years of simulation. This difference between flow line 1 and 2 suggests the need for spatially varying σ_max parameterization in addition to the temporally varying settings, which is discussed in detail in Sect. 6.

Figure 7 Annual movement of ice front observed by satellite imagery and modeled by ISSM with different σ_max parameterizations: optimal σ_max values from EGCN and constant σ_max values (0.7, 0.8, 0.9, and 1.0 MPa).

Figure 8 Movement of ice front location for (a) flow line 1 and (b) flow line 2 predicted by numerical models with different σ_max parameterizations; optimal σ_max values from EGCN and constant σ_max values (0.7, 0.8, 0.9, 1.0, and 1.1 MPa). The observed ice front movement is drawn as a black solid line.

6 Discussion

In this section, we discuss the scientific implications of σ_max variation and the significant role of GNN emulators in improving the physical modeling of glacier calving. Our results reveal that the fine-tuning of temporally varying σ_max leads to better performance in calving prediction compared to traditionally used constant σ_max parameterization. As a parameter indicating how difficult it is for calving to happen, a greater σ_max means a lower rate of calving, and a lower σ_max means a higher rate of calving. Previous studies have shown that the physical condition of mélange around the ice front, such as ice thickness and concentration, changes the frequency of calving events; a weaker or disappearing mélange area leads to an increase in calving events, which corresponds to lower σ_max (Wehrlé et al., 2023; Xie et al., 2019; Meng et al., 2025). Although this study is focused on finding optimal σ_max values using GNNs, exploring the linkage between optimal σ_max and mélange conditions can contribute to the accurate parameterization of σ_max for calving models. Additionally, given that mélange and sea ice conditions can vary along the ice front, assigning spatially varying σ_max can address the discrepancies in calving front migration for two flow lines shown in Figs. 7 and 8. To set spatially varying σ_max, we can consider finding the best σ_max for multiple flow lines and interpolate these values at the ice front. It would also be valuable to investigate how σ_max interacts with other environmental factors, including tidal forcing and deep water circulation, which are known as significant drivers for calving events (O'Neel et al., 2003; Bassis and Jacobs, 2013; Slater and Straneo, 2022). Consequently, finding the optimal calving parameterization setting and connecting it to external environmental drivers will contribute to the accurate modeling of glacier mass loss and sea level rise beyond the conventional constant calving calibration.

Given the scientific implication of temporally varying σ_max, the significance of this study lies in the (i) fidelity and (ii) computational efficiency of using GNN architectures, especially EGCN, in emulating numerical ice sheet models operating on unstructured meshes to find the optimal calving parameterization. By applying GNNs directly on raw unstructured meshes, we can keep high resolution around the fast ice area to delineate the calving front accurately. In particular, EGCN shows the best fidelity in predicting ice thickness due to its equivariance concept throughout the graph structures. Moreover, the use of GNN emulators results in a dramatic increase in computational efficiency, with speeds 30–34 times faster than the ISSM simulations. This speed-up allows us to find the temporally varying optimal σ_max values quickly. Given that running numerical simulations every time to determine the optimal σ_max setting is extremely time-consuming, our fast GNN emulators reduce the workload associated with calving calibration.

Besides these advantages of GNN architectures, we emphasize that GNN emulators have significant potential for further improvement, particularly through integration with various architectures, including recurrent neural networks (RNNs) (Wu et al., 2021). For example, while our GNN emulator only implements a mapping between a previous time step (t−1) to the next time step (t), embedding recurrent units into the GNN architecture can make it possible to find the sequential relationships between ice sheet dynamics and calving parameters (i.e., mapping between multiple time steps). This recurrent GNN architecture will be able to predict how the historical context of calving parameters affects ice dynamics and ice front migration.

However, it is also worth noting several limitations of our approach. Firstly, despite the computational efficiency of GNN emulators, they should be trained using numerical simulations. Since the model performance is highly dependent on the quality of training datasets, simulation data should be collected carefully with appropriate parameters. Second, the collection of training datasets from various climatological scenarios can be helpful for better generalizability and reliability of emulators to predict future glacier behaviors; however, this process can be exceptionally time-consuming. This generalizability issue can also be associated with the applicability of our GNN emulators for other glaciers beyond Helheim Glacier. To assess the fidelity and transferability of our GNN emulators trained on Helheim Glacier, we apply them to the Pine Island Glacier (PIG), Antarctica, which has different ice velocity and thickness distributions. The ice velocity RMSE in PIG ranges from 200 to 400 m a⁻¹ (2–4 times greater RMSE than Helheim Glacier), and the ice thickness RMSE ranges from 120 to 200 m (5–10 times greater RMSE than Helheim Glacier) (Figs. S1, S2, and Table S1 in the Supplement). These results suggest the importance of training GNN emulators on datasets encompassing a diverse range of ice velocity and thickness distributions to ensure broad generalizability and applicability across different glacier settings. Finally, although we employ the VM calving law to determine calving front migration, it is important to recognize that calving mechanisms are not yet fully understood and may be more complex than those represented by the VM law. While the VM method has been validated for many glaciers in Greenland and Antarctica (Choi et al., 2018; Wilner et al., 2023), detailed calving mechanisms remain elusive. Thus, although our GNN emulators can provide valuable insights into how to select appropriate σ_max for VM calving law, they rely on our current imperfect physical understanding of calving.

7 Conclusions

This study develops three graph neural networks – GCN, GAT, and EGCN – as surrogate models to reproduce finite-element ice dynamics and calving retrieved from the ISSM. After training these GNNs with 13-year transient simulations of Helheim Glacier, they demonstrate significant spatiotemporal agreement with ISSM simulations in predicting ice velocity, ice thickness, and ice front location. The GNN emulators successfully reproduce the retreat of the ice front for a lower calving stress threshold σ_max, as well as the stable condition of the ice front for a higher σ_max. Additionally, the GNN emulators show only <10 %–20 % of uncertainties in predicting out-of-sample σ_max values, implying a better capability in extrapolation and generalizability than previous statistical emulators. Among the three GNN architectures, the EGCN shows the best robustness in the prediction of ice thickness and velocity by preserving the equivariance of graph structures. By leveraging 30–34 times faster computational time of the GPU-based GNN emulators compared to numerical simulations, we efficiently find temporally varying optimal σ_max parameterization for the von Mises calving law. By applying satellite-derived observations to the EGCN, we find that the temporal variations of optimal σ_max have a significant correlation with the retreat and advance of the calving front from 2007 to 2020. Additionally, the numerical simulations based on the fine-tuned σ_max calibration show a better agreement with observations compared to conventional constant σ_max setting. Therefore, this result highlights the importance of setting σ_max values appropriately to improve the reliability of numerical models. As the first attempt to use GNNs for calibrating calving parameters, GNN emulators can further contribute to improving the prediction accuracy of ice sheet mass loss and resulting sea level rise.

Appendix A: T tests for training and test datasets

We conduct statistical t tests to examine (i) whether test error is significantly different from training error (Table A1) and (ii) whether a certain GNN architecture shows significantly better performance than others (Table A2).

Table A1 Reduction in test error relative to training error in percentage (%) for each training/test set and GNN emulator. The positive values mean lower error in test datasets than in training datasets. If the mean difference is not statistically significant (p value >0.01; i.e., statistically they have the same mean), we highlight it in bold.

Table A2 Performance improvement by GAT and EGCN relative to GCN (i.e., error reduction relative to GCN error) in percentage (%). The positive values mean a better performance by GAT or EGCN compared to GCN. If the mean difference between the two models is not statistically significant (p value >0.01; i.e., statistically they have the same mean), we highlight it in bold.