A Python library for solving ice sheet modeling problems using physics-informed neural networks, PINNICLE v1.0

Predicting the future contributions of the ice sheets to sea-level rise remains a significant challenge due to our limited understanding of key physical processes (e.g., basal friction, ice rheology) and the lack of observations of critical model inputs (e.g., bed topography). Traditional numerical models typically rely on data assimilation methods to estimate these variables by solving inverse problems based on conservation laws of mass, momentum, and energy. However, these methods are not versatile and require extensive code development to incorporate new physics. Moreover, their dependence on data alignment within computational grids hampers their adaptability, especially in the context of sparse data availability in space and time. To address these limitations, we developed PINNICLE (Physics-Informed Neural Networks for Ice and CLimatE), an open-source Python library dedicated to ice sheet modeling. PINNICLE seamlessly integrates observational data and physical laws, facilitating the solution of both forward and inverse problems within a single framework. PINNICLE currently supports a variety of conservation laws, including the Shelfy-Stream Approximation (SSA), MOno-Layer Higher-Order (MOLHO) models, and mass conservation equations, for both time-independent and time-dependent simulations. The library is user-friendly, requiring only the setting of a few hyperparameters for standard modeling tasks, while advanced users can define custom models within the framework. Additionally, PINNICLE is based on the DeepXDE library, which supports widely used machine learning packages such as TensorFlow, PyTorch, and JAX, enabling users to select the backend that best fits their hardware. We describe here the implementation of PINNICLE and showcase this library with examples across the Greenland and Antarctic ice sheets for a range of forward and inverse problems.

1 Introduction

Ice sheet modeling is essential for projecting future sea-level rise and understanding the complex dynamics that drive ice sheet behavior under changing climate conditions (e.g., Larour et al., 2012; Aschwanden et al., 2019; Goelzer et al., 2020; Seroussi et al., 2020). These numerical models provide insights into fundamental glaciological processes, advancing our understanding of ice flow and its response to and interactions with the climate system (e.g., Schoof, 2007; Joughin et al., 2021; Morlighem et al., 2024). However, developing models that accurately capture current ice sheet behavior and mass change remains challenging, further complicated by limited observational data and an incomplete understanding of critical physical processes. Despite decades of development, significant uncertainties persist in ice sheet models, impacting the accuracy of sea-level projections (Robel et al., 2019; Edwards et al., 2019; Aschwanden et al., 2021). A major source of uncertainty in models comes from parameters that are challenging to measure directly in the field, such as the bed topography (Morlighem et al., 2017, 2020); basal conditions (Gagliardini et al., 2007; Morlighem et al., 2013; Wernecke et al., 2023); and material properties of ice, including its rheology and thermal structure (e.g., Furst et al., 2015; Colgan et al., 2021).

Traditionally, these parameters are estimated by solving inverse problems, often framed as partial differential equation (PDE)-constrained optimization problems (MacAyeal, 1993; Morlighem et al., 2010; Goldberg and Sergienko, 2011). These inverse problems are difficult to solve because they require customized numerical algorithms that depend on the specific PDEs governing the system (Goldberg and Heimbach, 2013; Cheng and Lotstedt, 2020) and often demand regularization to avoid overfitting and facilitate convergence due to their inherent ill-posedness (Tikhonov, 1943). As data availability and quality continue to improve, integrating these data into traditional inverse frameworks presents new challenges, both in terms of computational demands and the complexity of the required numerical methods. This increase in data complexity calls for new approaches that can flexibly integrate observational data while remaining computationally efficient.

Recent advances in machine learning offer promising alternatives for handling the shortcomings of traditional ice sheet models, particularly when dealing with sparse or noisy data (Karniadakis et al., 2021). In particular, Physics-Informed Neural Networks (PINNs) have gained attention as a powerful tool for combining physical laws with observational data (Raissi et al., 2019; Karniadakis et al., 2021; Lu et al., 2021). Unlike conventional approaches that require extensive customization for each specific problem, PINNs allow for the flexible integration of various physical constraints, making them highly adaptable to different modeling scenarios. In glaciology, PINNs have been applied to a range of challenging modeling tasks recently, including inferring basal conditions, simulating complex ice flow dynamics, and testing novel hypotheses about physical processes in ice shelves (Riel et al., 2021; Wang et al., 2025; Jouvet and Cordonnier, 2023; Iwasaki and Lai, 2023; Cheng et al., 2024). The incorporation of physical knowledge within the neural network structure has the potential to introduce a regularizing effect, which is particularly beneficial for handling sparse and noisy observational data, supporting mesh-free modeling, and enabling flexible constraints beyond standard boundary conditions (Seo, 2024). This framework provides a balance between model complexity and accuracy, making it a robust tool for advancing ice sheet modeling (Raissi et al., 2020; Cheng et al., 2024).

In this context, we describe here PINNICLE (Physics-Informed Neural Networks for Ice and CLimatE), an open-source Python library for ice sheet modeling designed to facilitate the solution of both forward and inverse problems using PINNs. PINNICLE provides a unified framework to integrate observational data directly with the governing physical laws. Leveraging the DeepXDE library (Lu et al., 2021), PINNICLE supports machine learning platforms like TensorFlow, PyTorch, and JAX, giving users the flexibility to choose the backend that best suits their hardware and computational resources. This paper outlines the methodological framework of PINNICLE and illustrates its capabilities through applications on glaciers in the Greenland and Antarctic ice sheets.

2 Physics-informed neural networks

In recent years, physics-informed machine learning techniques have become increasingly popular in the field of ice sheet modeling (e.g., Riel et al., 2021; Brinkerhoff, 2022; Jouvet and Cordonnier, 2023; Bolibar et al., 2023; He et al., 2023). Ice flow is governed by a set of PDEs derived from fundamental conservation laws. The extent to which these governing equations are satisfied in different machine learning frameworks varies across studies. For example, the approaches described in He et al. (2023), Bolibar et al. (2023), and Koo et al. (2024) employ neural networks as emulators or surrogate models, learning internal relationships derived from numerical solutions of the governing PDEs, with minimal enforcement of physical constraints. In Jouvet and Cordonnier (2023), the neural network also acts as an emulator but incorporates the PDE residual directly into the training loss function, embedding physical principles into the optimization process for a more physically consistent model. Similarly, Riel et al. (2021) combines observational data with physics-inspired constraints, such as smoothness and sign of the basal drag, though without explicitly enforcing physical laws on the model’s outputs.

This study takes an alternative approach by leveraging neural networks to directly learn the physical constraints and relationships governing ice dynamics, aiming for a more comprehensive integration of data and physics. The architecture of PINNICLE follows the same framework as in Raissi et al. (2019), Wang et al. (2025), Iwasaki and Lai (2023), and Cheng et al. (2024), where a neural network is trained under the constraints of data misfit and PDE loss. The inputs of the neural network are the independent variables of the PDEs, which can be the spatial coordinates (x, y, z) and/or the time variable (t), depending on the governing equations. The outputs of the neural network are all the dependent variables of the PDEs. This architecture enables the neural network, through backpropagation, to compute the required spatial and temporal derivatives of these dependent variables for evaluating the PDEs. The loss function of the PINN consists of two parts, data and physical loss, which are both constructed using these output variables. The data loss measures the misfit between the data and the corresponding output variable at the location and time of the data acquired. The physical loss is computed by evaluating the residual of the PDEs on a set of randomly generated collocation points. Then, the training procedure is to minimize the loss function, such that the output of the neural network satisfies the governing PDEs and also matches the data provided. However, it is important to note that this formulation represents an idealized scenario. In practice, the residuals typically do not reach zero, nor does the data misfit, as discussed in Cheng et al. (2024).

3 Physics

PINNICLE is designed to be flexible so that diverse physical laws along with relevant observational data can be easily integrated. In ice sheet modeling, conservation of mass and momentum form the governing equation of ice dynamics. PINNICLE currently supports these laws by implementing widely used stress balance approximations. Users can readily apply or extend the framework to include additional physical constraints as required by specific modeling objectives.

3.1 Conservation of mass

In ice sheet modeling, the large aspect ratio of the horizontal dimensions (x, y) to the vertical dimension (z) leads most dynamics to occur in the x–y plane. Therefore, it is common practice to model the ice sheet using a depth-averaged approach that effectively captures horizontal behavior, simplifying the model by considering a two-dimensional domain (x, y) and focusing on the horizontal ice velocity components, denoted by $\mathbf{u}=(u,v{)}^T$.

Mass conservation is fundamental to simulating ice sheet dynamics, as it controls the variations in ice thickness over time. Since ice behaves as an incompressible fluid, the rate of change in ice thickness, H, is equal to the sum of the flux divergence in the horizontal direction and vertical mass exchange processes:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(uH\right)}{\partial x}}+{\displaystyle \frac{\partial \left(vH\right)}{\partial y}}=a,\end{array}$$

where a represents the net mass balance field from surface and basal processes, such as surface accumulation (e.g., snowfall) and losses like surface or basal melting.

3.2 Conservation of momentum

We describe here two approximations of the conservation of momentum that are widely used in glaciology. First, the Shelfy-Stream Approximation (SSA), provides a simplified form of the full Stokes equations by neglecting vertical shear stresses (MacAyeal, 1989). SSA is an excellent approximation of ice sheet flow in regions of fast sliding, such as ice streams, and floating ice shelves. This reduction allows us to describe the horizontal motion of ice through a system of PDEs, which balances gravitational driving forces with internal stress gradients and basal drag as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{c}\frac{\partial }{\partial x}\left(\mathrm{4}\mathit{\mu }H\frac{\partial u}{\partial x}+\mathrm{2}\mathit{\mu }H\frac{\partial v}{\partial y}\right)+\frac{\partial }{\partial y}\left(\mathit{\mu }H\frac{\partial u}{\partial y}+\mathit{\mu }H\frac{\partial v}{\partial x}\right)-{\mathit{\tau }}_{\mathrm{b}}\frac{u}{\left|\mathbf{u}\right|}={\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial x},\\ \frac{\partial }{\partial x}\left(\mathit{\mu }H\frac{\partial u}{\partial y}+\mathit{\mu }H\frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mathrm{2}\mathit{\mu }H\frac{\partial u}{\partial x}+\mathrm{4}\mathit{\mu }H\frac{\partial v}{\partial y}\right)-{\mathit{\tau }}_{\mathrm{b}}\frac{v}{\left|\mathbf{u}\right|}={\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial y},\end{array}\end{array}$$

where s is the ice surface elevation, ρ_i denotes ice density, g is gravitational acceleration, τ_b is the basal shear stress, and μ denotes the ice viscosity, which follows Glen's flow law (Glen, 1958), reflecting the nonlinear behavior of ice deformation:

$$\begin{array}{}\text{(3)}& \mathit{\mu }={\displaystyle \frac{B}{\mathrm{2}}}{\left({\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial y}}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}},\end{array}$$

where n=3 is the flow-law exponent, and B is a temperature-dependent pre-factor (Cuffey and Paterson, 2010).

To characterize the relationship between basal shear stress and sliding velocity, we use Weertman’s friction law (Weertman, 1957; Fowler, 1981; Cuffey and Paterson, 2010) as an example:

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{\mathrm{b}}=C^{\mathrm{2}}|\mathbf{u}{|}^m,\end{array}$$

where C is a spatially varying friction coefficient, and $m=\mathrm{1}/\mathrm{3}$.

In regions where vertical deformation significantly influences ice flow, such as in the ice sheet interior, we incorporate a MOno-Layer Higher-Order (MOLHO) model into PINNICLE, following the approach in dos Santos et al. (2022). This formulation represents ice velocity u as the sum of a basal component u^b and a shear velocity u^sh, modulated by a normalized depth factor ζ as follows:

$$\begin{array}{}\text{(5)}& \mathbf{u}={\mathbf{u}}^{\mathrm{b}}+{\mathbf{u}}^{\mathrm{sh}}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\end{array}$$

where $\mathit{\zeta }\left(z\right)=\frac{s-z}{H}$ scales with depth with z varying between the ice base b and the ice surface s.

Both basal and shear velocities are defined on a 2D domain, with vertical variations accounted for by the polynomial in ζ. The MOLHO model extends the SSA equations by accounting for vertical deformation, expressed as

$$\begin{array}{}\text{(6)}& \begin{array}{c}\frac{\partial }{\partial x}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}H{\mathit{ϵ}}_{\mathrm{11}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H\mathit{ϵ}{\mathrm{sh}}_{\mathrm{11}}\right)+\frac{\partial }{\partial y}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{sh}}\right)-{\mathit{\tau }}_{\mathrm{b}}\frac{u_{\mathrm{b}}}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}={\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial x},\\ \frac{\partial }{\partial x}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{sh}}\right)+\frac{\partial }{\partial y}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}H{\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{sh}}\right)-{\mathit{\tau }}_{\mathrm{b}}\frac{v_{\mathrm{b}}}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}={\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial y},\\ \frac{\partial }{\partial x}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{11}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}H{\mathit{ϵ}}_{\mathrm{11}}^{\mathrm{sh}}\right)+\frac{\partial }{\partial y}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{sh}}\right)+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}H{u}^{\mathrm{sh}}=\frac{(n+\mathrm{1})}{(n+\mathrm{2})}{\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial x},\\ \frac{\partial }{\partial x}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}H{\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{sh}}\right)+\frac{\partial }{\partial y}\left({\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}H{\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{b}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}H{\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{sh}}\right)+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}H{v}^{\mathrm{sh}}=\frac{(n+\mathrm{1})}{(n+\mathrm{2})}{\mathit{\rho }}_{\mathrm{i}}gH\frac{\partial s}{\partial y},\end{array}\end{array}$$

where ${\mathit{ϵ}}_{ij}^{\mathrm{b}}$ and ${\mathit{ϵ}}_{ij}^{\mathrm{sh}}$ represent the strain rate components for basal and shear velocities, respectively. These components are detailed by

$$\begin{array}{}\text{(7)}& \begin{array}{ccc}{\mathit{ϵ}}_{\mathrm{11}}^{\mathrm{b}}=\mathrm{4}\frac{\partial u^{\mathrm{b}}}{\partial x}+\mathrm{2}\frac{\partial v^{\mathrm{b}}}{\partial y},& {\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{b}}=\frac{\partial u^{\mathrm{b}}}{\partial y}+\frac{\partial v^{\mathrm{b}}}{\partial x},& {\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{b}}=\mathrm{2}\frac{\partial u^{\mathrm{b}}}{\partial x}+\mathrm{4}\frac{\partial v^{\mathrm{b}}}{\partial y},\\ \\ {\mathit{ϵ}}_{\mathrm{11}}^{\mathrm{sh}}=\mathrm{4}\frac{\partial u^{\mathrm{sh}}}{\partial x}+\mathrm{2}\frac{\partial v^{\mathrm{sh}}}{\partial y},& {\mathit{ϵ}}_{\mathrm{12}}^{\mathrm{sh}}=\frac{\partial u^{\mathrm{sh}}}{\partial y}+\frac{\partial v^{\mathrm{sh}}}{\partial x},& {\mathit{ϵ}}_{\mathrm{22}}^{\mathrm{sh}}=\mathrm{2}\frac{\partial u^{\mathrm{sh}}}{\partial x}+\mathrm{4}\frac{\partial v^{\mathrm{sh}}}{\partial y},\end{array}\end{array}$$

and the vertically integrated viscosities are calculated as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{b}{\overset{s}{\int }}\mathit{\mu }\mathrm{d}z,{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{b}{\overset{s}{\int }}\mathit{\mu }\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\mathrm{d}z,\\ & {\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{b}{\overset{s}{\int }}\mathit{\mu }{\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)}^{\mathrm{2}}\mathrm{d}z,{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{b}{\overset{s}{\int }}\mathit{\mu }{\left({\displaystyle \frac{n+\mathrm{1}}{H}}{\mathit{\zeta }}^n\right)}^{\mathrm{2}}\mathrm{d}z,\end{array}\end{array}$$

where μ is the effective viscosity defined as in Eq. (3).

The basal shear stress τ_b is defined using Weertman’s friction law (Weertman, 1957), as in the SSA model. Alternative friction laws exist (Budd et al., 1979; Gagliardini et al., 2007), and PINNICLE can be easily adjusted to support these friction laws.

4 Data

PINNICLE supports a variety of data formats, including CSV,¹ NetCDF, and MATLAB data files. These files can contain scattered data or structured data, such as the model struct from the Ice-sheet and Sea-level System Model (ISSM). The spatial and temporal coordinates of the data can be different between different variables, offering flexibility to accommodate any complex situations. Depending on the underlying problem to solve, users can specify the amount of data needed for training. The selected data contribute to the calculation of the data misfit component within the loss function.

To evaluate the data misfit, PINNICLE provides multiple metrics, integrating both the built-in data misfit functions from the DeepXDE package (Lu et al., 2021) and commonly employed functions from other machine learning libraries. Table 1 summarizes the data misfit functions currently available in PINNICLE. In the table, d represents the predicted solution from PINNICLE, $\widehat{d}$ denotes the “reference solution” from the data, and $\mathit{\epsilon }=\mathrm{2.2204}\times {\mathrm{10}}^{-\mathrm{16}}$ corresponds to the machine epsilon for double precisions. For advanced users, PINNICLE also provides a flexible interface to define custom misfit functions to specific modeling needs.

Table 1 Data misfit functions implemented in PINNICLE.

5 Neural networks

The fundamental neural network architecture implemented in PINNICLE is the fully connected neural network (FNN), which takes spatial and temporal coordinates as input and the dependent variables of the PDEs as output, through one single neural network. This architecture allows the neural network to implicitly capture the internal relationship among the dependent variables, and it has been widely used in many applications (Raissi et al., 2019; Iwasaki and Lai, 2023; Wang et al., 2025; Lu et al., 2021; Karniadakis et al., 2021; Teisberg et al., 2021). However, in some cases, the training process may not converge, as the relationship between some dependent variables may be too complex for the neural network to capture (Cheng et al., 2024). To address these problems when they occur, we also provide a parallel fully connected neural network (PFNN) (Lu et al., 2021). The PFNN uses multiple FNNs for each dependent variable with the same independent variables as input. In this case, each neural network only needs to learn its corresponding dependent variable from the data and the PDE constraints.

By default, we use the hyperbolic tangent function as the activation function for the neural network. To accommodate different physical problems, PINNICLE applies min–max normalization before the input layer and a min–max denormalization after the output layer. These normalization processes are strictly limited to the inputs and outputs of the neural network, while the governing PDE residuals and data misfit terms are all defined and calculated in the International System of Units (SI). Consequently, there is no need for the user to scale the physical equations or data misfit functions themselves.

While glacier ice typically exhibits low-pass filtering behavior, some variables, such as surface elevation, contain high-frequency components. To capture these variations accurately, we implemented Fourier feature transform (FFT) (Tancik et al., 2020; Wang et al., 2021) in PINNICLE. Assuming the dimension of the input variables x is d, the FFT for m features is written as $f\left(\mathbf{x}\right)=\left[\cos \right(\mathbf{Bx}),\sin (\mathbf{Bx}){]}^T$, where $\mathbf{B}\in {\mathbb{R}}^{m\times d}$ is sampled from a Gaussian distribution 𝒩(0,σ); the functions cos (⋅) and sin (⋅) are applied coordinate-wise on the vector Bx. This transformation is applied between the min–max normalization and the input layer of the neural network.

6 Loss function

The loss function in PINNICLE integrates both data misfit and physical constraints. The general form of the loss function is given by

$$\begin{array}{}\text{(9)}& \mathcal{L}={\mathcal{L}}_{\mathrm{d}}+{\mathcal{L}}_{\mathit{\phi }},\end{array}$$

where ℒ_d and ℒ_φ represent the weighted sum of data misfit and PDE residual, respectively. The specific formulation of the loss function depends on the underlying physical problem, the data, and the misfit functions (Table 1) employed. Examples of detailed loss function formulations are provided in Sect. 8, specifically in Eqs. (10), (11), and (12).

Physical variables in glaciological applications span multiple orders of magnitude; for example, ice velocity is typically around 10⁻⁵ m s⁻¹, ice thickness is around 10³ m, and driving stress of a glacier is approximately 10⁵ Pa. These differences can lead to imbalanced contributions in the optimization process, where certain variables dominate the loss function while others are underrepresented. To address this issue, PINNICLE applies weights to each individual term in the loss function, ensuring balanced contributions from different data sources and governing equations. The impact of weighting strategies on ice sheet modeling problems has been examined in previous studies (Iwasaki and Lai, 2023; Cheng et al., 2024). While no theoretical guidelines exist for determining optimal weights, an extensive set of over 15 000 numerical experiments described in Cheng et al. (2024), including an L-curve analysis, demonstrated that the best performance is achieved when the contributions of all terms are weighted to approximately the same order of magnitude (i.e., around order 1). Based on this empirical finding, PINNICLE assigns default weights for each term in the loss function accordingly. The default weights and their corresponding typical values are provided in Table 2. For greater flexibility, users can modify the weights to meet their specific modeling requirements. Examples of loss function construction and weight selection strategies are provided in Sect. 8.

Table 2 Default weights and typical values of physical variables in PINNICLE.

7 Structure of the package

PINNICLE is configured by the user based on a nested dictionary to define the neural network architecture, governing equations, computational domains, data, and other experiment-specific settings. This design ensures flexibility for incorporating new functionalities, allowing the model to evolve alongside future advances in physics-informed machine learning. To ensure reproducibility, all user-defined configurations are saved in a JSON² file, which can be reloaded to replicate simulations with the same parameters.

The core structure of PINNICLE is built upon five key modules: Physics, Neural Network, Data, Domain, and PINN. The interconnected structure of these modules is depicted in Fig. 1. The Physics module constructs the PDE constraints, gathering independent and dependent variables from the governing equations. These variables are unified to establish a well-defined mapping between inputs and outputs for the PINN framework. The module assigns these variables to the Neural Network and Data modules, enabling integration of physics-based constraints into the computational pipeline.

Figure 1 Flowchart illustrating the structure of PINNICLE, with arrows representing the flow of information between modules.

The Neural Network module constructs the architecture for the neural network based on user configurations. It supports several popular machine learning libraries, including TensorFlow (Abadi et al., 2015), PyTorch (Paszke et al., 2019), and JAX (Bradbury et al., 2018). To enable seamless transitions between these libraries without altering the codebase, the module employs the DeepXDE framework (Lu et al., 2021) as its backend. PINNICLE supports Python 3.8 or higher and requires minimal library versions, including TensorFlow 2.11.0, PyTorch 2.5, and JAX 0.4, ensuring compatibility with current machine learning standards.

The Data module manages data integration by loading datasets and assigning them to the model with specific loss functions as described in Sect. 4. Users can control the volume of training data, enabling the PINNICLE to handle forward or inverse problems based on the combination of data and PDEs provided. To prevent overfitting and generalization, PINNICLE employs a random sampling strategy to automatically load data from the data file depending on their types, which also acts as a form of implicit regularization in the optimization problem. For mesh data (e.g., ISSM or NetCDF), it applies uniform random sampling, while scattered data are downsampled using Cartesian grids to maintain spatial coverage and avoid over-clustering before random sampling.

The Domain module defines the computational domain by generating a polygon based on a list of user-defined vertices. Using Hammersley sequence sampling (Wong et al., 1997), it produces quasi-random collocation points within the domain. These points are then utilized during training to evaluate residuals of the governing PDEs.

After configuring all four modules, the PINN module integrates them using the DeepXDE package, which provides a comprehensive framework for compiling and training PINNs. By leveraging DeepXDE's built-in capabilities, PINNICLE streamlines the training process, allowing users to focus on model setup and interpretation rather than low-level implementation details.

8 Example of applications

8.1 Example 1: an SSA inverse problem on Helheim Glacier

In this example, we use PINNICLE to reproduce results of the inverse problem described in Cheng et al. (2024). The problem involves solving the 2D SSA equation outlined in Eq. (2) for the fast-flowing region of Helheim Glacier in Southeast Greenland. We use a 6-layer fully connected neural network, with 20 neurons per layer. A schematic representation of the problem, including the data used for training, is shown in Fig. 2, and the corresponding Python implementation is provided in Listing 1. We configure PINNICLE to randomly sample $N_{\mathbf{u}}=N_H=N_{\mathrm{s}}=\mathrm{4000}$ data points for each input variables (lines 19 to 22), including ice velocities, surface elevation, and ice thickness, along with N_φ=9000 collocation points over the entire basin of Helheim Glacier for evaluating the PDE residuals (lines 12 to 13). The unknown friction coefficient is inferred by specifying "C":None in the data section (line 20), ensuring that only boundary data are used to constrain the inversion. It is important to note that the friction coefficient within the red box in Fig. 2 is not included in the training processes but is retained as the “reference solution” for validation purposes in Fig. 3b and f.

Figure 2 The PINNICLE framework for Example 1: inverse problem on Helheim Glacier. The training data include ice velocity, ice thickness, and surface elevation. The friction coefficient C, highlighted in the red box, represents the numerical solution from ISSM and is used only on the domain boundary as a boundary condition. It also serves as the “reference solution” for comparison. The governing equation in the box labeled Physics is the vector form of the SSA in Eq. (2).

Figure 3 The PINNICLE predictions and the misfit of Example 1. (a–d) The predictions of surface velocity, friction coefficient, ice thickness, and surface elevation. (e–h) The misfit between the “reference solution” in Fig. 2 and the corresponding PINNICLE predictions in (a)–(d).

Listing 1 Python code of Example 1: inverse problem on Helheim Glacier.

Once the setup is complete, PINNICLE assembles the loss function based on these user-defined inputs. In this case, the complete formulation of the loss function is written as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathcal{L}=& \underset{{\mathcal{L}}_{\mathrm{u}}}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_{\mathbf{u}}}{N_{\mathbf{u}}}}\sum _{i=\mathrm{1}}^{N_{\mathbf{u}}}\left({\left({\widehat{u}}_i-u_i\right)}^{\mathrm{2}}+{\left({\widehat{v}}_i-v_i\right)}^{\mathrm{2}}\right)}}\\ & +\underset{{\mathcal{L}}_H}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_H}{N_H}}\sum _{i=\mathrm{1}}^{N_H}{\left({\widehat{H}}_i-H_i\right)}^{\mathrm{2}}}}+\underset{{\mathcal{L}}_{\mathrm{s}}}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{s}}}{N_{\mathrm{s}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}{\left({\widehat{s}}_i-s_i\right)}^{\mathrm{2}}}}\\ & +\underset{{\mathcal{L}}_C}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_C}{N_C}}\sum _{i=\mathrm{1}}^{N_C}{\left({\widehat{C}}_i-C_i\right)}^{\mathrm{2}}}}\\ & +\underset{{\mathcal{L}}_{\mathit{\phi }}}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_{\mathit{\tau }}}{N_{\mathit{\phi }}}}\sum _{i=\mathrm{1}}^{N_{\mathit{\phi }}}{\left|\mathrm{\nabla }\cdot {\mathit{\sigma }}_{\mathrm{SSA}}-{\mathit{\tau }}_{\mathrm{b}}-{\mathit{\rho }}_{\mathrm{i}}gH\mathrm{\nabla }s\right|}^{\mathrm{2}}}},\end{array}\end{array}$$

where the weights γ_(⋅) are the default weights in PINNICLE as shown in Table 2, and N_C=541 is the number of boundary points used for constraining friction coefficient C. The terms ℒ_u, ℒ_H, and ℒ_s are the weighted MSEs of the data misfit for ice velocity $\left|\mathbf{u}\right|$, ice thickness H, and surface elevation s, respectively. The term ℒ_C accounts for the MSE of the boundary condition imposed on the friction coefficient. However, it is possible to treat C as a free parameter in the inversion. In such cases, users can exclude C from the dataset by removing the entry "C":None from the issm["data_size"] configuration (e.g., line 20 in Listing 1). PINNICLE will then automatically omit the ℒ_C term from the loss function in Eq. (10). The PDE residual, denoted as ℒ_φ, is expressed in the vector form of the SSA in Eq. (2) and is weighted according to the driving stress to ensure balanced contributions in the loss function.

After training for 100 000 epochs (line 5 in Listing 1), the solution is presented in Fig. 3. The root-mean-square error (RMSE) for each variable comparing to the “reference solution” is as follows: velocity at 273.86 m a⁻¹, surface elevation at 22.21 m, ice thickness at 29.07 m, and friction coefficient at 1101.96 Pa${}^{\mathrm{1}/\mathrm{2}}$ m${}^{-\mathrm{1}/\mathrm{6}}$ s${}^{\mathrm{1}/\mathrm{6}}$. These results closely align with those reported in Cheng et al. (2024), with similar spatial patterns of errors. However, our errors are approximately twice as large, due to training for only one-tenth of the epochs compared to Cheng et al. (2024). Despite this, the overall consistency demonstrates the capability of PINNICLE to reproduce and approximate inverse problem results effectively.

Listing 2 Python code of Example 2: inferring basal friction and ice rheology for Pine Island Glacier.

8.2 Example 2: simultaneous inference of basal friction and ice rheology for Pine Island Glacier

This second example demonstrates PINNICLE's capability to simultaneously infer two spatially varying parameters: the basal friction coefficient beneath grounded ice and the ice rheology, described by the flow-law pre-factor, within the floating ice shelf. We employ the SSA in Eq. (2) on Pine Island Glacier, Antarctica. To capture ice rheology variations, we incorporate a spatially varying pre-factor B in Eq. (3). The Python implementation of this example is provided in Listing 2. The neural network architecture consists of a 6-layer fully connected network with 40 neurons per layer. To effectively capture high-frequency variations, we apply a Fourier feature transformation with σ=10 and m=30, as described in Sect. 5 (lines 10 to 12 in Listing 2). Given the larger domain compared to Example 1 in Sect. 8.1, we utilize $N_{\mathbf{u}}=N_H=N_{\mathrm{s}}=\mathrm{8000}$ data points and N_φ=18 000 collocation points to ensure comprehensive spatial coverage and resolution. The training dataset comprises surface velocity, surface elevation, and ice thickness data across the entire domain from PIG.mat. To constrain the model, we impose Dirichlet boundary conditions: setting C=0 on the floating ice shelf with N_C=4000 and prescribing $B=\mathrm{1.41}\times {\mathrm{10}}^{\mathrm{8}}$ Pa s${}^{\mathrm{1}/\mathrm{3}}$ for grounded ice with N_B=4000, corresponding to ice at −10 °C (Cuffey and Paterson, 2010). These two boundary conditions are imposed using a separate data file BC.mat, which contains scatter data points of C and B, as well as their corresponding coordinates.

In this case, the loss function is formulated as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathcal{L}& ={\mathcal{L}}_{\mathbf{u}}+{\mathcal{L}}_H+{\mathcal{L}}_{\mathrm{s}}+\underset{{\mathcal{L}}_C}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_C}{N_C}}\sum _{i=\mathrm{1}}^{N_C}{\left({\widehat{C}}_i-C_i\right)}^{\mathrm{2}}}}\\ & +\underset{{\mathcal{L}}_B}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_B}{N_B}}\sum _{i=\mathrm{1}}^{N_B}{\left({\widehat{B}}_i-B_i\right)}^{\mathrm{2}}}}+{\mathcal{L}}_{\mathit{\phi }},\end{array}\end{array}$$

where ℒ_u, ℒ_H, ℒ_s, and ℒ_φ have the same formulation as in Eq. (10). The additional terms ℒ_C and ℒ_B account for the data misfit of the friction coefficient C on the floating ice shelf and the ice rheology pre-factor B on the grounded ice, respectively.

Figure 4 The comparison of the “reference solution”, PINNICLE predictions, and misfit of Example 2 to infer the basal friction coefficient and ice rheology simultaneously for Pine Island Glacier. (a–e) The “reference solution” of surface velocity, friction coefficient, ice thickness, surface elevation, and the temperature-dependent pre-factor. (f–j) The PINNICLE predictions of surface velocity, friction coefficient, ice thickness, surface elevation, and the temperature-dependent pre-factor. (k–o) The misfit between the “reference solution” and the corresponding PINNICLE predictions.

After 1 000 000 training epochs, the results are presented in Fig. 4. The RMSE for each variable, compared to the “reference solution”, is as follows: surface velocity at 173.13 m a⁻¹, surface elevation at 18.28 m, ice thickness at 24.38 m, basal friction coefficient at 1242.05 Pa${}^{\mathrm{1}/\mathrm{2}}$ m${}^{-\mathrm{1}/\mathrm{6}}$ s${}^{\mathrm{1}/\mathrm{6}}$, and ice rheology pre-factor at 3.81×10⁷ Pa s${}^{\mathrm{1}/\mathrm{3}}$. The magnitudes of these errors are consistent with those observed in Example 1. Larger misfits are observed in two areas: the friction coefficient in slow-moving regions and the rheology near the ice front. Despite these localized discrepancies, the overall performance demonstrates PINNICLE's capacity to simultaneously infer multiple variables.

8.3 Example 3: time-dependent forward modeling of Helheim Glacier

In this final example, we demonstrate PINNICLE's capability to solve a transient problem by modeling the evolution of ice thickness over time using the mass transport equation in Eq. (1). Specifically, we simulate the dynamics of Helheim Glacier from 2008 to 2009. The framework for this example is illustrated in Fig. 5, and the corresponding Python implementation is provided in Listing 3. Since the problem is defined as time-dependent (lines 8 to 10), PINNICLE automatically constructs a spatiotemporal domain and sets the input variables of the neural network to include spatial coordinates x, y, and time t.

Figure 5 The PINNICLE framework for Example 3: the forward time-dependent problem. The training data consist of the initial ice thickness, a time series of ice velocity, and apparent mass balance. Notably, the neural network inputs are automatically adapted to incorporate both spatial and temporal variables.

Listing 3 Python code of Example 3: time-dependent problem of Helheim Glacier.

To drive the forward model, we provide a time series of ice velocity and surface mass balance at 0.1-year intervals (N_t=11), with each time slice containing $N_{\mathbf{u}}=N_a=\mathrm{3000}$ data points for each variable. The initial condition of the ice thickness is specified at t=2008, also using N_H=3000 data points. These data are extracted from the transient simulation results of Cheng et al. (2022). Additionally, we randomly select N_φ=10 000 collocation points distributed across the spatiotemporal domain to enforce the governing equation. The loss function is written as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathcal{L}=& \underset{{\mathcal{L}}_{\mathbf{u}}}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_{\mathbf{u}}}{N_{\mathbf{u}}{N}_t}}\sum _{j=\mathrm{1}}^{N_t}\sum _{i=\mathrm{1}}^{N_{\mathbf{u}}}\left({\left({\widehat{u}}_{i,j}-u_{i,j}\right)}^{\mathrm{2}}+{\left({\widehat{v}}_{i,j}-v_{i,j}\right)}^{\mathrm{2}}\right)}}\\ & +\underset{{\mathcal{L}}_H}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_H}{N_H}}\sum _{i=\mathrm{1}}^{N_H}{\left({\widehat{H}}_i-H_i\right)}^{\mathrm{2}}}}+\underset{{\mathcal{L}}_a}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_a}{N_a{N}_t}}\sum _{j=\mathrm{1}}^{N_t}\sum _{i=\mathrm{1}}^{N_a}{\left({\widehat{a}}_{i,j}-a_{i,j}\right)}^{\mathrm{2}}}}\\ & +\underset{{\mathcal{L}}_{\mathit{\phi }}}{\underbrace{{\displaystyle \frac{{\mathit{\gamma }}_{H/t}}{N_{\mathit{\phi }}}}\sum _{i=\mathrm{1}}^{N_{\mathit{\phi }}}{\left|{\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{v}H\right)-a\right|}^{\mathrm{2}}}},\end{array}\end{array}$$

where ℒ_u and ℒ_a represent the weighted data misfit functions of ice velocity and mass balance across the entire spatial and temporal domain, ℒ_H accounts for the misfit in the initial ice thickness, and ℒ_φ corresponds to the PDE residual of the mass conservation in Eq. (1).

After training for 800 000 epochs using a fully connected neural network with 6 layers and 32 neurons per layer, the results at the initial time step (t=2008) are shown in Fig. 6, while those at the final time step (t=2009) are presented in Fig. 7. The RMSE values for each variable over the entire simulation period (2008–2009), compared to the “reference solution”, are as follows: surface velocity at 186.18 m a⁻¹, surface mass balance at 0.03 m a⁻¹, and ice thickness at 53.11 m. These errors remain within the same order of magnitude as the other two examples, demonstrating PINNICLE's performance in solving transient problems.

Figure 6 The comparison of the “reference solution”, PINNICLE predictions, and misfits of the transient simulation in Example 3 at the initial time step t=2008. (a–c) The “reference solution” of surface velocity, mass balance, and ice thickness. (d–f) The PINNICLE predictions. (g–i) The misfit between the “reference solution” and the corresponding PINNICLE predictions.

Figure 7 The comparison of the “reference solution”, PINNICLE predictions, and misfits of the transient simulation in Example 3 at the final time step t=2009. (a–c) The “reference solution” of surface velocity, mass balance, and ice thickness. Note that (c) is not exposed to the training. (d–f) The PINNICLE predictions. (g–i) The misfit between the “reference solution” and the corresponding PINNICLE predictions.

9 Performance

We present the performance of the examples in Table 3. All experiments were conducted on the Texas Advanced Computing Center Lonestar6 system, using nodes equipped with NVIDIA A100 GPUs with 40 GB of high-bandwidth memory each. The reported wall time (in hours) represents the average of five identical runs for each numerical experiment. The computational cost is influenced by the number of epochs, neural network parameters, data points, and collocation points, among other factors. It is essential to recognize that problems similar to Examples 1 and 3 typically can be resolved more quickly with traditional numerical models, such as ISSM (Larour et al., 2012), often within just a few minutes on multicore CPUs to achieve comparable accuracy (Cheng et al., 2024, 2022). However, solving complex mixed inverse problems, like the one presented in Example 2, using traditional adjoint methods is nontrivial, since it would include the derivation of new adjoint equations with respect to the two unknown variables, simultaneously. Furthermore, we highlight that PINNICLE's current computational performance has not been fully optimized yet, and rapid advances in hardware technology will help further increase its performance. For example, migrating our experiments from NVIDIA V100 GPUs (Cheng et al., 2024) to NVIDIA A100 GPUs has led to at least a twofold reduction in computation time.

Table 3 Performance of the examples in PINNICLE.

Additionally, direct comparisons between machine learning frameworks (on GPU) and traditional data assimilation methods (on CPU) remain limited even in the broader climate science community. Lai et al. (2024), for example, points out critical differences between PINNs and conventional ensemble-based data assimilation methods, such as ensemble Kalman filters (Bauer et al., 2015). PINNs offer the distinct advantage of simultaneously addressing forward and inverse problems, significantly reducing computational demands compared to ensemble methods, which require multiple simulation runs. Conversely, ensemble methods, despite being computationally intensive, provide explicit uncertainty quantification through ensemble spread. Relative to classical adjoint methods, PINNs demonstrate enhanced robustness to noisy observational data, effectively resolving inversion problems with smoother, physically consistent solutions (Wang et al., 2025; Cheng et al., 2024).

For the current version of PINNICLE, all spatial, temporal, and parameter derivatives are computed using reverse-mode automatic differentiation (AD) as implemented in the DeepXDE backend. Although it has been shown that using forward-mode AD for spatial derivatives (e.g., with respect to t, x, and y) and reverse-mode AD for network parameters can significantly improve memory and computational efficiency (Cho et al., 2023), the current version of DeepXDE does not yet support this separable AD approach. We are actively following the ongoing development of this feature and plan to integrate it into PINNICLE as soon as it becomes available. We anticipate that the incorporation of this strategy in future versions of PINNICLE will further enhance the scalability and efficiency, particularly for large-scale simulations.

10 Conclusions

In this study, we introduced PINNICLE, a flexible and robust framework designed to solve a wide range of glaciological problems using Physics-Informed Neural Networks. PINNICLE integrates observational data with physical laws, enabling both inverse and forward modeling across diverse spatial and temporal scales. The framework demonstrates consistent performance in inferring spatially and temporally varying parameters. PINNICLE's flexibility allows users to customize neural network architectures, incorporate various types of data, and apply different physical constraints, making it suitable for complex modeling tasks in glaciology. Additionally, we emphasize that the PINNICLE framework is not intended to replace traditional numerical methods for solving standard forward or inverse problems. Instead, it is best viewed as a complementary tool, especially useful in scenarios involving the integration of diverse physical processes and the exploration of innovative scientific concepts. PINNICLE's ability to seamlessly integrate observational data with physical laws within a unified framework enhances flexibility and can significantly simplify the implementation and evaluation of new physical models. Future developments will focus on incorporating more advanced neural network architectures, optimization techniques, and additional physical constraints to further enhance accuracy and computational efficiency. PINNICLE’s modular design ensures adaptability to evolving machine learning methodologies and glaciological challenges, positioning it as a valuable tool for the cryosphere community.