Non-Newtonian viscous fluid models with learned rheology accurately reproduce Lagrangian sea ice simulations

TL;DR

This work addresses the challenge of accurately modeling sea-ice rheology by learning a nonlinear, concentration-dependent viscosity directly from DEM velocity data. The authors represent the effective viscosity as a neural-network function $oldsymbol{ u}_{oldsymbol{ heta}}(|oldsymbol{ abla}oldsymbol{u}|, A)$ that respects isotropy and frame-indifference, and train it via PDE-constrained optimization using velocity data $oldsymbol{ ext{K}}$ derived from a DEM (SubZero). The learning proceeds in two steps to ensure well-posedness, enforcing nonnegativity and monotonicity, and yields a rheology that transitions from shear-thickening to shear-thinning as ice concentration increases, with viscosity changing by orders of magnitude for modest concentration shifts. The learned model generalizes to unseen forcing, time-dependent problems, and two-dimensional configurations, suggesting a scalable approach to data-driven continuum sea-ice models that better capture granular dynamics than existing Hibler-style formulations.

Abstract

Polar sea ice is crucial to Earth's climate system. Its dynamics also affect coastal communities, wildlife, and global shipping. Sea ice is typically modeled as a continuum fluid using a model proposed almost 50 years ago, which is moderately accurate for packed ice, but loses its predictive accuracy outside of the central ice pack. Discrete element methods (DEMs), which are commonly used for modeling granular media, offer an alternative by resolving the behavior of individual ice floes, including collisions, frictional contact, fracture, and ridging. However, DEMs are generally too costly for large-scale simulations. To address this, we present a framework for inferring rheological behavior from DEM velocity data. We characterize isotropic constitutive laws as scalar functions of the principal invariants of the strain-rate tensor. These functions are parameterized by neural networks trained on DEM data. By combining machine learning and finite element methods, we incorporate the governing partial differential equation (PDE) into the training, requiring to solve a PDE-constrained optimization problem for the network parameters. We focus on unidirectional parallel shear flows, which allow us to infer the effective shear viscosity. We find that, over a wide range of ice concentrations, the velocity fields observed in a complex sea ice DEM can be captured by a nonlinear rheology. Depending on the ice concentration, a shear-thinning or a shear-thickening behavior is observed. Moreover, the effective shear viscosity is found to increase by several orders of magnitude with changes as small as 5% in the sea ice concentration. We show that the learned rheology generalizes to different forcing scenarios, time-dependent problems, and settings in which compressibility is not a dominant factor.

Figures (9)

• Figure 1: (a) Diagram summarizing the training framework: with the DEM we compute the position $\boldsymbol{x}_i$ for floe $i$ with mass $m_i$ forced with ocean drag ($\boldsymbol{f}_o$), wind drag ($\boldsymbol{f}_w$), and collisions ($\boldsymbol{f}_{ij}$). We extract horizontal velocities $u_k$ averaged over each cell $k$ and minimize the mismatch $\mathcal{J}_v$ with the solution $u(\boldsymbol{\theta})$ to the continuum model, which is defined in terms of an NN-based rheology $\psi_{\boldsymbol{\theta}}$ that also depends on the sea ice concentration $A$. (b) Ice-floe field simulated in SubZero. For generating the training data, we drive the floes with the hat-shaped horizontal ocean velocity profile shown in blue and zero wind velocity. We average over the elongated cells with black dashed edges. Inset: realistic power-law floe size distribution satisfied by ice floes for areas between $10^{-4}L^2$ and $2\times 10^{-3}L^2$. (c) NN-based shear stress model inferred from training are shown using solid lines. Markers represent DEM data. (d) Two steady velocity profiles from training dataset. Sea ice NN model in red line and DEM data in markers (markers of same type and color in panel (c), circled, correspond with same steady states); blue lines represent ocean velocity. (e) Effective shear viscosity $\psi_{\boldsymbol{\theta}}$ inferred with training.
• Figure 2: Steady states computed with the DEM for $A = 0.85$ and $U_o = 1m\per s$ with five randomized floe-field initializations (runs 1-5). Initial and final floe fields for run 1 (a)-(b) and run 2 (c)-(d). (e) Velocity profiles computed with runs 1-5.
• Figure 3: Comparison between velocity profiles computed with learned continuum model (black lines) and DEM (markers) for training. The steady states are computed with the DEM for concentrations $A = 0.85$ (a-e) and $A = 0.95$ (f-j) and maximum ocean velocities $U_o$ between 0.05 and 2 m/s. The rheology of the continuum models is inferred in a two step optimization process: in step 1, we minimize the stress-strain misfit $\mathcal{J}_s$ (red dotted line) and, in step 2, using the fit found in step 1 as initial guess, the velocity misfit $\mathcal{J}_v$ (black line), yielding the final continuum model.
• Figure 4: Horizontal velocity profiles for a steady one-dimensional problem used for testing the model's generalizability. We plot the velocity computed with the DEM (markers) and with the learned model (black lines). The red line represents the wind velocity profile, which is given by a cosine profile $u_w(y) = U_w/2 (1 - \cos(2\pi y/L))$ with amplitude $U_w = 20m\per s$.
• Figure 5: Unsteady test with concentration set to $A=0.875$. (a) Normalized amplitude of the wind velocity. (b) Difference between DEM and FEM model (black) and between DEM and free drift velocity (purple) (c-e) Comparison between velocity fields computed with our continuum model (black lines), the DEM (circles) and the free drift velocity (purple) for the one-dimensional unsteady test problem, at three different time instants. The velocities are non-dimensionalized with the equilibrium velocity $U^\ast = \sqrt{C_a\rho_a/(C_o\rho_o)}\,U_w$. (f) Shape of wind velocity profile when $U_w = 1$. (g) Shape of ocean velocity profile.