Paper deep dive
Graph neural network for colliding particles with an application to sea ice floe modeling
Ruibiao Zhu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 7/21/2026, 1:58:52 AM
Summary
This paper introduces the Collision-captured Network (CN), a Graph Neural Network (GNN) model designed for simulating sea ice floe dynamics in the Marginal Ice Zone (MIZ). The model utilizes the natural graph structure of sea ice, where nodes represent ice floes and edges represent physical interactions like collisions. By integrating data assimilation techniques and using a one-dimensional framework, the CN model accelerates simulation trajectories compared to traditional Discrete Element Methods (DEM) while maintaining accuracy. The approach addresses the computational intensity of DEM by learning direct mappings from data, offering an efficient tool for climate modeling and forecasting.
Entities (7)
Relation Signals (5)
Collision-captured Network → uses → Graph Neural Network
confidence 95% · This paper introduces a novel approach to sea ice modeling using Graph Neural Networks (GNNs)... the proposed model, termed the Collision-captured Network (CN)
Collision-captured Network → appliedto → Sea Ice Floe
confidence 93% · This paper explores the integration of Graph Neural Networks (GNNs) for modeling discrete sea ice floes in the Marginal Ice Zone (MIZ).
Graph Neural Network → models → Sea Ice Floe
confidence 90% · utilizing the natural graph structure of sea ice, where nodes represent individual ice pieces, and edges model the physical interactions
Collision-captured Network → comparedto → Discrete Element Method
confidence 88% · Traditional numerical methods, while effective, are computationally intensive and less scalable. By utilizing GNNs, the proposed model... accelerates the simulation
Collision-captured Network → inspiredby → Interaction Network
confidence 85% · proposed Collision-captured Network (CN) as Figure 1 coupled with data assimilation (DA) techniques, inspired by the Interaction Network [35].
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper introduces a novel approach to sea ice modeling using Graph Neural Networks (GNNs), utilizing the natural graph structure of sea ice, where nodes represent individual ice pieces, and edges model the physical interactions, including collisions. This concept is developed within a one-dimensional framework as a foundational step. Traditional numerical methods, while effective, are computationally intensive and less scalable. By utilizing GNNs, the proposed model, termed the Collision-captured Network (CN), integrates data assimilation (DA) techniques to effectively learn and predict sea ice dynamics under various conditions. The approach was validated using synthetic data, both with and without observed data points, and it was found that the model accelerates the simulation of trajectories without compromising accuracy. This advancement offers a more efficient tool for forecasting in marginal ice zones (MIZ) and highlights the potential of combining machine learning with data assimilation for more effective and efficient modeling.
Tags
Links
- Source: https://arxiv.org/abs/2602.16213v2
- Canonical: https://arxiv.org/abs/2602.16213v2
Trouble viewing inline? Open PDF directly →
Full Text
101,532 characters extracted from source content.
Expand or collapse full text
Graph neural network for colliding particles with an application to sea ice floe modeling Ruibiao Zhu 1* 1* College of Systems and Society, The Australian National University, ACT, 2600, Australia. Corresponding author(s). E-mail(s): ruibiao.zhu@anu.edu.au; Abstract This paper introduces a novel approach to sea ice modeling using Graph Neural Networks (GNNs), utilizing the natural graph structure of sea ice, where nodes represent individual ice pieces, and edges model the physical interactions, including collisions. This concept is developed within a one- dimensional framework as a foundational step. Traditional numerical methods, while effective, are computationally intensive and less scalable. By utilizing GNNs, the proposed model, termed the Collision-captured Network (CN), integrates data assimilation (DA) techniques to effectively learn and predict sea ice dynamics under various conditions. The approach was validated using synthetic data, both with and without observed data points, and it was found that the model accelerates the simulation of trajectories without compromising accuracy. This advancement offers a more efficient tool for forecasting in marginal ice zones (MIZ) and highlights the potential of combining machine learning with data assimilation for more effective and efficient modeling. Keywords: Collision Simulation, Sea Ice Simulation, Data Assimilation, Graph Neural Network, Machine Learning 1 Introduction This paper explores the integration of Graph Neu- ral Networks (GNNs) for modeling discrete sea ice floes in the Marginal Ice Zone (MIZ). Sea ice acts as a critical regulator of the Earth’s energy balance. Sea ice’s high reflectivity (albedo) plays a crucial role in moderating the global climate by reflecting solar radiation back into space [1–3]. However, the diminishing extent of sea ice due to This preprint has not undergone peer review (when applica- ble) or any post-submission improvements or corrections. The Version of Record of this article is published in Arabian Jour- nal for Science and Engineering, and is available online at https://doi.org/10.1007/s13369-026-11188-z. warming leads to a decrease in albedo and a cor- responding increase in solar energy absorption by the Earth’s surface, accelerating global warming [4]. Thus, the simulation of sea ice floes is a crucial component in climate modeling and prediction, playing a vital role in enhancing understanding of the Earth’s climate system, especially in MIZ [5]. The importance of sea floe simulation in under- standing and addressing the challenges posed by a changing climate cannot be overstated [6–9]. Accurately predicted results from sea floe simu- lation would benefit the economy [10–12], ecosys- tems [13, 14], astronomy [15], and human safety [16]. The Discrete Element Method (DEM) is par- ticularly adept at exploring the complex behav- ioral dynamics and spatial-temporal variations of 1 arXiv:2602.16213v2 [cs.LG] 20 Apr 2026 sea ice, especially on smaller scales such as those found in the MIZ. DEM’s detailed parameteriza- tions allow for the simulation of mechanical inter- actions and responses to external forces like ocean drag and atmosphere drag, providing deep insights into the stress responses and movement dynamics of sea ice under various environmental conditions. Overall, DEM enhances our understanding and ability to predict sea ice behavior in these crit- ical zones [17–21]. One significant limitation of the DEM is its high computational demand [22], which is especially evident in the context of sea ice modeling, where simulations must cover exten- sive geographic areas and long-duration phenom- ena. Historically, constraints on high-performance computing (HPC) resources have favored con- tinuum models, with DEM being too costly to execute [5]. This computational load primarily stems from the necessity to detect contact between neighboring elements, a process that becomes particularly costly when employing polygonal ele- ments due to the complex algorithms required to calculate contact points [23]. Consequently, this restricts the model’s application over larger spatial and temporal scales. In general, solving forward problems typi- cally involves fewer computational resources than solving inverse problems [24, 25]. The GNNs pro- vide a powerful tool for efficiently solving inverse problems by learning direct mappings from data to parameters, and GNNs have proven to be highly effective for analyzing graph-structured data, which underpins their growing application in various domains, including social networks, drug discovery, and more [26–28]. This effectiveness arises from GNNs’ ability to model relationships and interactions directly within the graph data, capturing complex connectivity patterns that are otherwise difficult with traditional neural network approaches [29, 30]. GNNs operate on graphs by utilizing nodes and edges where the nodes represent entities and edges denote the interactions between these enti- ties. Through techniques like message passing, nodes update their states based on the informa- tion from their neighbors, allowing GNNs to learn and generalize from the graph structure efficiently [29, 30]. This process is particularly advantageous for tasks that involve rich relational data. The motivation for applying GNNs to sea ice modeling is supported by other proven successes in simi- larly complex systems in other scientific domains [26–28]. Given these capabilities, the application of GNNs to sea ice modeling presents a promising new avenue. GNNs could capture more insights by conceptualizing sea ice dynamics into graphs, where nodes represent sea ice features, including position and velocity, and edges capture interac- tions such as collisions. Previous efforts to simulate sea ice dynamics within the MIZ have predominantly employed Discrete Element Methods (DEM), focusing on the physical interactions of individual ice floes. These methods, while detailed, are computation- ally intensive and often struggle to scale up due to their high computational demands [5, 22]. Surprisingly, despite the rapid advancement in machine learning technologies [31–34], there has been a noticeable absence of research exploring the use of neural networks for sea ice simulation in the MIZ. This study contributes to this area by introducing a GNN-based sea ice simulation model: proposed Collision-captured Network (CN) as Figure 1 coupled with data assimilation (DA) techniques, inspired by the Interaction Network [35]. By leveraging the inherent graph structure of sea ice interactions in the MIZ, the proposed model offers an innovative solution that reduces computational overhead while enhancing simulation fidelity, providing a promising new direction for future climate modeling efforts. 2 Methods Problem setting This study employs a one-dimensional sea ice simulation as the primary dataset to facilitate a more focused analysis of the collision mechanisms between different sea ice floes. This simplification is strategically chosen to reduce the computational complexity associated with multidimensional sim- ulations while retaining the essential dynamics of sea-ice interactions. By constraining the simula- tion to one dimension, the critical factors that influence floe collisions, such as velocity and posi- tion can be isolated and more effectively captured. In a one-dimensional setting, individual floes do not undergo rotation, eliminating the presence 2 Fig. 1: The illustration of CN predictions with only initial states as inputs. In the visualization, different floes are distinguished by representing each with a uniquely colored circle. The predicted state at time t is denoted with a tilde ( ̃ X t ) over the head, while the ground truth state at time t is represented without the tilde (X t ). The white box represents the input of the model, and blue shaded box stand for the output of the model, and the arrow shows the data flows. The proposed model CN recursively utilizes information from the pre- vious two-time steps to forecast the state at the subsequent time step, predicting the states for all times, except for the initial two-time steps which serve as ground truth inputs. of tangential forces at points of contact. Addi- tionally, external forces such as ocean-induced and atmospheric drag are disregarded. This sim- plification focuses the study on elucidating the mechanisms of collision. Also, disks are used to represent floes in the sea ice modeling to balance geometric simplicity and physical accuracy, enhancing the simulation’s effectiveness. The disk shape simplifies the com- putational process, facilitating more straightfor- ward calculations for interactions such as collisions and contacts, which is essential in scenarios with numerous floes or intricate environmental inter- actions [36]. This form also provides a realistic representation of real ice floe behaviors under var- ious conditions, as evidenced in simulations that analyze the impact of pancake ice floes on cylin- drical structures [37]. Thus, the choice of using disks as a representation for floes strategically bal- ances the need for detailed physical representation with computational efficiency, which proves vital in comprehensive simulations of ice dynamics and interactions. Currently, there is a lack of benchmark datasets specifically designed for training and test- ing neural network models in the context of sea ice simulation. For generating the ground truth data, the DEM and some governing equations [38] are used as shown in the following sections, and the data with visualization shown in Supplemen- tary information F is validated. For example, a floe cannot pass another floe and jump onto the other side of the floe in the one-dimensional simulation setting. Assuming denoting floes from left to right sequentially (Fig.3), the system would maintain x i < x j , where i < j and x i is the x-axis position of floe i at the same time. The equations of motion for sea floes Let x i and v i represent the center position and velocity of the i-th floe, and linear velocities for sea ice floe i are defined as following: v i = dx i dt (1) • v i is linear velocity • x i is the position The right direction is defined as positive and the left direction as negative for force, accelera- tion, and velocity. Thus, x i and v i are treated as scalars in the plane. Translational momentum balance The translational momentum balance based on Newton’s laws of motion F = ma for an ice floe i is m i d 2 x i dt 2 = X j f ij n (2) where x i is the center position of the floe i, and m i is the mass of floe i. When floe j contacts with floe i, f ij n is contact-normal force between floe i and floe j. 3 Fig. 2: The illustration of computing the overlap distance δ ij n when the two floes contact with each other. Contact-normal force with hookean linear elasticity The stresses resulting from axial compressive strain can be modeled as nonlinear [39]. How- ever, a nonlinear model requires extremely small time steps due to the numerical instability. Thus, a common approach for capturing the jamming behavior of the interaction of two floes in sea ice simulation is adopted [18, 38, 40]. The resistive force to contact compression, based on overlap dis- tance δ ij n between floe i and floe j (Fig.2) rather than complete elastic collision, is computed by Hookean linear elasticity [40–42]. This contact force is non-zero when there is an overlapping area between two floes. Otherwise, the contact force is zero. This is a common method in discrete ele- ment simulations to model the contact between two different cylindrical ice floes [38, 41, 43, 44]. When δ ij n =|x i − x j |− (r i + r j ) < 0, f ij n = A ij E ij δ ij n = R ij min(h i ,h j )E ij δ ij n = 2r i r j r i + r j · min(h i ,h j )E ij δ ij n (3) where R ij = 2r i r j r i +r j is the harmonic mean of floe i radius r i and floe j radius r j . E ij is the Young’s modulus, and r i and r j , h i and h j are ice floe i radius and ice floe r radius, ice floe i thickness and ice floe j thickness, respectively [38, 45]. Experiment setup In the simulations, unitless variables for all param- eters and state variables are employed to enhance the generalizability and applicability of the model across various scenarios. This decision allows the model to adapt to different scales and environ- ments without needing specific unit conversions. Also, it is essential to set the initial velocities of the floes to relatively high values and to con- strain their positional domain. These parameters are crucial to ensure that collisions occur fre- quently enough to provide substantial training data (Fig.F4 and Fig.F3), and to facilitate the ver- ification of compliance with physical laws. High initial velocities increase the likelihood and fre- quency of collisions within a controlled spatial domain, thereby enriching the dataset with diverse interaction scenarios. Additionally, by limiting the extent of the positional domain, the accuracy and realism of the floe dynamics can be more easily monitored and evaluated, ensuring that the simu- lation remains within acceptable physical bounds. The Euler method is used in equation 4 com- bining mentioned governing equations 1, 2, and 3 to march into the next time step. In this case, v i t j is the velocity of i at time t j , x i t j is the x-axis posi- tion of i at time t j , and f i t j is the contact force of i at time t j if there is any. v i t j = v i t j−1 + f i t j−1 m i dt x i t j = x i t j−1 + v i t j dt (4) The initial conditions for one dimension set- ting are the following: the radius of each floe is 1; the thickness of each floe is 1; the position of each floe is randomly generated in the domain, the initial velocity is from 150 to 200, and the direction of velocity can be left as negative or right as positive along with x-axis, which can be interpret as approximately 150 to 200 meter per hour if compared to the real observation. This is on the same order of magnitude as observed floe drift in real sea-ice fields in certain areas [46]. The larger number of unitless velocities can make floes traverse the bounded domain more quickly and are more likely to collide with each other or with the domain boundaries, leading to more frequent collision events. This is desirable for a collision-capturing training scenario, as it provides a richer dataset of collision outcomes for model learning. Also, the boundary is 0 and 100, and 0 and 200 for 10 and 30 nodes simulation, respec- tively. Furthermore, sea floes perform the same 4 Fig. 3: The illustration of experiment set- tings. The boundaries are shown as black walls, and the floes are represented as circles. collision mechanical interactions with the bound- ary as they contact with other floes, except that the boundaries are stationary (Fig.3). The time step length and Young’s modulus are carefully selected to accurately simulate the dynamics of sea ice floes. The need for smaller time steps in complex models arises from the necessity to accurately capture the dynamic behaviors and interactions within the collision system. Smaller time steps ensure that the numerical methods can adequately resolve these interactions without missing critical dynamics [47, 48]. For example, larger time steps, such as dt = 10 −2 and dt = 10 −3 , have been tested in the data generation set- ting, and they cannot work since the overlapping area of floes would be larger. This would cause the contact force to be significantly large based on the formula, and the updated x-axis position can be pretty inaccurate based on the extremely large acceleration. For example, a floe can pass to the other side of its neighbor due to the signifi- cantly large contact force, which is impossible in this physical domain. Thus, the model is adjusted to smaller time steps such as dt = 10 −4 for the data generation. For Young’s modulus E, reducing the elastic stiffness is a known technique in DEM simulations to permit larger stable time steps and speed up computations [38, 45, 49]. Using a very high Young’s modulus E would thus yield extremely stiff contacts and large, rapidly vary- ing contact forces, necessitating extremely small time-step sizes to maintain numerical stability. In contrast, a lower Young’s modulus E softens the contacts, which caps the peak contact force and makes collisions less abrupt. Thus, E = 2×10 7 was chosen in this study to reach a reasonable balance between computational efficiency and elastic com- pressibility [38], although this Young’s modulus is lower than the observed value [50, 51]. By optimizing these parameters, the simula- tions realistically capture the mechanical prop- erties and interactions of the ice, and aid in effectively training and validating neural network models. Related models Sea-ice floe dynamics with collisions pose a learning problem that when contacts occur, the state can change abruptly, with post-collision velocities depending on the coupled states of the colliding bodies and contact parameters. This makes the data distribution strongly imbal- anced—dominated by no-contact steps, while the important behavior is concentrated in rare, high- impact events. A closely related collision bench- mark in the Interaction Network (IN) highlights exactly this difficulty. For more than 99% of steps, a ball is not in contact and its next-step velocity equals its current velocity, whereas the remaining steps require learning a complex, state-dependent collision response [35]. Hamiltonian Neural Net- work (HNN), which represents dynamics through differentiable equations and recovers time evolu- tion by taking gradients, is a typical example of a physics-learning machine learning approach that can struggle in such settings, since HNN requires differential equations to successfully describe the whole system within the space and time domains, and these settings cannot satisfy the requirement of HNN [52]. The Interaction Network (IN) treats objects as nodes and relations as edges, and the Graph Network-based Simulator (GNS) extends IN with an encoder-decoder mechanism, noise injection, and stacked graph layers [35, 53]. These models perform well in predicting particle movement tra- jectories. These two models can be directly aligned with sea-ice floe modeling, and these models are tested for comparison. The CN model framework The total time T can be discretized into n time steps with t j = j T n ,j ∈ [0,...,n− 1] . The state of floes at time t j is denoted as X t j . The position and velocity at time t j are denoted as x t j and v t j , respectively. X t j−2 :t j−1 = [x t j−2 ,x t j−1 ,v t j−1 ,r] is previous consecutive objects’ information that contains the second last predicted position x t j−2 , last predicted position x t j−1 , the most recent inferred velocity v t j−1 and radius r. When constructing the graph, the floes and two boundaries can be treated as nodes, and the bidirectional edges are constructed between nodes and their neighbors (including the boundaries) 5 Fig. 4: The illustration of graph construction for GNN. The boundaries are shown as black walls, and e represents the edge, and x stands for the node. (Fig.4). The edge set is defined by neighbor adja- cency rather than distance-based thresholds: edges are created only between designated neighboring nodes, irrespective of their instantaneous separa- tion. With this construction, the graph connectiv- ity remains constant over time, and therefore, the number of edges is fixed throughout the simula- tion. This fixed relational structure provides two practical benefits: it removes the computational overhead and potential instability associated with dynamically rebuilding graphs at every time step, and it ensures that the learned interaction func- tions operate on a consistent topology, so that changes in predicted dynamics arise from evolving node and edge features rather than from time- varying graph structure. The edge feature E t j = [e ij ;e ji ] is the concatenation of the displacement of floe i to floe j and the displacement of floe j to floe i. To identify which node is the sender or the starting point for each directed edge in the graph, R s is defined. The entries of R s are binary with R s [i,j] = 1 if the j-th edge originates from node i and R s [i,j] = 0 if the j-th edge does not originate from node i. Meanwhile, to identify which node is the receiver or the endpoint for each directed edge in the graph, R r is defined. The entries of R r are binary with R r [i,j] = 1 if the j-th edge terminates at node i, and R r [i,j] = 0 if the j-th edge does not terminate at node i. More details regarding relation and permutation matrices are in Supplementary information B. The proposed Collision-captured Network (CN) has been defined as, E f = φ θ 1 ([R T r X t j−2 :t j−1 ;R T s X t j−2 :t j−1 ;E t j−1 ]) X t j = γ θ 2 ([X t j−2 :t j−1 ;R r E f ]) (5) Since the message passing in GNN is simi- lar to multi-layer perception (MLP), φ θ 1 with θ 1 learnable parameters and γ θ 2 with θ 2 learnable parameters are set to be MLP with four 150- length hidden layers and one 100-length hidden layer, respectively [54]. When the simulation has N floes and 2 boundaries, the number of edges is N e = 2(N + 1). Thus, both R r and R s are (N + 2)× N e matrix. E t j−1 is N e -length vector. X t j−2 :t j−1 is (N + 2)× 4 matrix, since it contains 4 different features as mentioned before within it. The φ θ 1 serve as the aggregating function [55, 56] the edge effects by taking the concatenation of interaction terms and edge features. The output of φ θ 1 is (N + 2)× D φ matrix, where the D φ is the output dimension of the φ. Thus, R r E f is the aggregation edge effects working on the nodes with (N + 2)×N e dimension. Furthermore, γ θ 2 serve as updating function [55, 56] in GNN by taking the input of the concatenation of node features and aggregated edge effects working on the nodes. The proposed model differs from the standard IN formulation in three deliberate and practi- cally important ways: (i) the state representation used for prediction, (i) the edge features to bet- ter encode collision-relevant geometry, (i) the activation function. A standard Interaction Network predicts the next-step state by conditioning on the full state at the most recent time step, typically including both position and velocity as node features. This design is appropriate when velocities are directly available from a simulator or can be reliably mea- sured. In contrast, our proposed model is intended to align with more realistic observation conditions in sea-ice applications, where position is the pri- mary quantity that can be consistently observed, while velocity is often unavailable or substantially noisier. Accordingly, we do not require velocity as an explicit input. Instead, we use positions from the two most recent time steps to derive the most recent-step velocity through a finite- difference approximation. This formulation retains essential dynamical information while reducing reliance on velocity observations. A second distinction is the explicit incorpo- ration of collision-relevant geometry into edge features. In standard Interaction Networks, inter- actions are only represented via learned message functions, but edge features are not necessarily constructed. Our model implements the edge rep- resentation with a specialized feature: the instan- taneous inter-floe displacement. This displace- ment signal directly encodes proximity, which is a primary determinant of contact and collision behavior. By injecting the instantaneous inter-floe 6 as edge features, the model becomes more sensi- tive to contact configurations and is better posi- tioned to learn physically plausible responses, such as resolving collision and preventing unrealistic overlap. All activation functions are Mish functions Mish(x) [57] as equation C13 to improve perfor- mance, as compared to all ReLU functions [58] in IN. This Mish activation function can improve the performance of the model since this function is non-monotonic, non-monotonic derivative, not saturated, and it has infinite continuity and an approximating identity near the origin. Unlike the ReLU, the non-monotonic utilizes the negative derivative and improves expressiveness. Compared to the sigmoid function, Mish is not saturated, and it avoids near-zero gradients to make the model learn more efficiently and effectively [57]. The model’s output can be velocity or position at time t j , proving the velocity would be bet- ter based on the error analysis in Supplementary information C. Thus, the loss function for training is: L = 1 n n X j=2 (CN(X t j−2 :t j−1 )− v t j ) 2 (6) To train the proposed model, the training pairs (X t j−2 :t j−1 ,v t j ) are randomly sampled, and be computed above loss function between the pre- dicted velocity and ground truth velocity of each floe. Based on Equation 6, the optimization prob- lem is to minimize the loss with θ 1 and θ 2 . Then, the parameters of the model are optimized over this loss with the Adam optimizer for its opti- mization performance in accelerating convergence [59]. More details about training are shown in the Supplementary information E. A key feature of the proposed model is its parameter-sharing structure, which enables the model to generalize effectively while maintaining statistical efficiency. Parameter sharing means that every edge mes- sage is computed using the same φ θ 1 , and every node state is updated using the same γ θ 2 across all time steps, rather than learning separate func- tions for particular floes or particular interacting pairs, unlike the GNS that employs deep stacks of message-passing blocks, each layer can have its own learnable weights, even though the layer structure is repeated. This principle aligns with the core motivation of IN style models: sharing the relation and object reasoning mechanisms enables the learned dynamics to generalize across different numbers and configurations of entities, because the model learns reusable interaction rules rather than instance-specific behavior [35]. By sharing parameters across tasks or components within the model, the overall number of learnable parameters is reduced, thus mitigating overfitting risks and reducing training time. This approach leverages the benefits of shared information, where parame- ters trained in one context can transfer knowledge to another, enhancing the model’s adaptability across varying data distributions. Consequently, parameter sharing allows the model to make effi- cient use of available data, resulting in faster con- vergence and lower computational costs, without compromising predictive accuracy [60]. In order to initialize the proposed model in prediction, it is imperative that the initial states for the first two consecutive time steps, which are t 1 and t 2 , along with the relevant edge informa- tion R r , be fed into the model as ground truth. This foundational data acts as a crucial starting point, enabling the model to set initial conditions accurately. Once initialized, the model leverages these inputs to predict the object states in the next step. To facilitate subsequent predictions, the model is designed to recursively incorporate its own predictions as inputs. This recursive appli- cation allows the model to continually refine its predictions based on the evolving dynamics cap- tured by the updated object states. Thus, the model functions in a feedback loop, where the out- put from one step serves as the input for the next, enhancing the model’s predictive capability over time. Furthermore, observational data is combined with the proposed model forecasts to refine predic- tions by using data assimilation. Specifically, the Ensemble Kalman filter (EnKF) and the Ensem- ble Transform Kalman filter (ETKF) are used to better handle the uncertainties in the system in the long-time trajectories. More details about both Kalman filters are in the Supplementary information A. 7 3 Results The approach was verified using unseen truth data. The results show that CN not only accel- erates the simulation of sea ice trajectories but also generalizes well to much longer time scales than those included in the training set, predict- ing sea ice dynamics over extended periods with high accuracy. This improvement in efficiency, combined with the model’s ability to predict accu- rately over longer time scales, makes the proposed approach a promising tool for more effective and computationally efficient sea ice modeling, partic- ularly for applications in MIZ forecasts. Model performance Three root mean squared error (RMSE equation from D17) metrics are reported to separate short- horizon accuracy from long-horizon simulation stability. One-step position RMSE is defined as the root mean squared error of the average next- step position prediction when the model is con- ditioned on ground-truth inputs. Similarly, one- step velocity RMSE is the RMSE of the average next-step velocity prediction under ground-truth inputs. Also, simulation RMSE measures long- horizon error: it is computed as the time-average of position RMSE over the entire simulation when the model is input only with the initial ground- truth states and then evolves autoregressively. In addition to RMSE, the pattern correlation coefficient (PCC equation from D18) is used to quantify how well the predicted trajectories pre- serve the ground-truth spatiotemporal pattern; as with standard correlation measures, PCC takes values in [−1, 1], where negative values indicate an inverse (anticorrelated) relationship between predicted and target patterns. Table 1 shows that all CN variants outper- form the Interaction Network (IN) baseline. This improvement is consistent with two design choices. First, CN explicitly augments the interaction model with specialized edge features (including direct inter-floe displacement information), which facilitate the message-passing mechanism to pre- serve separation and maintain physically plausible spacing during close-contact events. Second, CN conditions its updates on two consecutive posi- tion states, using these to infer the most recent velocity internally, rather than requiring exter- nally provided velocities. By contrast, the origi- nal IN formulation predicts dynamics via shared relation- and object-centric functions applied over a graph, but does not intrinsically enforce domain- specific edge geometry. Empirically, the fact that CN with one step already improves over IN indicates that edge-feature implementation and replace ReLU with Mish activation function can shift performance from moderate PCC to high PCC, supporting the interpretation that explicit displacement-aware relational features are partic- ularly beneficial in floes simulation tasks. Across CN variants, Table 1 indicates that Mish yields slightly better performance than SiLU, while both outperform ReLU. This aligns with prior findings that smoother, self-gated non- linearities can improve optimization and represen- tational fidelity in deep networks. Since Mish is explicitly related to the Swish family, the sim- ilarity in performance between Mish and SiLU is expected, whereas ReLU can be slightly less expressive [57, 61, 62]. The Graph Network-based Simulator (GNS) performs poorly in Table 1, with negative PCC values indicating that the learned simulation cap- tures an oppositely varying pattern relative to the target. While GNS is a strong general-purpose learned simulator that evolves particle systems via stacked message-passing blocks, long-horizon performance is known to depend sensitively on architectural depth and the handling of error accu- mulation. In collision floe dynamics, interactions are fundamentally local (short-range proximity), and propagating information beyond immediate neighbors can be unnecessary or even harmful. In message-passing GNNs, k rounds of prop- agation yield representations that increasingly mix information over the k-hop neighborhood, and deeper stacks are therefore more prone to blending local identities [63]. This connects to the well-documented over-smoothing phenomenon in deep GNNs, where node embeddings become indistinguishable as depth increases, degrading performance on tasks that require sharp local distinctions [64, 65]. Additionally, GNS uses a five- time-step information to forecast the next state; however, in discrete time, a two-step position his- tory is already adequate to derive the velocity or position via finite differences [17, 22, 23]. In this 8 sense, conditioning on longer histories may intro- duce redundant or noisy temporal correlations that are not required by the underlying DEM- generated dynamics, potentially contributing to degraded PCC and inflated RMSE in simulation. The results from this study on sea ice sim- ulation, explicitly focusing on one-dimensional scenarios involving 10 and 30 floes collisions, demonstrate significant achievements in the pre- dictive accuracy of this newly developed model. The model’s effectiveness is underscored by high Pattern Correlation Coefficients (PCC equation from D18) and low Root Mean Squared Errors (RMSE equation from D17) across both sce- narios. Notably, all variables in the proposed model are unitless to ensure the generalizability of results across different scenarios and settings (More details in Supplementary information D). Visualization can be checked in Supplementary information F including snapshots for the predic- tion (Fig.F5 and Fig.F6). Further details about methodologies are comprehensively discussed in the Methods section. For the system with 10 floes, the PCC was high at 98.98% (Fig.5a and Table A1), indicat- ing an almost perfect alignment with the observed trends, suggesting that the model can accu- rately simulate the interactions and dynamics of a smaller number of ice floes. Moreover, the RMSE for this scenario was 1.16, which is relatively low considering the domain range of 100. This indi- cates that the average deviation of the simulated values from the actual data is minimal, at only about 1.16, translating to an average deviation per floe of approximately 0.12. In the more complex scenario involving 30 floes, the model also performed with a PCC of 91.06% (Fig.5b), reinforcing its capability to maintain high accuracy even as the number of interactions increases. The RMSE for this setup was 3.01, set against a domain range of 200. Similar to the 10-floe scenario, this error metric suggests that the model maintains a high degree of precision. The average deviation across all 30 floes is about 3.01, with an individual floe deviating by only about 0.10 from the ground truth. Additionally, the visualizations from the simu- lations confirm that the model adheres to physical rules, showcasing trajectories where sea ice floes do not pass through each other or cross defined boundary limits. This aspect of the simulation is critical as it reflects the model’s ability to not only predict positions and movements accurately but also to ensure that these predictions respect the spatial constraints and interactions dictated by physical laws. Generalization ability To assess generalization beyond the training regime, we compared the simulation of the model against trajectories generated by DEM as ground truth. Generalization was quantified using PCC and RMSE. In climate-model evaluation studies, PCC values exceeding 0.8 are commonly inter- preted as indicating strong replication of observed spatial patterns [66, 67]. Under this criterion, the proposed model maintains high accuracy in generating a 20,000 time steps simulation, achiev- ing a PCC of 0.871 and a simulation RMSE of 2.334, which indicates that the learned simula- tion preserves the dominant spatiotemporal struc- ture of the DEM trajectories over a substantially extended horizon. As the simulation range increases shown in Table 2, PCC declines to 0.791 in generating 20,000 time steps simulation and further to 0.576 in generating 20,000 time steps simulation, while simulation RMSE increases to 7.230 and 11.055, respectively. These values indicate error accumula- tion and reduced pattern fidelity with the increase of simulation range. The generalization capability of the proposed newly developed model for simu- lating sea ice dynamics demonstrates remarkable performance, particularly in its ability to predict beyond the scope of its training data. Despite the training dataset comprising only 10,000-time steps, the model successfully projects plausible trajectories for up to 20,000 steps in simula- tions involving both 10 and 30 floes (Fig.F7 and Fig.F8). This capacity not only highlights the robustness of the model, but also its high fidelity in long-term prediction scenarios. The model’s ability to adhere to physical rules over these extended simulations is critical. It ensures that the trajectories of the sea ice floes do not unrealistically intersect or breach predefined boundaries, thus maintaining realistic simulations over time. By proving to be capable of reliable predictions over durations much longer than the training period, the model shows its potential to deliver insights and forecasts for unseen situations. 9 Table 1: Model performance comparison ModelSimulation RMSE one-step position RMSE one-step velocity RMSEPCC GNS [53]8.8974760.0010891.536386-0.243754 Interaction Network [35]11.4903810.0000380.3738040.567590 CN with one step4.4829290.0000100.0705390.712879 CN with Mish2.3340240.0000060.0284100.934317 CN with ReLU3.0484210.0000150.1171860.894427 CN with SiLU2.4529330.0000090.0652830.927664 Table 2: Model generalization comparison Simulation time range PCC Simulation RMSE 100000.9342.334 200000.8714.256 300000.7917.230 400000.57611.055 Efficient inference In the experiments, the performance of the pro- posed model is compared against a traditional model, which is DEM used for generating the ground truth, in simulations involving varying numbers of sea ice floes. The simulations were con- ducted first with 10 floes and then with 30 floes, recording the total execution time for each sce- nario. The results were highly favorable for the proposed model. When the simulation involved 10 floes, the running time for the proposed model was only 7.3 seconds compared to 7.7 seconds for the traditional model. Notably, as the complexity increased to 30 floes, the proposed model exhib- ited a minimal increase in running time, clocking in at 8.9 seconds. In contrast, the traditional model’s running time ballooned to 24 seconds under the same conditions. The proposed model demonstrated signifi- cant performance advantages over the traditional model in simulations with different numbers of floes. When tested with 10 floes, the improvement was modest, outperforming the traditional model by 5%. This slight edge can be attributed to the overhead involved in initializing and setting up the GPUs, although minimal, which affects per- formance when the computational scale is small [68]. However, as the complexity of the simu- lation increased to 30 floes, the benefits of the proposed model became more pronounced, with the proposed model achieving a 63% improve- ment in performance over the traditional model on the hardware configuration of the Intel i7- 11800H CPU and NVIDIA RTX 3080 Laptop GPU. The DEM baseline is executed on CPU NumPy (Numerical Python) [69], while the pro- posed method is accelerated by GPU and PyTorch [70, 71]. To assess performance on a more recent platform, the 30 floes case is also tested using the Intel Core i9-13900K CPU and an NVIDIA RTX 4090 GPU. Under this configuration, the proposed method requires 3.9 seconds on average, whereas the traditional DEM baseline requires 20 seconds. This difference indicates the GNN-based model’s ability to handle more complex interactions effi- ciently, making the proposed model particularly effective in scaled-up settings. However, the pro- posed model requires about 16 hours of training on a single NVIDIA RTX 4090 GPU to achieve high performance in simulation. Data assimilation To further enhance the accuracy of this model, the next step involves integrating data assimila- tion techniques, which will leverage observational data to refine the predictions and improve overall model performance. Incorporating data assim- ilation into the proposed model is crucial to address the inherent challenge of error accumu- lation during simulations. Even minor inaccura- cies in model states can magnify as simulations progress, leading to significant deviations from ground truth data. By integrating data assimila- tion, specifically through Ensemble Kalman filter (EnKF)[72, 73] or Ensemble transform Kalman filter (ETKF)[74, 75], observational data can sys- tematically be combined with model predictions. This method adjusts the model state by weight- ing the differences between observed values and 10 (a) PCC for 10 floes situation with the average of 98.98% (b) PCC for 30 floes situation with the average of 91.06% Fig. 5: The PCC plots for CN in inference. The x-axis is the index of floes, and the y-axis is the PCC value. (a) PCC for 10 floes situation with ENKF and ETKF (b) PCC for 30 floes situation with ENKF and ETKF Fig. 6: The PCC plots for CN with ENKF and ETKF in inference. The x-axis stands for different noise levels for model and observation, and the y-axis is the PCC value. Suffix stands for the observation frequency. For example, ENKF100 means the model couples the ENKF, and observes the half states with a specific noise level every 100-time steps. model predictions, effectively reducing the fore- cast error. This approach not only enhances the accuracy and reliability of the model’s predictions but also extends its applicative value in prac- tical scenarios, ensuring that it remains robust against the underlying uncertainties inherent in environmental modeling. The simulation framework developed in this study incorporates data assimilation techniques to address the inevitable accumulation of predic- tion error in long-term forecasting. The proposed model generates future states in a recurrent man- ner: only the first two time-step states are required as input for initialization, and for each subse- quent step, the model takes the two most recent 11 (a) RMSE for 10 floes situation with ENKF and ETKF (b) RMSE for 30 floes situation with ENKF and ETKF Fig. 7: The RMSE plots for CN with ENKF and ETKF in inference. The x-axis stands for different noise levels for model and observation, and the y-axis is the RMSE value. Suffix stands for the observation frequency. For example, ENKF 100 means the model couples the ENKF, and observes the half states with a specific noise level every 100-time steps. predicted states to forecast the next. While this design provides flexibility and computational effi- ciency, it also introduces the potential for error propagation over extended time horizons, as any modeling inaccuracies can compound recursively. To mitigate this error, the model integrates EnKF and ETKF methodologies into the simu- lation loop. These data assimilation techniques enable the system to incorporate partial observa- tions (positional data of sea ice floes) at specified intervals to recalibrate the internal model state. which reflects both observational availability and practical constraints in real sea-ice monitoring. In many operational and research settings, remote sensing provides snapshots of ice patterns, and ice motion is inferred by comparing the displace- ment of identifiable structures between successive images, rather than directly measuring complete dynamical states for every floe [76, 77]. Different model noises and observation noises are injected to represent uncertainty in both the learned simulator and the measurements to better reflect the reality. In the Kalman filtering formula- tion, the model is assumed imperfect with additive process noise, and observations are assumed noisy with additive observation noise; the analysis step then combines model forecasts and measurements using weights determined by their relative uncer- tainties. In this probabilistic framework, a larger specified model error increases the degree to which the analysis corrects toward observations, whereas a larger observation error yields analyses that remain closer to the model forecast. In the data assimilation approach for sea ice dynamics, the choice to use position data as the observational input for the ENKF and ETKF is primarily driven by a decision to align more closely with real-world conditions. The harsh and remote nature of polar environments imposes sig- nificant logistical challenges, which often restrict the types of data that can be consistently and reli- ably gathered. Position data are predominantly accessible through satellite imagery and aerial photography, which are effective for monitoring vast areas of sea ice [76, 77]. Thus, the obser- vation model is a linear model, and it observes only the position information and cannot observe any velocity information. For the total 10000- time steps, the observation occurs every 100-time steps, or 500-time steps. When there is no obser- vation, each ensemble member only applies the prediction model. Also, the well-trained proposed model is used as the prediction model for ENKF and ETKF. This ensures that the simulations are both scientifically sound and practically relevant, 12 adhering closely to the realities of polar research and observation capabilities. The results of data assimilation experiments using both the ENKF and the ETKF within the proposed model clearly demonstrate the critical impact of observation frequency on model perfor- mance. As illustrated in the results, a decrease in observation frequency from every 100 to 500-time steps significantly compromises the performance, evidenced by a dramatic increase in error and a sharp decline in the PCC for scenarios involving both 10 and 30 floes (Fig.6 and Fig.7). Fur- ther, even when some level of noise is introduced into the model and observation data, adjusting the observation frequency appropriately allows the model to manage these inaccuracies and still achieve high levels of accuracy (Fig.6 and Fig.7). This capability is particularly evident in more complex scenarios involving 30 floes, where the appropriate observation frequency setting enables the model to effectively maintain reliability and accuracy against noises. The underlying strength of the data assimilation process provides a robust framework that enhances the model’s ability to predict in long time trajectory against noises. 4 Conclusions This research introduces a novel application of proposed CN combined with DA techniques to improve the simulation of sea ice dynamics, mainly focusing on collision processes in one-dimensional settings. The findings indicate that the GNN- based model is capable of accurately simulating the dynamics of sea ice, adhering to physical laws, and handling data from various sources, thus confirming the robustness and versatility of the proposed approach. Moreover, the model significantly reduces computational demands com- pared to traditional numerical methods, facilitat- ing faster simulations without sacrificing accuracy. However, a notable limitation of the current imple- mentation is its restriction to one-dimensional simulations. This constraint simplifies the mod- eling but does not capture the full complexity of sea ice dynamics, which naturally occur in a two-dimensional space. The simplification to one dimension might overlook certain interactions and behaviors pertinent to the realistic depiction of sea ice processes. Extending the present one-dimensional frame- work to two dimensions introduces rotational degrees of freedom and frictional contact physics that fundamentally change what a learned sim- ulator must represent: in two dimensional set- tings, floes can rotate, and collisions generate not only normal impulses but also tangential con- tact forces that drive sliding–sticking behavior and produce torques, which cannot be captured if the model predicts translation alone. This struc- ture appears explicitly in established floe-scale DEM formulations, where translational momen- tum includes both contact-normal and contact- tangential forces, and angular momentum includes contact torques proportional to the lever arm crossed with the tangential contact force. Relative velocity and constraint forces have both normal and tangential components and are coupled to angular velocities and torques, making rotation and friction inseparable from collision modeling [38, 45, 49]. Consequently, a GNN-based model for two-dimensional floe collisions should model rotation and translation at the same time by aug- menting node states with orientation and angu- lar velocity and by learning edge-to-node torque transfer in addition to force transfer. Given the promising results achieved with the one-dimensional model, a logical extension of this research would be to adapt and test the GNN- based model within a two-dimensional framework. This expansion would allow the model to han- dle more complex interactions and provide a more comprehensive understanding of sea ice dynam- ics. A two-dimensional model would be able to simulate a broader range of physical phenomena and could potentially offer more precise predic- tions in more realistic situations. In summary, while the current model represents a fundamen- tal step toward applying GNNs to the simulation of sea ice dynamics, its extension to two dimen- sions represents a crucial next step in developing a comprehensive tool that can further effectively contribute to the field of climate research. Declarations Data availability All codes used to produce the raw data are available at the following repository: figshare[78]. 13 The details of implementation can be found in “Methods” and “Supplementary materials”. Code availability All codes used to produce the analysis are avail- able at the following repository: figshare[78]. The details of implementation can be found in “Meth- ods” and “Supplementary materials”. Author contributions statement R.Z. contributed to all aspects of the article, including the conception of the research idea, the development and design of the study, data analysis, and the drafting and revision of the manuscript. Correspondence Correspondence and requests for materials should be addressed to Ruibiao Zhu. Acknowledgements R. Zhu gratefully acknowledges the financial sup- port from the Australian Government Research Training Program (RTP) scholarship, as well as the computation support from project zv32 in National Computational Infrastructure (NCI). References [1] Massonnet, F., Vancoppenolle, M., Goosse, H., Docquier, D., Fichefet, T., Blanchard- Wrigglesworth, E.: Arctic sea-ice change tied to its mean state through thermody- namic processes. Nature Climate Change 8(7), 599–603 (2018) https://doi.org/10. 1038/s41558-018-0204-z [2] Letterly, A., Key, J., Liu, Y.: Arctic cli- mate: changes in sea ice extent outweigh changes in snow cover. The Cryosphere 12(10), 3373–3382 (2018) https://doi.org/10. 5194/tc-12-3373-2018 [3] Wunderling, N., Willeit, M., Donges, J.F., Winkelmann, R.: Global warming due to loss of large ice masses and arctic summer sea ice. Nature communications 11(1), 5177 (2020) [4] Notz, D., Stroeve, J.: Observed arctic sea- ice loss directly follows anthropogenic co2 emission. Science 354(6313), 747–750 (2016) https://doi.org/10.1126/science.aag2345 [5] Blockley, E., Vancoppenolle, M., Hunke, E., Bitz, C., Feltham, D., Lemieux, J.-F., Losch, M., Maisonnave, E., Notz, D., Ram- pal, P., Tietsche, S., Tremblay, B., Turner, A., Massonnet, F., ́ Olason, E., Roberts, A., Aksenov, Y., Fichefet, T., Garric, G., Iovino, D., Madec, G., Rousset, C., Melia, D.S., Schroeder, D.: The future of sea ice model- ing: Where do we go from here? Bulletin of the American Meteorological Society 101(8), 1304–1311 (2020) https://doi.org/10.1175/ BAMS-D-20-0073.1 [6] Mori, M., Kosaka, Y., Watanabe, M., Naka- mura, H., Kimoto, M.: A reconciled esti- mate of the influence of Arctic sea-ice loss on recent Eurasian cooling. Nature Climate Change 9(2), 123–129 (2019) https://doi. org/10.1038/s41558-018-0379-3 [7] Ogawa, F., Keenlyside, N., Gao, Y., et al.: Evaluating impacts of recent arctic sea ice loss on the northern hemisphere winter cli- mate change. Geophysical Research Letters 45(7), 3255–3263 (2018) https://doi.org/10. 1002/2017GL076502 [8] Stroeve, J.C., Markus, T., Boisvert, L., Miller, J., Barrett, A.: Changes in arctic melt season and implications for sea ice loss. Geophysical Research Letters 41(4), 1216–1225 (2014) https://doi.org/10.1002/ 2013GL058951 [9] Cunningham, C.X., Williamson, G.J., Bow- man, D.M.: Increasing frequency and inten- sity of the most extreme wildfires on earth. Nature Ecology & Evolution, 1–6 (2024) [10] Cr ́epin, A.S., Karcher, M., Gascard, J.C.: Arctic climate change, economy and soci- ety (access): Integrated perspectives. Ambio 46(Suppl 3), 341–354 (2017) https://doi.org/ 10.1007/s13280-017-0953-3 [11] Alvarez, J., Yumashev, D., Whiteman, 14 G.: A framework for assessing the eco- nomic impacts of arctic change. Ambio 49(2), 407–418 (2020) https://doi.org/10. 1007/s13280-019-01211-z [12] Hou, Z., Liu, J., Pang, J.: Green transition and climate adaptation: empirical evidence of corporate resilience in the low-carbon economy. GLOBAL NEST JOURNAL 27(6) (2025) https://doi.org/10.30955/gnj.07421 [13] Ani, C.J., Robson, B.: Responses of marine ecosystems to climate change impacts and their treatment in biogeochemical ecosys- tem models. Marine Pollution Bulletin 166, 112223 (2021) https://doi.org/10.1016/ j.marpolbul.2021.112223 [14] Lubitz, N., Daly, R., Smoothey, A.F., Vianello, P., Roberts, M.J., Schoeman, D.S., Sheaves, M., Cowley, P.D., Dagorn, L., For- get, F.G., et al.: Climate change-driven cool- ing can kill marine megafauna at their distri- butional limits. Nature Climate Change, 1–10 (2024) [15] Tollenaar, V., Zekollari, H., Kittel, C., Farinotti, D., Lhermitte, S., Debaille, V., Goderis, S., Claeys, P., Joy, K.H., Pattyn, F.: Antarctic meteorites threatened by cli- mate warming. Nature Climate Change, 1–4 (2024) [16] Rockstr ̈om, J., Gupta, J., Qin, D., et al.: Safe and just earth system boundaries. Nature 619(7968), 102–111 (2023) https://doi.org/ 10.1038/s41586-023-06083-8 [17] Xu, Z., Tartakovsky, A.M., Pan, W.: Discrete-element model for the interaction between ocean waves and sea ice. Phys. Rev. E 85, 016703 (2012) https://doi.org/10. 1103/PhysRevE.85.016703 [18] Herman,A.:Discrete-elementbonded- particle sea ice model design, version 1.3 a–model description and implementation. GeoscientificModelDevelopment 9(3), 1219–1241 (2016) [19] West, B., O’Connor, D., Parno, M., Krackow, M., Polashenski, C.: Bonded discrete ele- ment simulations of sea ice with non-local failure: Applications to nares strait. Jour- nal of Advances in Modeling Earth Sys- tems 14(6) (2022) https://doi.org/10.1029/ 2021MS002614 [20] Deng, Q., Stechmann, S.N., Cheng, N.: Particle-continuum multiscale modeling of sea ice floes. SIAM Multiscale Modeling and Simulation (MMS) (2023) [21] Montemuro, B.P., Manucharyan, G.E.: Sub- zero: a discrete element sea ice model that simulates floes as evolving concave polygons. Journal of Open Source Software 8(88), 5039 (2023) https://doi.org/10.21105/joss.05039 [22] Hopkins, M.A.: A discrete element lagrangian sea ice model. Engineering Computations 21(2/3/4), 409–421 (2004) https://doi.org/ 10.1108/02644400410519857 [23] Turner, A.K., Peterson, K.J., Bolintineanu, D.: Geometric remapping of particle distri- butions in the discrete element model for sea ice (demsi v0.0). Geoscientific Model Devel- opment 15(5), 1953–1970 (2022) https://doi. org/10.5194/gmd-15-1953-2022 [24] Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM, Philadelphia (2005). https://doi.org/ 10.1137/1.9780898717921 [25] Isakov, V.: Inverse Problems for Partial Dif- ferential Equations. Springer, Philadelphia (2018) [26] Ying, R., He, R., Chen, K., Eksombatchai, P., Hamilton, W.L., Leskovec, J.: Graph convo- lutional neural networks for web-scale recom- mender systems. In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. KDD ’18, p. 974–983. Association for Comput- ing Machinery, New York, NY, USA (2018). https://doi.org/10.1145/3219819.3219890 [27] Gligorijevi ́c, V., Renfrew, P.D., Kosciolek, T., Leman, J.K., Berenberg, D., Vatanen, T., Chandler, C., Taylor, B.C., Fisk, I.M., 15 Vlamakis, H., et al.: Structure-based pro- tein function prediction using graph con- volutional networks. Nature communications 12(1), 3168 (2021) [28] Merchant, A., Batzner, S., Schoenholz, S.S., Aykol, M., Cheon, G., Cubuk, E.D.: Scaling deep learning for materials discovery. Nature 624(7990), 80–85 (2023) [29] Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., Philip, S.Y.: A comprehensive survey on graph neural networks. IEEE transactions on neural networks and learning systems 32(1), 4–24 (2020) [30] Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., Sun, M.: Graph neural networks: A review of methods and applications. AI open 1, 57–81 (2020) [31] Liu, J., Hou, Z., Liu, B., Zhou, X.: Math- ematical and machine learning innovations for power systems: Predicting transformer oil temperature with beluga whale optimization- based hybrid neural networks. Mathematics 13(11), 1785 (2025) https://doi.org/10.3390/ math13111785 [32] Hou, Z., Liu, J., Shao, Z., Ma, Q., Liu, W.: Machine learning innovations in renewable energy systems with integrated nrbo-txad for enhanced wind speed forecasting accuracy. Electronics 14(12), 2329 (2025) https://doi. org/10.3390/electronics14122329 [33] Xie, J., He, J., Gao, Z., Wang, S., Liu, J., Fan, H.: An enhanced snow abla- tion optimizer for uav swarm path plan- ning and engineering design problems. Heliyon 10(18) (2024) https://doi.org/10. 1016/j.heliyon.2024.e37819 [34] Liu, J., Duan, Z., Hu, X., Zhong, J., Yin, Y.: Detracking autoencoding conditional genera- tive adversarial network: Improved generative adversarial network method for tabular miss- ing value imputation. Entropy 26(5), 402 (2024) https://doi.org/10.3390/e26050402 [35] Battaglia,P.,Pascanu,R.,Lai,M., Jimenez Rezende, D., et al.: Interaction net- works for learning about objects, relations and physics. Advances in neural information processing systems 29 (2016) [36] Hopkins, M.A., Tuhkuri, J.: Compres- sion of floating ice fields. Journal of Geophysical Research: Oceans 104(C7), 15815–15825 (1999) https://doi.org/10.1029/ 1999JC900127 [37] Sun, S., Shen, H.H.: Simulation of pancake ice load on a circular cylinder in a wave and current field. Cold Regions Science and Tech- nology 78, 31–39 (2012) https://doi.org/10. 1016/j.coldregions.2012.02.003 [38] Damsgaard, A., Adcroft, A., Sergienko, O.: Application of discrete element methods to approximate sea ice dynamics. Journal of Advances in Modeling Earth Systems 10(9), 2228–2244 (2018) https://doi.org/10.1029/ 2018MS001299 [39] Obermayr, M., Dressler, K., Vrettos, C., Eberhard, P.: A bonded-particle model for cemented sand. Computers and Geotechnics 49, 299–313 (2013) https://doi.org/10.1016/ j.compgeo.2012.09.001 [40] Ergenzinger, C., Seifried, R., Eberhard, P.: A discrete element model to describe failure of strong rock in uniaxial compression. Granu- lar Matter 13(4), 341–364 (2011) https://doi. org/10.1007/s10035-010-0230-7 [41] Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. G ́eotechnique 29(1), 47–65 (1979) https:// doi.org/10.1680/geot.1979.29.1.47 [42] Damsgaard, A., Egholm, D.L., Beem, L.H., Tulaczyk, S., Larsen, N.K., Piotrowski, J.A., Siegfried, M.R.: Ice flow dynamics forced by water pressure variations in subglacial granular beds. Geophysical Research Letters 43(23), 12165–12173 (2016) https://doi.org/ 10.1002/2016GL071579 [43] Ergenzinger, C., Seifried, R., Eberhard, P.: A discrete element model to describe failure of strong rock in uniaxial compression. Granular 16 Matter 13(4), 341–364 (2011) [44] Luding, S.: Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory. European journal of environmental and civil engineering 12(7-8), 785–826 (2008) [45] Damsgaard, A., Sergienko, O., Adcroft, A.: The effects of ice floe-floe interactions on pressure ridging in sea ice. Journal of Advances in Modeling Earth Systems 13(7), 2020–002336 (2021) [46] Lund, B., Graber, H.C., Persson, P., Smith, M., Doble, M., Thomson, J., Wadhams, P.: Arctic sea ice drift measured by ship- board marine radar. Journal of Geophysical Research: Oceans 123(6), 4298–4321 (2018) [47] J. Hanley, K., O’Sullivan, C.: Analytical study of the accuracy of discrete element simulations. International Journal for Numer- ical Methods in Engineering 109(1), 29–51 (2017) [48] O’Sullivan, C., Bray, J.: Selecting a suitable time step for discrete element simulations that use the central difference time inte- gration scheme. Engineering Computations 21, 278–303 (2004) https://doi.org/10.1108/ 02644400410519794 [49] Tuhkuri, J., Poloj ̈arvi, A.: A review of dis- crete element simulation of ice–structure interaction. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376(2129), 20170335 (2018) [50] Schulson, E.M.: The structure and mechani- cal behavior of ice. JOM 51(2), 21–27 (1999) https://doi.org/10.1007/s11837-999-0206-4 [51] Petrovic, J.J.: Review mechanical properties of ice and snow. Journal of Materials Science 38(1), 1–6 (2003) https://doi.org/10.1023/A: 1021134128038 [52] Greydanus, S., Dzamba, M., Yosinski, J.: Hamiltonian neural networks. Advances in neural information processing systems 32 (2019) [53] Sanchez-Gonzalez, A., Godwin, J., Pfaff, T., Ying, R., Leskovec, J., Battaglia, P.: Learn- ing to simulate complex physics with graph networks. In: International Conference on Machine Learning, p. 8459–8468 (2020). PMLR [54] Elman, J.L.: Finding structure in time. Cog- nitive Science 14(2), 179–211 (1990) https: //doi.org/10.1207/s15516709cog1402 1 [55] Merkwirth, C., Lengauer, T.: Automatic generation of complementary descriptors with molecular graph networks. Journal of Chemical Information and Modeling 45(5), 1159–1168 (2005) https://doi.org/10.1021/ ci049613b . PMID: 16180893 [56] Scarselli, F., Gori, M., Tsoi, A.C., Hagen- buchner, M., Monfardini, G.: The graph neu- ral network model. IEEE Transactions on Neural Networks 20, 61–80 (2009) [57] Misra, D.: Mish: A self regularized non- monotonicactivationfunction.BMVC (2020) [58] Hahnloser, R.H., Sarpeshkar, R., Mahowald, M.A., Douglas, R.J., Seung, H.S.: Digital selection and analogue amplification coex- ist in a cortex-inspired silicon circuit. nature 405(6789), 947–951 (2000) [59] Kingma, D.P.: Adam: A method for stochas- tic optimization. ICLR (2015) [60] Pham, H., Guan, M., Zoph, B., Le, Q., Dean, J.: Efficient neural architecture search via parameters sharing. In: International Con- ference on Machine Learning, p. 4095–4104 (2018). PMLR [61] Hendrycks, D.: Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415 (2016) [62] Ramachandran, P., Zoph, B., Le, Q.V.: Searching for activation functions. arXiv preprint arXiv:1710.05941 (2017) 17 [63] Nikolentzos, G., Dasoulas, G., Vazirgiannis, M.: K-hop graph neural networks. Neural Networks 130, 195–205 (2020) [64] Wu, X., Ajorlou, A., Wu, Z., Jadbabaie, A.: Demystifying oversmoothing in attention- based graph neural networks. Advances in Neural Information Processing Systems 36, 35084–35106 (2023) [65] Kelesis, D., Fotakis, D., Paliouras, G.: Ana- lyzing the effect of residual connections to oversmoothing in graph neural networks. Machine Learning 114(8), 184 (2025) [66] Salazar, ́ A., Thatcher, M., Goubanova, K., Bernal, P., Guti ́errez, J., Squeo, F.: Cmip6 precipitation and temperature projections for chile. Climate Dynamics 62(3), 2475–2498 (2024) [67] Christophersen, J.A., Rydbeck, A., Flatau, M., Janiga, M., Reynolds, C.A., Jensen, T., Smith, T.: Oceanic rossby wave predictabil- ity in ecmwf’s subseasonal-to-seasonal refore- casts. Quarterly Journal of the Royal Meteo- rological Society 150(759), 1048–1067 (2024) [68] Dubey, S.P., Kumar, M.S., Balaji, S.: Gpu computing for compute-intensive scientific calculation. In: Das, K.N., Bansal, J.C., Deep, K., Nagar, A.K., Pathipooranam, P., Naidu, R.C. (eds.) Soft Computing for Prob- lem Solving, p. 131–140. Springer, Singapore (2020) [69] Harris, C.R., Millman, K.J., Walt, S.J., Gom- mers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., Kerkwijk, M.H., Brett, M., Haldane, A., R ́ıo, J.F., Wiebe, M., Peterson, P., G ́erard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (2020) https://doi.org/ 10.1038/s41586-020-2649-2 [70] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.: Pytorch: An imperative style, high-performance deep learning library. Advances in neural informa- tion processing systems 32 (2019) [71] Paszke,A.,Gross,S.,Chintala,S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch. In: NIPS 2017 Workshop on Autodiff (2017). https://openreview.net/forum?id=BJJsrmfCZ [72] Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans 99(C5), 10143–10162 (1994) https: //doi.org/10.1029/94JC00572 [73] Burgers,G.,Leeuwen,P.J.,Evensen, G.:Analysisschemeintheensem- blekalmanfilter.MonthlyWeather Review126(6),1719–1724(1998) https://doi.org/10.1175/1520-0493(1998) 126⟨1719:ASITEK⟩2.0.CO;2 [74] Bishop, C.H., Etherton, B.J., Majum- dar,S.J.:Adaptivesamplingwith theensembletransformkalmanfilter. parti:Theoreticalaspects.Monthly Weather Review 129(3), 420–436 (2001) https://doi.org/10.1175/1520-0493(2001) 129⟨0420:ASWTET⟩2.0.CO;2 [75] Harlim, J., Hunt, B.R.: A non-gaussian ensemble filter for assimilating infrequent noisy observations. Tellus A: Dynamic Mete- orology and Oceanography 59(2), 225–237 (2007) [76] Lavergne, T., Pi ̃nol Sol ́e, M., Down, E., Donlon, C.: Towards a swath-to-swath sea- ice drift product for the copernicus imag- ing microwave radiometer mission. The Cryosphere 15(8), 3681–3698 (2021) [77] Wang, X., Chen, R., Li, C., Chen, Z., Hui, F., Cheng, X.: An intercomparison of satel- lite derived arctic sea ice motion products. Remote Sensing 14(5), 1261 (2022) [78] Zhu, R.: Graph neural network for colliding particles with an application to sea ice floe modeling. figshare https://doi.org/10.6084/ 18 m9.figshare.28658180 (2025) [79] Liu, J., Hou, Z., Xu, Y., He, Y., Wang, B.: Enhanced beluga whale optimization meets gru and adaptive cubature kalman filter: a novel approach for state of charge estima- tion in lithium-ion batteries. Ionics 31(10), 10643–10669 (2025) https://doi.org/10.1007/ s11581-025-06578-6 [80] Barab ́asi, A.-L.: Network science. Philosoph- ical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sci- ences 371(1987), 20120375 (2013) [81] Cai, M., Gu, H., Hong, P., Li, J.: A combina- tion of physics-informed neural networks with the fixed-stress splitting iteration for solving biot’s model. Frontiers in Applied Mathemat- ics and Statistics 9, 1206500 (2023) [82] Zhang, Y., Ruan, W., Wang, F., Huang, X.: Generalizing universal adversarial pertur- bations for deep neural networks. Machine Learning 112(5), 1597–1626 (2023) [83] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial exam- ples. CoRR abs/1412.6572 (2014) [84] Scher, S., Messori, G.: Generalization proper- ties of feed-forward neural networks trained on lorenz systems. Nonlinear Processes in Geophysics 26(4), 381–399 (2019) https:// doi.org/10.5194/npg-26-381-2019 [85] Lu, L., Meng, X., Mao, Z., Karniadakis, G.E.: Deepxde: A deep learning library for solving differential equations. SIAM Review 63(1), 208–228 (2021) https://doi.org/10. 1137/19M1274067 [86] D’Angelo, F., Andriushchenko, M., Varre, A.V., Flammarion, N.: Why do we need weight decay in modern deep learning? Advances in Neural Information Processing Systems 37, 23191–23223 (2024) [87] Senior, A., Heigold, G., Ranzato, M., Yang, K.: An empirical study of learning rates in deep neural networks for speech recognition. In: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, p. 6724–6728 (2013). IEEE [88] Ren, Y., Ma, C., Ying, L.: Understanding the generalization benefits of late learning rate decay. In: International Conference on Artifi- cial Intelligence and Statistics, p. 4465–4473 (2024). PMLR [89] Li, Y., Wei, C., Ma, T.: Towards explain- ing the regularization effect of initial large learning rate in training neural networks. Advances in neural information processing systems 32 (2019) [90] Hawkins, D.M.: The problem of overfitting. Journal of chemical information and com- puter sciences 44(1), 1–12 (2004) [91] Ying, X.: An overview of overfitting and its solutions. In: Journal of Physics: Conference Series, vol. 1168 (2019). IOP Publishing Appendix A Supplementary information: kalman filter (KF) The observation operator H is intentionally par- tial in both Ensemble kalman filter (EnKF) and Ensemble transform Kalman filter (ETKF). First, only half of the state variables were observed at each assimilation time. This choice reflects a standard geophysical data-assimilation setting where the observation dimension is substantially smaller than the full system state, i.e., only a subset of variables and locations are measured, and the filter must infer the remainder through the forecast–analysis coupling. Second, we assimi- lated positions only (floe locations), while treating velocities as unobserved. This design is moti- vated by observational availability and practical constraints in sea-ice monitoring: satellite remote sensing provides widespread spatial coverage but does so primarily through sequential imagery, where motion (and hence velocity) is typically derived by tracking ice features between images rather than directly measured as a state variable [76]. Consequently, using positions as the primary 19 assimilated quantity better matches the informa- tion content and reliability of common observa- tional pipelines, while allowing the filter to infer velocity implicitly through the dynamical model and ensemble cross-covariances. This position- only, partially observed setup therefore empha- sizes practical realism and provides a more faithful test of whether the proposed model remains effec- tive when the observation stream is sparse, noisy, and incomplete—conditions that are typical in operational polar data assimilation [77]. All data assimilation experiments are con- ducted using an ensemble size of I = 100 for both EnKF and ETKF, with identical settings to ensure a controlled comparison. Covariance infla- tion is 1 for ETKF. In the experiments, model and observation error covariance matrices Q and R are constructed as diagonal (square) covari- ance matrices: all off-diagonal entries are set to zero, and diagonal entries are set to σ 2 with σ ∈ 0.1, 1, 2, corresponding to low-, moderate-, and high-noise regimes. For each ensemble member i ∈ I, the model perturbation q i and observation perturbation r i are independently sampled from multivariate normal distributions q i ∼ N (0,Q) and r i ∼N (0,R) with 0 mean. A.1 Ensemble kalman filter (EnKF) The Ensemble Kalman Filter is a Monte Carlo method to handle the state estimation problem in high-dimensional systems where the traditional KF becomes computationally impractical [72, 73]. It uses a collection (ensemble) of system states to represent the probability distribution of the state estimates. Each member of the ensemble is updated based on the measurement and the mean and covariance of the ensemble. Compared to Extended Kalman Filter (EKF) [79], EnKF does not need any tangent approximation of the nonlinear prediction function. This approach allows for better handling of the uncertainties and nonlinearity in the system by using multiple par- allel instances of the filter [72, 73]. The Ensemble Kalman filter (EnKF) works as follows: Algorithm 1 Ensemble Kalman filter Require: Known prediction model g and known observation model H Q,R,T,I,n← Initialized t j = j T n ,j ∈ [1,...,n] while i≤ I do u i ← Initialized end while while t j ≤ T do if There is an observation then i← 1 while i≤ I do q i ∼ Q u i t j ←u i t j +q i end while ̄u t j ← 1 I P I i=1 u i t j U ← u 1 t j − ̄u t j ,u 2 t j − ̄u t j , ...u I t j − ̄u t j if The observation model H is nonlinear then V ← H(u 1 t j )− H( ̄u t j ), ... H(u I t j )− H( ̄u t j ) else if The observation model H is linear then V ← HU end if K ← 1 I−1 UV T ( 1 I−1 V V T + R) −1 y t j ← observation data i← 1 while i≤ I do r i ∼ R y i t j ←y t j +r i u i t j ←u i t j + K(y i t j − H(u i t j )) end while else if There is no observation then i← 1 while i≤ I do q i ∼ Q u i t j ← g(u i t j ) +q i end while end if end while ̄u← 1 I P I i=1 u i A.1.1 Prediction step Letu t j be one of the ensemble members at time t j . If there are I ensemble members,u i t j =u i t j +q i whereq i samples from Q, and the ensemble mean ̄u t j would be ̄u t j = 1 I P I i=1 u i t j . The matrix U ∈ R N×I and V ∈ R M×I are formed, where N is the number of states in the 20 prediction model, and M is the number of states in the observation model, U = u 1 t j − ̄u t j ,u 2 t j − ̄u t j , ...u I t j − ̄u t j , and V = H(u 1 t j )− H( ̄u t j ), H(u 2 t j )− H( ̄u t j ), ... H(u I t j )− H( ̄u t j ) , where H is the nonlinear observation operator. When observation operator is linear called H, V = HU A.1.2 Update step Compute the Kalman Gain K = 1 I − 1 UV T ( 1 I − 1 V V T + R) −1 (A1) Perturb the observation to approximate the posterior error covariance matrix [73], y i t j =y t j +r i (A2) wherey t j is the observation at time t j andr i samples from R. Update the estimate via measure- ment u i t j =u i t j + K(y i t j − H(u i t j ))(A3) A.2 Ensemble transform kalman filter (ETKF) The Ensemble Transform Kalman Filter is a vari- ant of the ensemble Kalman filter [74, 75]. It focuses on transforming the ensemble members to better align with the measurement update. This transformation ensures that the ensemble main- tains the correct statistical properties (mean and covariance) after the update step. The ETKF uses a deterministic approach to generate the ensemble members, making it more efficient in specific appli- cations, particularly those involving geophysical and meteorological forecasting [74, 75]. TheEnsembletransformKalmanfilter (ETKF) works as follows: A.2.1 Prediction step Letu t j be one of the ensemble members at time t j . If there are I ensemble members,u i t j =u i t j +q i Algorithm 2 Ensemble transform Kalman filter Require: Known prediction model g and known observation model H Q,R,T,I,n← Initialised t j ← j T n ,j ∈ [1,...,n] while i≤ I do u i ← Initialized end while while t j ≤ T do if There is an observation then i← 1 while i≤ I do q i ∼ Q u i t j ←u i t j +q i end while ̄u t j ← 1 I P I i=1 u i t j U ← u 1 t j − ̄u t j ,u 2 t j − ̄u t j , ...u I t j − ̄u t j if The observation model H is nonlinear then V ← H(u 1 t j )− H( ̄u t j ), ... H(u I t j )− H( ̄u t j ) else if The observation model H is linear then V ← HU end if J ← I−1 1+r I + V T (R) −1 V XΣX T ← svd(J )▷ Singular value decomposition K ← U (J ) −1 V T (R) −1 T ← √ I − 1XΣ − 1 2 X T ▷ The transformation matrix y t j ← observation data ̄u t j ← ̄u t j + K(y t j − H( ̄u t j )) U ← UT i← 1 while i≤ I do u i t j ←u i t j + ̄u t j end while else if There is no observation then i← 1 while i≤ I do q i ∼ Q u i t j ← g(u i t j ) +q i end while end if end while ̄u← 1 I P I i=1 u i 21 whereq i samples from Q, and the ensemble mean ̄u t j would be ̄u t j = 1 I P I i=1 u i t j . The matrix U ∈ R N×I and V ∈ R M×I are formed, where N is the number of states in the prediction model, and M is the number of states in the observation model, U = u 1 t j − ̄u t j ,u 2 t j − ̄u t j , ...u I t j − ̄u t j , and V = H(u 1 t j )− H( ̄u t j ), H(u 2 t j )− H( ̄u t j ), ... H(u I t j )− H( ̄u t j ) , where H is the nonlinear observation operator. When observation operator is linear called H, V = HU A.2.2 Update step Compute the singular value decomposition of the I × I matrix J = I − 1 1 + r I + V T (R) −1 V = XΣX T (A4) where r is the inflation number. Compute the Kalman Gain K = U (XΣ −1 X T )V T (R) −1 (A5) Update the ensemble mean ̄u t j = ̄u t j + K(y t j − H( ̄u t j ))(A6) Compute the transformation matrix T = √ I − 1XΣ − 1 2 X T (A7) Update posterior ensemble U = UT(A8) u i t j =u i t j + ̄u t j (A9) A.3 Results for CN with DA in 10 floes and 30 floes situation The following tables are the PCC and RMSE for ENKF and ETKF in 10 floes and 30 floes situ- ations with different noise levels. Suffix 100 and 500 means the proposed model observes the half states every 100 and 500 time steps, respectively. Appendix B Supplementary information: sender and receiver relation matrices B.1 Definition of sender and receiver relation matrices The sender relation matrix R s and the receiver relation matrix R r are both binary matrices that represent the connections between nodes in a graph through directed edges. Here is a detailed definition of each, including the indexing of nodes and edges: B.1.1 Indexing of nodes and edges Before defining the matrices, it is crucial to estab- lish the indexing for nodes and edges: • Nodes in the graph are indexed arbitrarily from 1 to n, where n is the total number of nodes in the graph. • Edges are also indexed arbitrarily from 1 to e, where e is the total number of directed edges in the graph. • The ordering of these indices does not impact the structural or functional properties of the matrices; it only affects the representation. B.1.2 Sender relation matrix • Dimensions: R s has dimensions n× e, where n is the number of nodes and e is the number of edges in the graph. • Entries: The entries of R s are binary: – R s [i,j] = 1 if the j-th edge originates from node i. – R s [i,j] = 0 if the j-th edge does not originate from node i. This matrix is used to identify which node is the sender or the starting point for each directed edge in the graph. B.1.3 Receiver relation matrix • Dimensions: R r also has dimensions n× e. • Entries: The entries of R r are binary: – R r [i,j] = 1 if the j-th edge terminates at node i. 22 Table A1: PCC for ENKF and ETKF in 10 floes with different noise levels PCC Model noiseObservation noiseENKF 100ETKF 100ENKF 500ETKF 500 0.10.10.993480.993130.901780.89945 0.11.00.992200.992070.897650.89774 0.12.00.990420.990730.899080.89312 1.00.10.990400.990570.888160.88732 1.01.00.990050.990170.890260.89120 1.02.00.990420.990480.874180.88799 2.00.10.986720.986210.868180.87926 2.01.00.987220.987000.886940.88400 2.02.00.986960.986950.873850.88544 Table A2: RMSE for ENKF and ETKF in 10 floes with different noise levels RMSE Model noiseObservation noiseENKF 100ETKF 100ENKF 500ETKF 500 0.10.10.871840.890753.479523.51363 0.11.00.988090.996653.567903.56005 0.12.01.108621.096383.575053.64760 1.00.11.081391.070643.733823.83360 1.01.01.094881.083253.704203.68827 1.02.01.076191.081633.978733.85944 2.00.11.295091.289584.063473.97697 2.01.01.269321.260453.979173.98605 2.02.01.278061.261984.162283.91707 – R r [i,j] = 0 if the j-th edge does not termi- nate at node i. This matrix indicates which node is the receiver or the endpoint for each directed edge in the graph. B.2 Definition and proof of the column permutation operator in bi-directional graph Given the sender relation matrix R s and the receiver relation matrix R r for a graph, where each directed edge in the graph is represented by an entry in these matrices, a column permutation operator P is defined that reorders the columns of one matrix to transform it into the other, under the condition that each relation between two nodes is represented by two directional edges. B.2.1 Setup • Graph Definition: Assume a graph with nodes and directed edges such that each rela- tionship between two nodes u and v is repre- sented by two directed edges: one from u to v and the other from v to u. • Matrix Dimensions: R s and R r are n × e matrices, where n is the number of nodes and e is the number of directed edges, e is not neces- sarily equals to (n− 1)× 2 since there are some nodes might not be connected [80]. B.2.2 Column permutation operator • Definition: Column permutation operator P is an e× e permutation matrix, where each row and each column contains exactly one entry of 1, and all other entries are 0. 23 Table A3: PCC for ENKF and ETKF in 30 floes with different noise levels PCC Model noiseObservation noiseENKF 100ETKF 100ENKF 500ETKF 500 0.10.10.993480.991350.897350.89360 0.11.00.989320.989510.897130.89262 0.12.00.987490.987130.891870.89227 1.00.10.987550.987780.877430.88166 1.01.00.987000.987390.873680.88840 1.02.00.986860.987080.890460.87620 2.00.10.983610.983470.868260.86095 2.01.00.982420.982840.870360.86763 2.02.00.982100.982390.862200.86103 Table A4: RMSE for ENKF and ETKF in 30 floes with different noise levels RMSE Model noiseObservation noiseENKF 100ETKF 100ENKF 500ETKF 500 0.10.10.847000.849203.075653.11070 0.11.00.978520.970213.079693.14298 0.12.01.088281.100313.167823.15995 1.00.11.040121.035043.375143.36270 1.01.01.074271.065593.445873.25165 1.02.01.086451.064423.224793.44050 2.00.11.218801.222603.602533.69672 2.01.01.265261.247983.543383.62936 2.02.01.269871.253113.728763.71761 • Action: P acts on the columns of R s such that R s P reorders the columns of R s . The goal is to find a permutation P such that R s P = R r . B.2.3 Proving the existence of permutation operator To demonstrate that such a permutation matrix P exists, the prove as follows: 1. Bijective mapping of edges: Each edge i in R s representing a directed edge from node u to node v has a corresponding edge j in R r repre- senting the directed edge from node v to node u, given the assumption that each relationship is bidirectional. 2. Construction of P : • For each column i in R s (representing edge i : u→ v), there exists a unique column j in R r (representing edge j : v → u). • Set P ji = 1 if and only if the edge represented by column i in R s corresponds to the edge represented by column j in R r . 3. Validity of P : • Since each directed edge in the graph is unique and bidirectional, the mapping from edges in R s to R r is bijective (one-to-one and onto). Thus, P constructed in this manner is a valid permutation matrix. 4. Permutation Result: Applying P to R s , R s P is constructed. By the construction of P , R s P will have its columns rearranged to match exactly the columns of R r . 24 B.3 Properties of the permutation matrix • Each row and each column of P contains exactly one entry of 1. • The matrix P rearranges (permutes) the columns of any matrix it multiplies. • P is symmetric, meaning that P = P T . This property arises because ,for every entry P ij = 1, where i is permuted to j, there is a correspond- ing entry P ji = 1 reflecting the inverse of that permutation, which, for a permutation matrix that swaps pairs of indices, means i and j swap roles symmetrically. B.4 Multiplication of permutation matrices The product P×P results in another matrix where each entry (i,j) is computed by: (P × P ) ij = n X k=1 P ik P kj This formula considers the intermediate step k, where P ik = 1 if the permutation moves column i to column k, and P kj = 1 if it subsequently moves column k to column j. B.5 Identity matrix The identity matrix I is defined as: I ij = ( 1 if i = j 0 if i̸= j It leaves any vector or matrix unchanged when multiplied. B.6 Proof of relationship between identity matrix and permutation matrix By definition, the property has P n k=1 P ik P kj = 0 if and only if i = j. Therefore, 1. Diagonal elements: For each i, (P×P ) i = 1 because the column i moves to some column k and then back to i due to the permutation, ensuring it returns to its original position. 2. Off-diagonal elements: For each i ̸= j, (P × P ) ij = 0 because there is no intermedi- ate column k such that column i would move to column j through k. The only contributing k would have to coincide with both i and j, which cannot happen if i̸= j. Thus, P × P = I. B.6.1 Theoretical result Thus, the column permutation matrix P exists such that R s P = R r , under the assumption of bidirectional relationships represented as pairs of directed edges. This shows that P correctly trans- lates the sender matrix into the receiver matrix by simply reordering the edges according to their reverse directions. Lemma 1 Given relation matrices R s and R r for a graph where a relation between two nodes is repre- sented as two directional edges, there exists a column permutation of one of the matrices to represent the other one. The permutation would be column i and col- umn j when two directional edges are indexed as i and j between two nodes. Utilizing column permutation to represent R s is particularly advantageous in terms of mem- ory storage. Instead of storing a completely new matrix R s , the matrix R r can be simply stored, and R s matrix can be derived by column permu- tation. For example, if the graph is constructed as Figure B1b, the relation matrices R s and R r would be as follows: R s = edge 1 edge 2 edge 3 edge 4 edge 5 edge 6 node 1 101000= 2 node 2 010001= 2 node 3 000110= 2 = 1= 1= 1= 1= 1= 1 R r = edge 1 edge 2 edge 3 edge 4 edge 5 edge 6 node 1 010100= 2 node 2 100010= 2 node 3 001001= 2 = 1= 1= 1= 1= 1= 1 It can be identified that the columns 2, 4 and 6 permutations of R r into columns 1, 3 and 5 would 25 (a) A graph(b) A bi-directional graph Fig. B1: Illustration of constructing graphs (Circle represents node, and arrow represents edge) be identical with R s . Additionally, the column permutation matrix can be obtained by definition as follows: P = 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 The result can be verified by checking R s × P = R r , as follows: 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 × 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 = 0 1 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 1 Remark 1 Since P is an involutory matrix, both R s × P = R r and R r × P = R s holds. Also, P ij = 1 if and only if labeling a pair of edges as i and j. Thus, applying the permutation requires far less memory than storing an entirely new matrix. For the 10 floes setup as Fig.4, the same defi- nitions are used from B.1.2 and B.1.3, and the R s and R r can be constructed as the following: For more realistic two-dimensional case, both R s and R r can be constructed with the same matrix definition from B.1.2 and B.1.3. If there are n floes in the two dimensional space, then each floe might interact with at most n− 1 other floes. Thus, n(n−1) edges are formed and indexed. Both R s and R r have n rows and n(n − 1) columns, and their values are either 1 or 0, where 1 rep- resents the edge that exists, and 0 represents the edge that does not exist. The sum of each row of the matrix is n− 1, and the sum of each column of the matrix is 1. Appendix C Supplementary information: Error analysis Quantifying the errors in neural networks is an intricate challenge that continues to be an active area of research [81]. For example, adding a small noise to perturb the input can result in different classifications in the image [82, 83]. The accuracy of neural network approximations is influenced by multiple factors, including network architecture, the characteristics of the training data, and the specific implementations of learning algorithms. These networks operate within complex, often non-convex optimization landscapes, where the choices of activation functions, network topology, and the training regime all contribute significantly to the network’s performance and generalization capabilities outside the training set [84, 85]. Based on the mean value theorem, if f (x) is continuous on x,y and differentiable on (x,y), then ∃c∈ (x,y) such that 26 R s = edge 1 edge 2 edge 3 edge 4 edge 5 edge 6 edge 7 edge 8 edge 9 edge 10 edge 11 edge 12 edge 13 edge 14 edge 15 edge 16 edge 17 edge 18 edge 19 edge 20 edge 21 edge 22 left boundary1000000000000000000000 x 1 0110000000000000000000 x 2 0001100000000000000000 x 3 0000011000000000000000 x 4 0000000110000000000000 x 5 0000000001100000000000 x 6 0000000000011000000000 x 7 0000000000000110000000 x 8 0000000000000001100000 x 9 0000000000000000011000 x 10 0000000000000000000110 right boundary0000000000000000000001 R r = edge 1 edge 2 edge 3 edge 4 edge 5 edge 6 edge 7 edge 8 edge 9 edge 10 edge 11 edge 12 edge 13 edge 14 edge 15 edge 16 edge 17 edge 18 edge 19 edge 20 edge 21 edge 22 left boundary0100000000000000000000 x 1 1001000000000000000000 x 2 0010010000000000000000 x 3 0000100100000000000000 x 4 0000001001000000000000 x 5 0000000010010000000000 x 6 0000000000100100000000 x 7 0000000000001001000000 x 8 0000000000000010010000 x 9 0000000000000000100100 x 10 0000000000000000001001 right boundary0000000000000000000010 f ′ (c) = f (y)− f (x) y− x (C10) If a function f has a continuous derivative f ′ on [x,y], by the mean value theorem, ||f (y)− f (x)||≤ C||y− x|| (C11) satisfies a Lipschitz continuous on [a,b] where C ≥ max|∇f (x)|, which means there exists a constant C that is greater than or equal to the maximum rate of change of function f in its domain. This constant C is called the Lipschitz constant, and this function f is called a Lipschitz function. Thus, the activation function would be chosen as a smooth and differentiable function across the domain to satisfy Lipschitz contin- uous. The Mish function Mish(x) [57] is used for all activation layers in the proposed model for non-monotonic, non-monotonic derivative, not saturated, and infinite continuity. tanh(x) = e x − e −x e x + e −x softplus(x) = ln(1 + e x ) Mish(x) = xtanh(softplus(x)) = x e ln(1+e x ) − e −ln(1+e x ) e ln(1+e x ) + e −ln(1+e x ) (C12) Mish ′ (x) = e x ω δ 2 (C13) where ω = 4(x + 1) + 4e 2x + e 3x + e x (4x + 6) and δ = (e x + 1) 2 + 1. Based on the chain rule, the composition of smooth functions is also smooth. The proposed model would be smooth after applying all smooth activation functions. For predicting the future state, velocity is chosen to be predicted first; then the position is inferred, rather than predicting position directly, based on the following analysis. When directly predicting the node position, an error term of the predicted position is denoted as ε. Then, v i t = x i t +ε−x i t−1 dt can be derived. Since the setting dt = 1e −4 , this would significantly increase the error of predicted position x i t to v i t with 1e 4 magnitude. Thus, the error would be carried in the next iteration for the velocity feature input as v i t+1 = v i t + ε∗ 1e 4 + F i t m i dt for the neural net- works, and this causes the predicted position to 27 be inaccurate with more errors. After more itera- tions, the results would be unreasonable, and the experiments can show the model would simulate unjustifiable results after around a few hundred iterations, which is consistent with the theory. The error can be reduced by predicting the velocity rather than the x-axis position in the model since fewer errors in the input of the model can bind tighter output of the model based on the Lipschitz continuous. Let denote ε t is a one-step error for the output of the proposed network with input ̃ x t−1 , ̃ x t , ̃ v t at t, where v t+1 is ground truth velocity value and ̃ v t+1 is the predicted velocity value. ε t = CN( ̃ x t−1 , ̃ x t , ̃ v t )− v t+1 = ̃ v t+1 − v t+1 (C14) Table C5 indicates that the model trained to directly predict positions tend to learn a short- horizon mapping from the current state to the next states, rather than learning the underlying physical movement patterns that govern how floe states evolve through time. The CN predicting position even fails to simulate across the whole time domain, thus, the simulation RMSE and PCC are not applicable. Empirically, this dis- tinction becomes most apparent in simulation. While position-predicting variants may achieve good one-step position errors but much worser one-step velocity errors , the errors compound rapidly when the model is run autoregressively, leading to substantially larger simulation error than velocity-predicting counterparts. For the predicted position at time t, the equation can be expanded into the ground truth position at time 1, since ̃ x 2 = x 1 + ̃ v 2 dt that is predicted by the CN with initial conditions x 0 ,x 1 . Also, ε max is denoted as the maximum one-step prediction error from the proposed model. ̃ x t+1 = ̃ x t + ̃ v t+1 dt = ̃ x t + (v t+1 + ε t )dt = ̃ x t−1 + (v t + ε t−1 )dt + (v t+1 + ε t )dt ... = ̃ x 2 + (v 3 +· + v t+1 )dt + (ε 2 +· + ε t )dt = x 1 + (v 2 +· + v t+1 )dt + (ε 1 +· + ε t )dt ≤ x 1 + (v 2 +· + v t+1 )dt + tε max dt (C15) Therefore, the model’s error would be linearly bounded by the maximum one-step prediction error, which is preferable to exponentially increas- ing the error. Appendix D Supplementary information: Skill scores Mean Squared Error (MSE) is a widely used sta- tistical measure for evaluating the performance of predictive models. It has been adopted to eval- uate the performance of the model, since MSE is preferred for its ability to heavily penalize larger errors and its suitability for mathematical optimization. MSE = 1 n n X i=1 (x i −y i ) 2 (D16) wherex is predicted value,y is the ground truth value, and n is the size of the sample. Root Mean Squared Error (RMSE) is com- puted by taking the square root of MSE. RMSE = √ MSE = v u u t ( 1 n ) n X i=1 (y i −x i ) 2 (D17) The pattern correlation coefficient is used to measure the overall strength of the relationship between the predicted trajectories and ground- truth trajectories. PCC = P n i=1 (x i − ̄ x)(y i − ̄ y) p P n i=1 (x i − ̄ x) 2 p P n i=1 (y i − ̄ y) 2 (D18) 28 Table C5: CN performance comparison with different predicting states ModelSimulation RMSE one-step position RMSE one-step velocity RMSEPCC CN predicting velocity2.3340240.00000600.0284100.934317 CN predicting position-0.00000033.022292- Different evaluation methods have different attributes and emphasize different traits of the performance. For a more comprehensive and accu- rate performance analysis, all mentioned tech- niques will be used to compensate for the limita- tions of using only a specific method. 29 (a) RMSE for velocity(b) RMSE for position Fig. E2: Training and validation loss in different epochs. Appendix E Supplementary information: Training configuration All models were trained with 10 million ground truth and prediction pairs data with a batch size of 100 using the Adam optimizer [59], which is widely used due to its computational efficiency and adaptive moment estimates that stabilize stochastic optimization in high-dimensional parameter spaces. The initial learning rate is set to 1× 10 −4 and an exponential learning-rate schedule is applied with a multiplicative decay factor (gamma) of 0.99 per epoch. This choice is motivated by the well-established sensitivity of neural training dynamics to the learning rate. Learning-rate decay is a practice in modern deep learning because it enables larger exploratory updates early in training while progressively shifting toward smaller, fine-grained updates as optimization approaches a solution. Moreover, theoretical and empirical analyses show that maintaining a comparatively larger learning rate early in training, while progressively shifting toward smaller updates, can have a regularization effect and can reduce premature memorization of spurious patterns, thereby improving generalization [86–89]. As shown in Figure E2, the training and validation losses for both position and velocity stabilize after approximately 180 epochs, with no evidence of divergence between the two curves, and training required approximately 16 hours on a single NVIDIA RTX 4090 GPU. This behavior is consistent with good generalization: overfitting is typically diagnosed when the training loss continues to decrease while the validation loss stops improving and begins to increase, producing a clear separation between training and validation loss trajectories [90, 91]. In the experiments, the absence of this characteristic divergence indicates that the selected optimization setting yields stable convergence without observable overfitting under the training duration. Appendix F Supplementary information: Visualization To enhance the interpretability of model performance and provide qualitative evidence of physical consis- tency, supplementary videos are provided with the simulation results. These visualizations are designed to convey critical aspects of the predicted floe dynamics, as the following: • No Unphysical Overlap or Boundary Breach: Throughout the simulations, individual sea ice floes are clearly shown to maintain spatial separation, never unphysical overlapping or passing through one another. This behavior demonstrates that the model correctly maintain physical constraints. Similarly, floes are consistently confined within the simulation domain, and no floe crosses the defined boundary, further validating the adherence to physical domain restrictions. 30 • Realistic Rebound Dynamics: When collisions occur, floes exhibit realistic elastic rebound behavior rather than unphysical sticking or interpenetration. The rebound trajectories are visually consistent with expected outcomes, as defined in the DEM framework. • Color-Coded Floe Identification: To support visual tracking, each floe is assigned a distinct color that remains fixed throughout the simulation. This coloring facilitates frame-by-frame analysis of individual floe motion and makes it easier to assess whether floes interact and move correctly. F.1 Dataset visualization Fig. F3: 10 sea floes dataset snapshot The video for dataset can be found at figshare[78]. This shows the generated dataset follows physical rules. 31 Fig. F4: 30 sea floes dataset snapshot 32 F.2 Prediction visualization Fig. F5: 10 sea floes prediction snapshot Figure F6 provides a qualitative snapshot of the predicted dynamics in the 30-floe experiment and illustrates two representative interaction modes captured by the model: floe–floe collisions and 33 floe–boundary collisions. For floe–floe interaction, the first two floes on the left demonstrate a clear approach–contact–separation sequence. Between t = 0.071 and t = 0.082, the two floes move closer as their distance decreases, and they make contact at t = 0.082. Immediately after contact, the subsequent frames show that the floes rebound and move in opposite directions, consistent with physical laws. For the floe and boundary interaction, it is visible for the first floe on the right. From t = 0.071 to t = 0.085, this floe translates toward the right boundary, reaches contact at t = 0.085, and then rebounds, reversing direction and moving away from the boundary thereafter. Together, these sequences demonstrate that the model produces physically plausible interaction patterns across both internal collisions and bound- ary interactions, supporting the qualitative validity of the learned collision dynamics in addition to the quantitative performance metrics reported. Fig. F6: 30 sea floes prediction snapshot The videos for prediction can be found at figshare[78]. These visualizations vividly illustrate the inter- actions modeled within the proposed model in different node size settings, effectively capturing various collision dynamics. Blue nodes represent individual sea ice floes, while black lines delineate the bound- aries of the simulation area. This visualization demonstrates the model’s capability to accurately simulate different scenarios: the collision between two floes, a single floe’s interaction with the boundary, and cases where no collisions exist. Each scenario is distinctly represented, showcasing the model’s comprehensive and precise handling of physical interactions in the dynamic environment. This accuracy substantiates that the proposed model’s predictions are visually feasible and do not contravene any physical laws. 34 F.3 Generalization visualization Fig. F7: 30 sea floes generalization snapshot The video for generalization can be found at figshare [78]. The visualization for generalization outside the ground truth data shows the proposed model’s capability of generalization outside the original time domain. 35 Fig. F8: 10 sea floes generalization snapshot 36