WIce-FOAM 1.0: coupled dynamic and thermodynamic modelling of heterogeneous sea ice and waves using OpenFOAM-v2306

We present WIce-FOAM 1.0, a numerical model built on OpenFOAM that couples the dynamics and thermodynamics of heterogeneous sea ice to analyse waves' response in marginal ice zone regions composed of consolidated ice floes and interstitial grease ice. The model represents prototypical conditions on the 5 km scale, where each 10 m grid cell classified as ice floe or grease ice may contain both ice types, but are predominantly occupied by one. Our model aims to study the mean shear viscosity of heterogeneous sea ice to bridge the gap with larger-scale ocean-sea ice models in which sub-grid details and wave effects are neglected. We tested the model in the Southern Ocean using a realistic sea-ice field from a SAR satellite image and complemented our analysis by idealised simulations. The thermodynamic model was coupled online to optimize the stiffness of the process scales and to explicitly account for the distinct characteristics of different ice types. We first investigated the dynamic response of sea ice to one-way wave forcing across a range of wave periods and directions. The results show that the domain-averaged sea-ice viscosity is scale invariant from approximately 800 m to 5 km and is primarily governed by the relative proportion of ice floes to grease ice, with less sensitivity to wave periods and directions. While the wave direction affects the local strain rate and viscosity, and the presence and orientation of narrow connections between the larger ice floes significantly influence the mean viscosity, these effects do not break the observed scale invariance. Finally, we demonstrate that, despite the different time scales, the mean viscosity responds nonlinearly to the inclusion of thermodynamic sea-ice growth. This model represents a first step towards a mechanistic understanding and description of heterogeneous sea ice, which is common in the Antarctic and is increasing in the warming Arctic. It can be used to design field experiments and to derive parametrisations of waves-in-ice response for large-scale sea-ice models.

1 Introduction

Antarctica has experienced relatively stable sea ice levels over recent decades and an apparent change of regime since 2017 (Diamond et al., 2024; Gilbert and Holmes, 2024; Wang et al., 2023). To understand the intricate response of Antarctic sea ice to climate change (Golden et al., 2020), new-generation models should address variability in sea ice across smaller spatial and temporal scales (Iovino et al., 2022; Selivanova et al., 2024).

Sea ice is a dynamic heterogeneous medium, which is typically described as a mixture of ice constituents and open water areas (Barthélemy et al., 2016). However, away from the thicker pack ice regions, and especially in the Antarctic marginal ice zone, each type of ice has distinct material properties with pronounced thermal gradients (Tersigni et al., 2023).

During the winter season, the Antarctic sea ice cover expands over the Southern Ocean and reaches its maximum in September and October. In this phase, the ice constituents vary from newly formed ice in open water, which under the action of ocean waves transitions from a thin layer of frazil ice into a thicker, slushy layer known as grease ice, and eventually grows into more solid-like floes known as pancake ice floes (Wadhams et al., 2018; Nose et al., 2021). These pancake ice floes merge into aggregated ice floes and a more consolidated sea ice cover in the interior when the external forcing from ocean waves subsides.

During sustained warming conditions, ice-albedo feedback accelerates the melting of sea ice (Golden et al., 2020; Squire, 2020). Consequently, the heterogeneous sea ice cover weakens and disaggregates, and becomes more responsive to wind and ocean currents, also allowing for the propagation of waves deeper into the sea ice (Iovino et al., 2022; Squire, 2020; Thomson and Rogers, 2014). This creates a feedback loop, where larger ice floes break into smaller, mobile floes, increasing the transfer of momentum and energy into the sea ice (Thomson, 2022; Asplin et al., 2012). This process hinders the formation of consolidated sea ice, keeping the sea ice cover in a more heterogeneous condition (Day et al., 2024).

Most existing large-scale sea ice models adopt the continuum approach originally proposed by Hibler (1979) and Thorndike et al. (1975), which is generally considered valid over spatial scales of hundreds of kilometres (Mehlmann and Richter, 2017). Advancements in computing power in recent decades have enabled a shift toward higher-resolution models. This increased computational capacity has sparked renewed interest in exploring the properties of sea ice on spatial scales of tens of kilometres or less (Zhang, 2021), allowing us to apply and assess the continuum approach at finer resolutions.

Persistently challenging in continuum models is the necessity of integrating a suitable rheology for sea ice, i.e. the relationship between sea ice stress and deformation. Given the diverse spatial and temporal characteristics of the sea ice cover, there is no straightforward approach to accurately model sea ice dynamics using an effective large-scale rheology that accounts for all relevant processes (Åström et al., 2024). The need for introducing appropriate rheology for the various ice types has recently been highlighted (Herman, 2016; Golden et al., 2020; Skatulla et al., 2022). Different sea ice thicknesses and types experience specific growth and melting rates, i.e. thinner ice experiences faster growth and melting compared to thicker ice. Thinner ice is also more susceptible to mechanical deformation (Golden et al., 2020). Moreover, the presence and properties of frazil and grease ice and the concurrent action of waves are barely addressed (Dumont, 2022), largely due to limited research on frazil ice in terms of field observations (Paul et al., 2021).

Small-scale models are useful for informing the development of parametrisations for large-scale sea ice models. Marquart et al. (2021, 2023) introduced a two-dimensional small-scale computational fluid dynamics (CFD) model implemented in OpenFOAM, designed to represent heterogeneous sea ice conditions on the metre scale. The model distinguishes between two distinct ice types: ice floes and interstitial grease ice, as shown in Fig. 1a, each characterised by unique material properties, and captures sea ice dynamics under wavy conditions. Thermodynamic effects were excluded, as the simulations focused on short time scales of less than a minute.

Figure 1 Comparison of the heterogeneous sea ice cover at two spatial scales: (a) the metre scale (10–100 m²), based on a photograph taken during the winter cruise in July 2017, where individual ice floes and grease ice are clearly distinguishable (Marquart et al., 2021, 2023), and (b) the kilometre scale (1–10 km²), based on a photograph taken during the SCALE winter cruise in July 2022, where regions identified as either ice floes or grease ice primarily consist of one type, though both may be present.

Here, we extend the model developed by Marquart et al. (2021, 2023) by expanding the spatial domain to the kilometre scale. At this scale, individual ice floes and grease ice patches can no longer be explicitly resolved, as illustrated in Fig. 1b. Nonetheless, representing the heterogeneous nature of the sea ice cover remains essential. To accommodate this, we reinterpret the definitions of “ice floe” and “grease ice” within the model: a region is classified as “ice floe” if it predominantly consists of ice floes, even if grease ice is present. Conversely, the same applies to grease ice.

Additionally, at this kilometre scale, we incorporate thermodynamic processes, which require extending the time window to the order of days. The aim is to capture the evolution of key sea ice variables, such as thickness and viscosity, across space and time, which are particularly relevant for upscaling to large-scale, global sea ice models.

The remainder of this work is structured as follows. Section 2.1 describes the dynamic model and its implementation within the OpenFOAM framework. Section 2.2 introduces the thermodynamic model, based on the formulation by Tedesco et al. (2009) and the frazil ice production equation from Haarpaintner et al. (2001). The coupling approach between the dynamics and thermodynamics models is presented in Sect. 2.3. Section 2.4 outlines the model configuration used in the simulations, and results are discussed in Sect. 3. A detailed discussion of the key findings is provided in Sect. 4. Finally, Sect. 5 summarises the work.

2 Methodology

2.1 Dynamic OpenFOAM model

OpenFOAM is an open-source computational fluid dynamics (CFD) software toolbox used to solve a wide range of problems, including complex fluid flows, heat transfer and solid mechanics (Medina et al., 2015). The software is highly customizable and extensively used in both research and industry.

Our proposed two-dimensional continuum framework, WIce-FOAM 1.0, is implemented using OpenFOAM-v2306 and formulated in the xy-plane using the finite volume method (FVM) and the volume of fluid (VOF) method. The FVM discretizes the domain into a finite number of control volumes, or cells, and applies conservation laws to each cell (Ferziger et al., 2019). The VOF method is a numerical FVM technique to describe the interface between two immiscible and incompressible materials. In our model, we use the IsoAdvector method as the VOF implementation (Roenby et al., 2016). Both methods are key to the model and are important for coupling the dynamics and thermodynamics models in Sect. 2.3.

Within this framework, we consider a domain with 100 % sea ice concentration and investigate the impact of sea ice heterogeneity on both the dynamic response to wave forcing (in a one-way coupled setup) and on thermodynamic processes. Phase transitions between different ice types are not included at this stage. We model a heterogeneous sea ice cover composed of ice floes and grease ice. Each cell is classified as either “ice floe” or “grease ice” according to the dominant ice type within that cell, determined from the initial condition field derived from a synthetic aperture radar (SAR) image (Sect. 2.4), even though both types may be present within a single cell in the real system.

Grease ice is treated as a highly viscous fluid that dissipates energy during ice floe collisions. As a result, interactions between ice floes are represented as continuous, churning contact of varying intensity rather than brief, forceful impacts. This justifies the exclusion of ice floe failure and fracture due to floe-floe collisions. Other potential fracture mechanisms are also excluded, such as those driven by ice convergence, shear deformation, or wave-induced bending. Lastly, we focus our analysis on one-way wave-ice interactions: a harmonic propagating wave forcing is imposed on the sea ice domain, which excludes wave dissipation effects.

2.1.1 Momentum balance equation

The momentum transport of incompressible sea ice is formulated in the momentum balance equation as

$$\begin{array}{}\text{(1)}& \mathit{\rho }\left({\displaystyle \frac{\partial h\mathit{U}}{\partial t}}+(\mathit{U}\cdot \mathrm{\nabla })h\mathit{U}\right)={\mathit{\tau }}_{\mathrm{w}}+{\mathit{\tau }}_{\mathrm{a}}+\mathrm{\nabla }\cdot \mathit{\sigma },\end{array}$$

where U represents the two-dimensional sea ice velocity vector, h is the sea ice thickness, and t is time. External forces, in the form of the in-plane oceanic wave and wind stresses, act on the sea ice cover and are denoted as τ_w and τ_a, respectively. The Cauchy stress tensor, σ, characterizes the stress state of the sea ice and differs by ice type, denoted σ_f for ice floes and σ_g for grease ice.

The sea ice mass per unit area, m, is

$$\begin{array}{}\text{(2)}& m=\mathit{\rho }h,\end{array}$$

where ρ denotes the sea ice density. The evolution of sea ice thickness (Hibler, 1979; Hutchings et al., 2004), h, is given as

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial h}{\partial t}}+\mathrm{\nabla }\cdot \left(h\mathit{U}\right)=\mathrm{0},\end{array}$$

which extends the model proposed by Marquart et al. (2021, 2023), wherein a constant sea ice thickness was assumed.

The external in-plane oceanic wave stress, denoted as τ_w, is derived from the linear wave theory (Holthuijsen, 2010; Herman, 2018). This stress is characterized by the presence of two distinct components

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{\mathrm{w}}={\mathit{\tau }}_{\mathrm{sd}}+{\mathit{\tau }}_{\mathrm{fk}},\end{array}$$

where the skin drag, τ_sd, represents the viscous effects across the ice-ocean interface, and the Froude–Krylov force, τ_fk, results from the wave pressure field acting on the submerged surface of the ice floes.

The quadratic skin drag (Hutchings, 2000) is given as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\tau }}_{\mathrm{sd}}& ={\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{w}}|{\mathit{U}}_{\mathrm{w}}-\mathit{U}|\left(\left({\mathit{U}}_{\mathrm{w}}-\mathit{U}\right)\cos {\mathit{\theta }}_{\mathrm{w}}+\left({\mathit{U}}_{\mathrm{w}}-\mathit{U}\right)\right.\\ & \left.\times \mathit{k}\sin {\mathit{\theta }}_{\mathrm{w}}\right),\end{array}\end{array}$$

where ρ_w represents the water density, and θ_w is the ice-ocean turning angle. The unit normal vector to the ice surface is denoted by k, and C_w is the ice-ocean drag coefficient, which varies between the two ice types. The orbital wave velocity of the water, U_w, for monochromatic waves is defined as

$$\begin{array}{}\text{(6)}& {\mathit{U}}_{\mathrm{w}}=\left(\begin{array}{c}{U}_{w_x}\\ U_{w_y}\\ U_{w_z}\end{array}\right)=\left(\begin{array}{c}a\mathit{\omega }\sin \left(\mathit{\omega }t-k_{x}x\right)\\ a\mathit{\omega }\sin \left(\mathit{\omega }t-k_{y}y\right)\\ a\mathit{\omega }\cos \left(\mathit{\omega }t-k_{x}x-k_{y}y\right)\end{array}\right),\end{array}$$

where x, y and z represent the Cartesian coordinates within the sea ice domain, with z denoting the vertical direction. The model is two-dimensional, the vertical component of the wave velocity does not enter the governing momentum balance equation, but is included here for completeness. The parameters a, ω, and k correspond to the wave amplitude, wave frequency, and wave number, respectively. The wave frequency, ω, and wave number, k, are given by $\mathit{\omega }=\mathrm{2}\mathit{\pi }/T$ and $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$, derived from the wave period, T, and wavelength, λ, using the deep water dispersion relation, ω²=gk, where g denotes the gravitational acceleration. The wave numbers in x and y direction are formulated as

$$\begin{array}{}\text{(7)}& k_x=k\cos \left({\mathit{\theta }}_{\mathrm{wa}}\right),k_y=k\sin \left({\mathit{\theta }}_{\mathrm{wa}}\right),\end{array}$$

where θ_wa denotes the wave direction angle, measured relative to the x-axis and defined as positive in the counter-clockwise direction. A value of θ_wa=0° corresponds to wave propagation exclusively along the x-axis.

The in-plane form drag, acting on the ice floe circumference due to velocity differences between the floes and surrounding grease ice, is implicitly captured by the continuum approach, which includes both ice constituents and enforces velocity continuity at the interface.

The Froude–Krylov force, τ_fk, acts in the basal plane (the lower ice surface in contact with the ocean parallel to the xy-plane) and represents the horizontal surge force generated by the wave-induced pressure (Herman, 2018) at the interface between the ice floe and the water. This is:

$$\begin{array}{}\text{(8)}& {\mathit{\tau }}_{\mathrm{fk}}=-\underset{h_{\mathrm{w}}}{\int }p\mathit{n}\mathrm{d}z,\end{array}$$

where h_w is the height of the submerged portion of the ice floe thickness, and n is the outward pointing unit vector, normal to the ice floe circumference. The wave-induced pressure, denoted as p, is

$$\begin{array}{}\text{(9)}& p={\mathit{\rho }}_{\mathrm{w}}ga\sin (\mathit{\omega }t-kx).\end{array}$$

The atmospheric forcing is represented by a wind stress (Hutchings, 2000), τ_a, applied to the upper (apical) surface of the ice exposed to the atmosphere (also parallel to the xy-plane), and is expressed as

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{a}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{a}}\left|{\mathit{U}}_{\mathrm{a}}\right|\left({\mathit{U}}_{\mathrm{a}}\cos {\mathit{\theta }}_{\mathrm{a}}+\mathit{k}\times {\mathit{U}}_{\mathrm{a}}\sin {\mathit{\theta }}_{\mathrm{a}}\right),\end{array}$$

where U_a is the wind velocity, ρ_a denotes the air density, C_a is the ice-air drag coefficient, and θ_a is the wind turning angle. Note that this term is analogous to the ocean drag formulation. However, since sea ice drift velocities are typically much smaller than wind speeds, the difference is generally negligible (Leppäranta, 2011).

2.1.2 Sea ice rheology

The sea ice rheology in the model defines the relationship between internal ice stresses, denoted by σ, and ice deformation, expressed in terms of the strain rate, $\dot{\mathit{ϵ}}$. In accordance with the infinitesimal small strain theory, the strain rate tensor can be expressed in relation to the sea ice velocity gradient, ∇U, as

$$\begin{array}{}\text{(11)}& \dot{\mathit{ϵ}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{\nabla }\mathit{U}+(\mathbf{\nabla }\mathit{U}{)}^T\right).\end{array}$$

Each component is characterized by its own rheology, with ice floes displaying solid-like behaviour and grease ice exhibiting fluid-like behaviour. Marquart et al. (2021, 2023) used a “Hookean-like” flow rule to describe the constitutive law for the solid-like behaviour of ice floes. While the model effectively reproduced solid-like behaviour, it did not account for elastic unloading.

In this study, we replace the “Hookean-like” ice floe rheology with one that incorporates viscous behaviour. Therefore, the ice floes follow the viscous-plastic rheology proposed by Hibler (1979):

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\mathrm{f}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{f}}\dot{\mathit{ϵ}}+\mathit{I}\left(\left({\mathit{\zeta }}_{\mathrm{f}}-{\mathit{\eta }}_{\mathrm{f}}\right)\mathrm{tr}\left(\dot{\mathit{ϵ}}\right)-{\displaystyle \frac{P_{\mathrm{f}}}{\mathrm{2}}}\right),\end{array}$$

where σ_f denotes the Cauchy stress tensor for ice floes. The shear and bulk viscosities of ice floes are represented by η_f and ζ_f, respectively. The internal ice floe strength is P_f and I represents the identity tensor. The two strain rate-dependent viscosities are coupled via

$$\begin{array}{}\text{(13)}& {\mathit{\zeta }}_{\mathrm{f}}={\displaystyle \frac{P_{\mathrm{f}}}{\mathrm{2}\mathrm{\Delta }}},{\mathit{\eta }}_{\mathrm{f}}={\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{f}}}{e^{\mathrm{2}}}},\end{array}$$

where e indicates the eccentricity, the ratio between the in-plane principal axes of the elliptical yield curve (Hibler, 1979). The effective strain rate, Δ, is

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }=\sqrt{\left({\dot{\mathit{ϵ}}}_{\mathrm{11}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}^{\mathrm{2}}\right)\left(\mathrm{1}+e^{-\mathrm{2}}\right)+\mathrm{4}{e}^{-\mathrm{2}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{22}}\left(\mathrm{1}-e^{-\mathrm{2}}\right)},\end{array}$$

where ${\dot{\mathit{ϵ}}}_{\mathrm{11}}$, ${\dot{\mathit{ϵ}}}_{\mathrm{22}}$ and ${\dot{\mathit{ϵ}}}_{\mathrm{12}}$ denote the Cartesian components of the symmetric strain rate tensor. As the strain rate approaches zero, the viscosity becomes unbounded. To prevent this singularity, we impose a lower bound on the effective strain rate, ${\mathrm{\Delta }}_{\min }=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹. The modified internal ice floe strength, P_f, can be written as

$$\begin{array}{}\text{(15)}& P_{\mathrm{f}}=P_{\mathrm{f}}^{*}h,\end{array}$$

where $P_{\mathrm{f}}^{*}$ represents an empirical constant. Equation (15) differs from Hibler (1979) by excluding the compactness parameter, which is unnecessary in our model due to the explicit treatment of the two ice phases.

In Marquart et al. (2021, 2023) grease ice was governed by a viscous-plastic (VP) material law, which is similar to the rheology developed by Hibler (1979) and Thorndike et al. (1975). However, Marquart et al. (2021, 2023) observed a singularity in the viscosity and strain rate of grease ice associated with the passage of waves. This singularity is characterized by locally very high viscosity values, linked to strain rate values approaching zero. To address unnatural behaviour due to this singularity, we present here a revised rheology for the grease ice.

In this study, grease ice is assumed to behave as an incompressible fluid, based on literature indicating that it is primarily composed of water (Smedsrud, 2011; Mackie et al., 2020). Consequently, its rheology is represented as that of an incompressible, non-Newtonian viscous fluid (Newyear and Martin, 1997; Eberhard et al., 2019):

$$\begin{array}{}\text{(16)}& {\mathit{\sigma }}_{\mathrm{g}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{g}}\dot{\mathit{ϵ}},\end{array}$$

where the Cauchy stress tensor for grease ice is indicated as σ_g. The shear viscosity, η_g, follows the Cross model incorporating shear thinning, wherein viscosity is constrained for extremely high and low strain rate values (Hauswirth et al., 2020; Eberhard et al., 2019; Galindo-Rosales et al., 2010). The shear viscosity, η_g, is

$$\begin{array}{}\text{(17)}& {\mathit{\eta }}_{\mathrm{g}}={\mathit{\eta }}_{\mathrm{\infty }}+{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{0}}-{\mathit{\eta }}_{\mathrm{\infty }}}{\mathrm{1}+\left(C\right|\dot{\mathit{\gamma }}|{)}^m}},\end{array}$$

where C represents the Cross time constant, $\left|\dot{\mathit{\gamma }}\right|=\sqrt{\mathrm{2}\dot{\mathit{ϵ}}:\dot{\mathit{ϵ}}}$ the shear rate magnitude, and m the degree of shear thinning. Inertia is likely to contribute to the effective viscosity of grease ice (de Carolis et al., 2005). Therefore, we assume that η_∞ and η₀ are thickness-dependent infinite and zero shear viscosities of grease ice, and follow a power law:

$$\begin{array}{}\text{(18)}& {\mathit{\eta }}_{\mathrm{0}}={\mathit{\eta }}_{{\mathrm{0}}_{\mathrm{c}}}{\left({\displaystyle \frac{h}{h_{{\mathrm{g}}_{\mathrm{0}}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{0}}},{\mathit{\eta }}_{\mathrm{\infty }}={\mathit{\eta }}_{{\mathrm{\infty }}_{\mathrm{c}}}{\left({\displaystyle \frac{h}{h_{{\mathrm{g}}_{\mathrm{0}}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{\infty }}},\end{array}$$

where ${\mathit{\eta }}_{{\mathrm{\infty }}_{\mathrm{c}}}$ and ${\mathit{\eta }}_{{\mathrm{0}}_{\mathrm{c}}}$ are the reference infinite and zero shear viscosities of grease ice, respectively, at the reference thickness of grease ice, $h_{{\mathrm{g}}_{\mathrm{0}}}$. The power law exponents α_∞ and α₀ describe how η_∞ and η₀ change with thickness.

Parameter values associated with the dynamic model are summarized in Table 1, which is organised into three groups, namely ice floe rheology parameters, grease ice rheology parameters, and wave-related parameters. Ice floe rheology parameters are based on the model by Hibler (1979) and related formulations (e.g. Mehlmann and Richter, 2017). Notably, the limit of the effective strain rate (Leppäranta and Hibler, 1985) is deliberately reduced by one order of magnitude to ensure numerical stability. Although this adjustment affects the solidity of the ice floes by modifying their viscosity, the overall simulation results remain largely unaffected. Values for the grease ice rheology are derived from literature sources such as Paul et al. (2021), and further refined through empirical tuning via iterative simulations. Wave-related parameters include the ice-ocean turning angle, set to zero to reflect that the water drag on the ice acts purely along the flow direction. Drag coefficient values associated with wave characteristics are taken from Smedsrud (2011) and Alberello et al. (2020). Finally, a range of wave parameters is selected to conduct a sensitivity analysis in Sect. 3.

Table 1 Parameters related to the dynamic model.

2.2 Thermodynamic model

As with the rheologies of ice floes and grease ice, discussed in Sect. 2.1.2, this study applies two distinct thermodynamic models to the two ice types.

For the ice floes, the one-dimensional thermodynamic model in the z-direction, developed by Tedesco et al. (2009), is applied to OpenFOAM cells associated with ice floes to simulate thermodynamic variations in snow and ice thickness. The model accounts for multiple layers of sea ice (including columnar ice, snow ice, and superimposed ice) and/or snow, as well as surface heat fluxes. These layers are assumed to be in thermal equilibrium, with interface temperatures determined by the continuity of heat fluxes:

$$\begin{array}{}\text{(19)}& {\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{i}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial z}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left[I_{\mathrm{pen}}\left(z\right)\right],\end{array}$$

where ρ is the density, c the specific heat and T the temperature of the layers (the subscript “i” denotes either snow or the sea ice layers). The thermal conductivity is defined by K, and I_pen is the flux of penetrating solar radiation through each layer. The vertical coordinate is defined as positive in the downward direction, where z=0 corresponds to the top surface. The penetrating radiation can be written as

$$\begin{array}{}\text{(20)}& I_{\mathrm{pen}}\left(z\right)=I_{\mathrm{0}}\exp \left(-{\mathit{\kappa }}_{\mathrm{i}}z\right),\end{array}$$

where I₀ indicates the penetrating solar flux at the top ice/snow surface and κ_i the extinction coefficient. The new temperatures for each layer are calculated using a finite-difference formulation of Eq. (19), with the calculation being performed through an iterative process. Growth and melting rates are determined from expressions for the enthalpy of snow and sea ice. The enthalpy of snow (or fresh ice) per unit volume (Hunke et al., 2015), q_s, is

$$\begin{array}{}\text{(21)}& q_{\mathrm{s}}\left(T\right)=-{\mathit{\rho }}_{\mathrm{s}}\left(-c_{\mathrm{0}}T+L_{\mathrm{0}}\right),\end{array}$$

where ρ_s is the snow density, c₀ the specific heat of fresh ice, and L₀ the latent heat of fusion of fresh columnar ice. The enthalpy of sea ice is not straightforward due to the presence of brine pockets, where salinity varies inversely with temperature. Assuming a predetermined value for salinity, a direct relationship between temperature and the enthalpy of sea ice per unit volume can be obtained:

$$\begin{array}{}\text{(22)}& q_{\mathrm{i}}\left(T\right)=-{\mathit{\rho }}_{\mathrm{i}}\left(c_{\mathrm{0}}\left(T_{\mathrm{m}}-T\right)+L_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{m}}}{T}}\right)-c_{\mathrm{w}}{T}_{\mathrm{m}}\right),\end{array}$$

where ρ_i is the sea ice density, T_m the temperature at which the ice is completely melted, and c_w the specific heat of seawater. Once the enthalpy is obtained, the surface growth/melting is:

$$\begin{array}{}\text{(23)}& q_{\mathrm{s},\mathrm{i}}\mathrm{\Delta }{h}_{\mathrm{s},\mathrm{i}}=\left\{\begin{array}{ll}\left(F_{\mathrm{0}}-F_{\mathrm{ct}}\right)\mathrm{\Delta }t& \mathrm{if}{F}_{\mathrm{0}}>F_{\mathrm{ct}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.,\end{array}$$

where Δh_s,i represents the thickness change in snow or sea ice, respectively, over the time step Δt. The net surface heat flux from the atmosphere to the ice, F₀, represents the resultant of all fluxes included in the thermodynamic model (Hunke et al., 2015). These comprise the sensible heat flux, latent heat flux, incoming and outgoing long wave radiation, and incoming short wave radiation. F_ct denotes the conductive flux from the top surface to the interior of the ice.

The growth and melting occurring at the bottom layer of ice is

$$\begin{array}{}\text{(24)}& q_{\mathrm{i}}\mathrm{\Delta }{h}_{\mathrm{i}}=\left(F_{\mathrm{cb}}-F_{\mathrm{bot}}\right)\mathrm{\Delta }t,\end{array}$$

where F_cb is the conductive heat flux at the bottom surface and F_bot is the net downward heat flux from the ice to the ocean.

For the cells in OpenFOAM associated with grease ice, a different thermodynamic approach is applied. We implemented the frazil ice production rate proposed by Haarpaintner et al. (2001) in the thermodynamic model (Tedesco et al., 2009):

$$\begin{array}{}\text{(25)}& \mathrm{\Delta }{h}_{\mathrm{f}}={\displaystyle \frac{F_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{f}}{L}_{\mathrm{f}}}}\mathrm{\Delta }t,\end{array}$$

where the rate of change of frazil ice thickness depends on the net surface heat flux, F₀, a constant value of frazil ice density, ρ_f, and latent heat of fusion specific for frazil ice, L_f.

In the thermodynamic model (Tedesco et al., 2009), most variables, including changes in snow and sea ice thickness, depend on the state at the previous time step, thereby necessitating the tracking of sea ice variables across time steps. This temporal tracking is particularly important for the coupling approach discussed in Sect. 2.3. The model is forced with atmospheric input variables for the selected location, including total cloud cover fraction, specific humidity of air and surface, air temperature, wind velocity components, downward surface solar radiation flux and the mean total precipitation rate, all sourced from ERA5 (Hersbach et al., 2023a, b), although alternative datasets can also be used.

Figure 2 Results from the thermodynamic model, illustrating (a) the thickness evolution from June to November 2022, and (b) a comparison of snow and sea ice thickness evolution with frazil ice development in July 2022.

Thickness evolution results from the thermodynamic model, based on the formulations by Tedesco et al. (2009), are presented in Fig. 2a, which depicts a complete seasonal cycle within 2022. The evolution of snow and sea ice thickness are shown in blue and red, respectively. Figure 2b provides a detailed view of July, comparing the evolution of sea ice thickness with frazil ice thickness, derived from Haarpaintner et al. (2001). Figure 2 shows that the thickness of both sea ice and snow is zero, i.e. no ice, during the summer months as the air temperature exceeds the threshold for ice formation and growth. Ice begins to form and its thickness to increase in winter (June) when the air temperature drops. Note that in the model, snow growth commences only once the minimum snow threshold ($h_{\mathrm{s},\min }=\mathrm{0.02}$ m) has been exceeded, see Fig. 2b. Sea ice thickness returns to zero in spring, after the sea ice melting. As expected, frazil ice exhibits a higher growth rate than sea ice. This is primarily because the net surface heat flux, F₀, increases for thinner ice, accelerating the rate of thickness change. Additionally, the lower latent heat of fusion of frazil ice compared to ice floes further enhances its growth. We observe that the thickness of frazil ice exceeds that of sea ice. While this outcome is physically unrealistic due to the expected phase transition, it does not affect our study, as our simulation results in Sect. 3 will be limited to a one-day period rather than an entire seasonal cycle.

The main parameter values associated with the thermodynamic model are summarized in Table 2. Only the constant parameters in the equations presented in this work are included; all other equations, parameters, as well as initial conditions, are available in the supplementary code.

Table 2 Parameters related to the thermodynamic model.

2.3 Coupling between the dynamics and thermodynamics models

Thermodynamic contributions can be seamlessly integrated into the dynamic model within OpenFOAM by manipulating the thickness variable and utilizing the VOF approach in the dynamic model to explicitly differentiate between ice types. The coupling between dynamics and thermodynamics is achieved by alternately running the dynamics and thermodynamics models within a for-loop.

The main challenge in the coupling arises from integrating both models that operate on distinct temporal scales, for which conventional up- or down-sampling methods to achieve a common timescale are considered impractical. The dynamic model represents a fast process, capturing the interaction between sea ice and a harmonically propagating wave. To accurately resolve the wave characteristics, the dynamic model requires a time step shorter than the wave period, typically in the order of seconds. By contrast, the thermodynamic model evolves on a much slower timescale, as processes of sea-ice melting and growth occur gradually over hours rather than seconds, and a time step of a few minutes is sufficient. To further improve computational efficiency, a well-established approach based on thickness categories is adopted (Sun and Solomon, 2024). In this method, sea ice is divided into a discrete set of thickness ranges, each representing ice of similar physical properties. The cells within the xy-plane are then grouped accordingly, and the average thickness for each category is computed. These category-averaged values are then used as inputs to the thermodynamic model, which is executed once for each thickness range, plus an additional run for grease ice, thereby reducing the number of required model runs.

Figure 3 Schematic of the coupling between the dynamics and thermodynamics models illustrating the first two iterations in the for-loop.

We note that the motion of sea ice subjected to a harmonically propagating wave exhibits periodic behaviour (Marquart et al., 2023). Therefore, we can assume that the dynamic response becomes periodic after one full wave period. Figure 3 illustrates a schematic of the coupling approach, showcasing the first two iterations in the for-loop. The simulation begins with a spin-up phase involving only the dynamics, during which the system evolves toward equilibrium under periodic wave forcing. Equilibrium for the dynamic component is considered reached when oscillations in sea ice velocity repeat over one wave period without any net change in velocity. Typically, the required spin-up time corresponds to approximately four to five wave periods (denoted t_n). After the spin-up phase, coupling between dynamics and thermodynamics begins at t_n, initiated by a for-loop, with the two models being executed alternately.

The first iteration begins with the execution of the thermodynamic model over one thermodynamic interval, defined as the time between T₁ and T₂. The resulting thickness change, $\mathrm{\Delta }{h}_{T_{\mathrm{1}}}$, is initially stored. Subsequently, the dynamic model resumes from the end of the spin-up phase and runs for one dynamic interval – equivalent to a single wave period (indicated in boldface in Fig. 3) – as successive, identical wave periods within the thermodynamic interval are not simulated (represented by dashed oscillations and marked with a thick cross in Fig. 3). The thickness in the dynamic model is then updated using the stored output from the thermodynamic model (indicated by $\mathrm{\Delta }{h}_{T_{\mathrm{1}}}$ in Fig. 3), resulting in a change in sea ice velocity. This concludes the first iteration. The second iteration begins with the thermodynamic model advancing over the next thermodynamic interval, from T₂ to T₃, and the sequence repeats until the for-loop is completed.

One of the objectives of the present work is to demonstrate that coupling between the dynamics and thermodynamics models is necessary. This is assessed by evaluating the linearity of their relationship with respect to sea ice thickness and sea ice viscosity using the following equations:

$$\begin{array}{}\text{(26)}& \underset{\mathrm{decoupled}}{\underbrace{{\stackrel{\mathrm{‾}}{h}}_{\mathrm{D}}+\sum \mathrm{\Delta }{h}_T}}=\underset{\mathrm{coupled}}{\underbrace{{\stackrel{\mathrm{‾}}{h}}_{\mathrm{C}}}},\underset{\mathrm{decoupled}}{\underbrace{{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{D}}+\sum \mathrm{\Delta }{\mathit{\eta }}_T}}=\underset{\mathrm{coupled}}{\underbrace{{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{C}}}}.\end{array}$$

If the dynamics and thermodynamics act independently in a decoupled manner, then the domain-averaged sea ice thickness and viscosity derived from the dynamic model, ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{D}}$ and ${\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{D}}$, combined with the cumulative thermodynamic changes per time step, Δh_T and Δη_T, should be equivalent to the domain-averaged thickness and viscosity of the coupled model, ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{C}}$ and ${\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{C}}$. Any deviation between the left- and right-hand sides of Eq. (26) would indicate nonlinear behaviour in the evolution of these variables, highlighting the significance of the coupling.

The thermodynamic changes per time step, Δh_T and Δη_T, can be obtained by running the thermodynamic model separately. However, we employ the coupled model with the wave amplitude set to zero for simplicity. This ensures that sea ice dynamics are excluded, while preserving the correct category proportions across the domain.

2.4 Model configurations

We design and test a few configurations to showcase the dynamics of the ice floe-grease ice heterogeneous system and to demonstrate the effects of the coupling.

Figure 4 Sea ice distribution derived from a SAR image, illustrating (a) SAR image acquired by the COSMO-SkyMed (CSK) satellite on 22 July 2022, at 07:02 UTC, (b) sea ice concentration from AMSR2 on the same day, (c) selected 1024×1024 pixel region corresponding to 100 % sea ice concentration, and (d) binarised sea ice domain used to initialize the full-field configuration.

To realistically distribute the sea ice types within the computational domain, we derive the full-field configuration from a synthetic aperture radar (SAR) image (Fig. 4a) acquired by the COSMO-SkyMed (CSK) satellite on 22 July 2022, at 07:02 UTC, in support of the SCALE-WIN22 research expedition (Vichi, 2023). The SAR image provides the intensity of the reflected radar signal, which can be assimilated to derive surface properties. The open ocean and ice-covered regions are clearly distinguishable, with the open ocean appearing darker. For our analysis, we select a subregion within the ice-covered area, outlined by the green rectangle in the figure. A comparison with the sea ice concentration derived from AMSR2 satellite on the same day (Fig. 4b) confirms that this subregion corresponds to an area of 100 % sea ice concentration at the 25 km scale. In the region of interest (Fig. 4c), distinct wave patterns are clearly visible that may confound the retrieval of heterogeneity. To remove these wave signatures, a mask is applied in Fourier space. The resulting intensity variations are then interpreted as differences in ice type. Pixels with a filtered amplitude greater than the median are classified as ice floes, while those with lower amplitude are identified as interstitial grease ice. We observed that using a different threshold would change the relative distribution of ice floes and grease ice. Therefore, we tested this by comparing model results across different subregions of the domain, each containing varying proportions. Figure 4d shows the upper-left corner of the selected region following the binarisation process. The final result is a 504×504 grid with a 10 m pixel resolution.

Thickness information cannot be obtained from Fig. 4, therefore, reference visual observations collected during the SCALE-WIN22 research expedition (Vichi, 2023) were used to supplement the analysis.

Figure 5a shows the initial thickness prescribed at t=0 s for the full-field case. The sea ice in the region of the SAR image had a variable thickness ranging from a few centimetres for grease ice to 0.4–0.5 m for pancake ice. Since most sea-ice models do not simulate ice formation starting from frazil ice aggregation (Smedsrud, 2011; Smedsrud and Martin, 2015), a minimum thickness threshold of 0.1 m is typically used (e.g. Rousset et al., 2015). Accordingly, we randomly initialised three ice floe thicknesses just above this threshold, h_f=0.15 m, h_f=0.13 m, and h_f=0.11 m, along with a grease ice thickness slightly below the smallest ice floe thickness, set to h_g=0.08 m. The use of thinner ice is also more compatible with the simplified experimental design in which waves' propagation is not affected by the ice medium. The domain, measuring 5040×5040 m², is discretised using a uniform grid with cell dimensions of 10×10 m².

Figure 5 Initial layout along with additional information for (a) the full-field case, and (b–e) the test cases.

The complementary test cases, shown in Fig. 5b–e, represent an idealised version of Fig. 5a, featuring large ice floes with narrow connections, a characteristic frequently observed in the full-field case. They are designed to clarify the behaviour observed in the full-field case in a more controlled and simplified setting. Two circular floes of different sizes and thicknesses are linked by a narrow connection, with thickness spatially varying from one floe to the other. In this configuration, we investigate the effect of the narrow connection and analyse the domain-averaged viscosity by varying its orientation with respect to the imposed wave (always from the west). The geometry scales were chosen in relation to the wave characteristics (see Table 1). The wiggles were included in Fig. 5e to test whether the irregular shape induces a significant difference in the domain-averaged viscosity results. The circular ice floes have a radii of r=60 m and r=120 m, each with initial thicknesses of h_f=0.11 m and h_f=0.15 m, respectively. The grease ice has the same thickness as in the full-field configuration. The domain dimensions are 720×720 m², discretised using a uniform grid with cell sizes of 2×2 m².

Periodic boundary conditions are applied to all boundaries in the full-field and test cases. Both domains are subjected to a harmonic propagating wave forcing, characterized by a wave amplitude, direction and period, as shown in Table 1. The chosen wave period (or wavelength) ensures that multiple wavelengths fit exactly within the domain length, avoiding potential numerical issues associated with periodic boundary conditions. Wind forcing is not considered, allowing us to isolate the effects of wave-ice interactions.

All simulations in this study are conducted over a 24 h period to allow sufficient time for potential thermodynamic processes to develop and to facilitate a direct comparison between dynamic simulations and coupled dynamic and thermodynamic simulations. It is important to note that the current simulations do not account for the phase transition between grease ice and ice floes, which will be the focus of future works.

3 Results

In this section, we present simulation results to demonstrate the dynamic model and the coupling between the dynamics and thermodynamics models. Three simulations are compared: dynamics only, thermodynamics only, and fully coupled. The analysis focuses on two key variables: sea ice thickness and shear viscosity.

3.1 Full-field case

The propagating waves cause spatially heterogeneous variations in the shear viscosity field of grease ice and ice floes, as shown in the snapshots after 24 h (Fig. 6). These variations reflect both the intrinsic differences in their rheological laws and the local thickness changes induced by wave forcing, which also implicitly modulates viscosity. The influence of wave direction and period on shear viscosity is also visually evident, and the time series of the domain-averaged values presented in Fig. 7 quantify these differences.

Figure 6 Sea-ice shear viscosity (in kg s⁻¹) after 24 h from the full-field configuration for waves with different directions (from the north, east, south or west) and two selected wave periods: (a–d) T=8.8 s, and (e–h) T=15.2 s.

The shorter wave period, T=8.8 s, results in a higher number of regions with lower sea ice shear viscosity ($<\mathrm{1.4}\times {\mathrm{10}}^{\mathrm{8}}$ kg s⁻¹) compared to the longer wave period. This is attributed to higher strain rates, consistent with the sea ice rheology (Eqs. 12 and 13). Lower viscosities are predominantly observed in regions where the ice floes are narrow, while larger ice floe regions remain largely unaffected by variations in wave period. The largest reduction in shear viscosity is observed with changing the wave direction; waves from the north and south are similar to each other but different from the east-west directions. This highlights the role of the orientation of narrow connections between ice floes relative to the wave direction in determining the sea ice viscosity, which is discussed in detail with the test cases in Sect. 3.2.

The time series of the domain-averaged sea ice shear viscosity, shown in Fig. 7, are dominated by the larger ice floe values. The curves show a small negative trend with high-frequency oscillations in the short-wave cases; which are also further analysed in Sect. 3.2. The simulations in the north–south directions show the highest values, with less sensitivity to wave period. In contrast, the different wave periods result in approximately a 20 % difference in the east–west case.

Figure 7 Domain-averaged sea ice shear viscosity (in kg s⁻¹) for three different wave periods, with the wave direction coming from (a) the north, (b) the south, (c) the east, and (d) the west.

The role of heterogeneity, indicated by different percentages of grease ice and ice floes, is analysed in Fig. 8, where we partitioned the domain into 36 subdomains, each identified by a unique number, and calculated spatial and temporal averages of viscosity in relation to the ratio between ice floes and the grease ice. We focused on the lowest wave period because it shows the highest variability, and selected the south and west directions because of the similarity of the results for opposite directions. The time series for every subdomain are shown in Figs. A1 and A2, in which we observe that the shorter wave periods exhibit a highly dynamic response that changes substantially between the subdomains and with respect to the domain-averaged results in Fig. 7.

The colour distributions in the heat maps of the ice floes percentage and viscosities (Fig. 8b–d) reveal the same pattern. The variation in sea ice viscosity is primarily determined by the percentage of ice floes within each subdomain, determining the overall magnitude of the mean viscosity.

Figure 8 (a) Subdomain partitioning and reference numbers, with colours indicating one example combination for each grid size; (b) prescribed percentage of ice floes in each subdomain; (c) the mean sea ice shear viscosity (in kg s⁻¹) in each subdomain over a 24 h period for waves propagating from the south with a wave period of T=8.8 s, (d) same as (c) but for waves propagating from the west; (e) percentage difference in viscosity between wave periods of T=15.2 s and T=8.8 s for waves propagating from the south, and (f) same as (e) but for waves propagating from the west.

Figure 8e and f illustrates the relative difference in shear viscosity between the highest and lowest wave periods. The viscosity values in the case with waves from the west are significantly higher than those from the south, indicating an increased sensitivity to wave periods. As previously mentioned, this discrepancy is attributed to the presence and orientation of the narrow connections and not just to the percentage. This behaviour is further supported by the subdomains with low relative differences, which appear in the rightmost panels of Fig. 8e and the top panels of Fig. 8f. These subdomains are dominated by larger ice floes with fewer narrow connections, making them less responsive to the different wave periods.

The results of these simulations can be further summarized in Fig. 9, where we observe an emergent linear response of the mean viscosity of each subdomain to the percentage of ice floes. The domain is limited to 69 %, as this is the maximum value observed in the SAR image (Fig. 8b), and the intercept at 0 % represents the viscosity of grease ice (≈440 kg s⁻¹). The north-south orientation of the incoming wave describes a linear relationship between sea ice viscosity and the percentage of ice floes in each subdomain (see the equation in Fig. 9). This relationship slightly deviates when all data points are included, due to the increased spread at higher floe percentages, which is clearly dependent on both wave direction and period. This indicates that there is a further relationship between the pattern of the simulated heterogeneous field and the wave direction, which adds to the influence of the floe percentage.

Figure 9 Linear regressions of the mean shear viscosity (in kg s⁻¹) for the subdomains in Fig. 8a against the percentage of floes in each subdomain. The symbols represent different wave conditions. The intercept at 0 % ice floes corresponds to grease-ice viscosity (≈440 kg s⁻¹), as indicated by the regression equations shown in the top-left corner.

3.1.1 Scales of heterogeneity

We have performed a scaling analysis on the full-field results to evaluate the presence of a scaling law, which is often found in sea ice kinematics (e.g. Weiss, 2013). In this case, we are interested in the response of viscosity to the spatial scales of heterogeneity, here represented by the combination of interstitial grease ice and ice floes. Therefore, we considered grid sizes of increasing size from 1×1 to 6×6, with the latter corresponding to the full domain (see Fig. 8a) and calculated the mean shear viscosity and standard deviation in all possible group combinations.

The results are presented in Table 3, showing a strong scale invariance of the mean viscosity from 840×840 m² for a 1×1 grid up to 5040×5040 m² for a 6×6 grid. The 840 m scale is already sufficient to capture the heterogeneity of the ice cover, and variations in ice type patterns do not affect the mechanical response at the larger scales up to approximately 5 km. While variance increases with length scale, the average viscosity is well reconstructed, as the 10 m resolution provides a good estimate of the percentage of heterogeneity.

Table 3 Scaling analysis, where size: tile group size, combinations: possible combinations, N: total number of subdomains, viscosity S: mean sea ice viscosity (south), std S: standard deviation (south), viscosity W: mean sea ice viscosity (west), std W: standard deviation (west).

This behaviour is independent of the wave direction, as also shown in Table 3. The direction changes the absolute value, but we can conclude that the domain-averaged sea ice viscosity is primarily controlled by the ratio of ice floes to grease ice, which is also scale-invariant in this configuration.

Based on the inclusion of smaller scale processes that we assume realistic, the emergence of the linear relationship presented in Fig. 9 and the strong scale invariance of the mean viscosity of sea ice, we are confident that our results can be used to extract properties at larger scales as further discussed in Sect. 4.

3.1.2 Dynamics and thermodynamics coupling

The coupling between the dynamics and thermodynamics models is examined using the full-field case, with the initial layout illustrated in Fig. 5a. We classify both ice floes and grease ice into six distinct thickness categories:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{Ice}\mathrm{floes}:\mathrm{0}\text{–}\mathrm{0.10},\mathrm{0.10}\text{–}\mathrm{0.12},\mathrm{0.12}\text{–}\mathrm{0.14},\mathrm{0.14}\text{–}\mathrm{0.16},\\ {\displaystyle }& {\displaystyle }\mathrm{0.16}\text{–}\mathrm{0.18},\mathrm{0.18}\text{–}\mathrm{0.28}\mathrm{m}.\\ {\displaystyle }& {\displaystyle }\mathrm{Grease}\mathrm{ice}:\mathrm{0}\text{–}\mathrm{0.06},\mathrm{0.06}\text{–}\mathrm{0.09},\mathrm{0.09}\text{–}\mathrm{0.12},\mathrm{0.12}\text{–}\mathrm{0.15},\\ {\displaystyle }& {\displaystyle }\mathrm{0.15}\text{–}\mathrm{0.18},\mathrm{0.18}\text{–}\mathrm{0.24}\mathrm{m}.\end{array}$$

The upper limits of 0.28 m for ice floes and 0.24 m for grease ice are considered sufficient, as significant thickening beyond these thresholds is unlikely within a 24 h period.

Figure 10 presents the spatial distribution of sea ice thickness and shear viscosity at t=24 h for waves propagating from the west with a period of T=8.8 s, comparing results from the dynamics-only model, the thermodynamics-only model, and the fully coupled dynamics and thermodynamics model (see also the videos in the Supplement in the Appendix).

Figure 10 (a–d) Sea ice thickness (in m), and (e–h) sea ice shear viscosity (in kg s⁻¹) for waves originating from the west with a wave period of T=8.8 s, comparing the outcomes of (a, e) the dynamics-only model, (b, f) the thermodynamics-only model, and (c, g) the fully coupled dynamics and thermodynamics model. Panels (d) and (h) show the difference in sea ice thickness and shear viscosity between the fully coupled model and the dynamics-only model.

As previously mentioned in Sect. 2.3, the spatial distribution of the thermodynamics-only model is obtained by using the coupled model with wave amplitude set to zero, thereby excluding sea ice dynamics while preserving the correct category proportions. As a result, Fig. 10b and f exhibit a uniformly increased ice thickness and shear viscosity compared to the dynamics-only results in Fig. 10a and e. This uniformity arises from the absence of ice motion, which prevents interaction and, consequently, the redistribution of sea ice thickness and viscosity between the two ice types. The increase in thickness, as illustrated by the difference between the fully coupled model and the dynamics-only model in Fig. 10d, is most pronounced in the grease ice regions due to the higher growth rate of grease ice compared to ice floes. However, this increase is not visible in the viscosity distributions due to the chosen colour bar, which emphasizes the ice floes that dominate the overall viscosity within the domain.

When comparing the shear viscosity fields from the dynamics-only model in Fig. 10e, and the fully coupled dynamics and thermodynamics model in Fig. 10g, we observe that the spatial distribution of lower-viscosity regions, highlighted in green, remains nearly identical. However, the difference plot in Fig. 10h shows that, in some of these regions, the dynamics-only model exhibits higher viscosities (indicated by negative values), despite the increased thickness in the fully coupled model. This highlights the importance of dynamics in shaping the viscosity distribution of sea ice, as the interaction between dynamics and thermodynamics is inherently nonlinear.

Figure 11a and b illustrates the comparison of the domain-averaged sea ice thickness and shear viscosity time series over a 24 h period between the dynamics-only model and the fully coupled dynamics and thermodynamics model for waves propagating from the west with a period of T=8.8 s.

Figure 11 Comparison between the dynamics-only model and the fully coupled dynamics and thermodynamics model, showing the domain-averaged of (a) sea ice thickness (in m), and (b) sea ice shear viscosity (in kg s⁻¹) for waves originating from the west with a wave period of T=8.8 s. Panel (c) displays the percentage difference between the curves in panel (b), as well as for wave periods, T=12.4 s and T=15.2 s, which are not shown in panel (b).

Initially, both models produce identical results for the domain-averaged sea ice thickness and shear viscosity, reflecting the spin-up time required for the dynamic model to reach equilibrium conditions. The spin-up time, equivalent to five times the wave period, marks the phase before the thermodynamic model is incorporated, as illustrated in Fig. 3 in Sect. 2.3.

In the dynamics-only model, sea ice thickness remains constant over time, as mass conservation is inherently preserved. In contrast, the dynamics and thermodynamics model shows an increase in sea ice thickness due to thermodynamic ice growth. Over a 24 h period, the domain-averaged sea ice thickness, accounting for both ice floes and grease ice, increases by over 1 cm. The domain-averaged sea ice thickness results are nearly identical for the two longer wave periods. This similarity arises from the identical heat fluxes applied across all scenarios.

A similar trend is observed in the domain-averaged shear viscosity, as depicted in Fig. 11b, with viscosity increasing due to thermodynamic sea ice growth. Despite the difference in viscosity values, both viscosity curves exhibit a comparable shape, likely attributed to the similar spatial distribution of lower-viscosity regions in both cases, as discussed previously. Figure 11c illustrates the relative difference, expressed as a percentage, between the domain-averaged shear viscosity for the dynamics-only model and the dynamics and thermodynamics model across three wave periods, T=8.8 s, T=12.4 s and T=15.2 s. The relative difference increases slightly with the wave period.

The relationship between the dynamics and thermodynamics models is subsequently examined using the method described in Sect. 2.3 and Eq. (26), with the same configuration as above. Figure 12a and b presents a comparison of the domain-averaged sea ice thickness and shear viscosity, contrasting the decoupled dynamics and thermodynamics models with the fully coupled dynamics and thermodynamics model.

Figure 12 Comparison between the decoupled dynamics and thermodynamics models and the fully coupled dynamics and thermodynamics model, showing the domain-averaged of (a) sea ice thickness (in m), and (b) sea ice shear viscosity (in kg s⁻¹) for waves originating from the west with a wave period of T=8.8 s. Panel (c) displays the percentage difference between the curves in panel (b), as well as for wave periods, T=12.4 s and T=15.2 s, which are not shown in panel (b).

A negligible discrepancy is observed in the domain-averaged sea ice thickness, as depicted in Fig. 12a. In contrast, the difference in domain-averaged sea ice shear viscosity between the decoupled and coupled models (Fig. 12b) grows over time, with a relative difference of 1.3 % at t=24 h (Fig. 12c). This suggests that the decoupled configuration overestimates the domain-averaged sea ice viscosity. Notably, the relative difference diminishes with increasing wave period, decreasing from 1.3 % to 0.1 % at t=24 h. Over longer timescales, on the order of five to ten days, these differences will accumulate, potentially resulting in a more substantial difference between the decoupled and coupled models.

3.2 Test cases

In this section, we present a detailed analysis of the idealised configuration introduced in Sect. 2.4 describing two larger floes of different thicknesses connected by a narrow bridge. This configuration is designed to help interpret and better understand the results obtained in the full-field case. This case also shows the response of the grease ice component to the wave that was not visually evident in Sect. 3.1 due to the larger scale. The focus will be on the role of floes' orientation and connectivity in relation to the wave direction.

The model is able to describe the features of grease ice thickness and its interaction with the ice floes under different wave conditions, as indicated by the darker shading – which changes position depending on the orientation and wave period – near the interface between the grease ice and ice floes (see Fig. 13a–h and the supplementary videos in the Appendix).

Figure 13 (a–d) Sea ice thickness (in m; the colour bar ranges are chosen to emphasize the wave motion in grease ice), and (i–l) sea ice shear viscosity (in kg s⁻¹) for waves with a wave period of T=8.8 s. (e–h) Sea ice thickness, and (m–p) sea ice shear viscosity for waves with a wave period of T=15.2 s. The black box in the top-right corner highlights the narrow connection between two ice floes.

Figure 13i–p focuses exclusively on the shear viscosity of the two ice floes and the narrow connection between them. The viscosity of the two ice floes is unaffected by the wave period, as both sets of figures display similar colour patterns within the ice floes. In contrast, the viscosity of the narrow connection depends on both the wave period and its orientation relative to the wave. When the connection is aligned perpendicular to the wave front, the viscosity is minimally affected, whereas it reaches its minimum when the connection is parallel to the wave front.

Figure 14 Domain-averaged sea ice shear viscosity (in kg s⁻¹) for three different wave periods, with the narrow connection in (a) horizontal, (b) vertical, (c) diagonal, and (d) zigzag direction.

If we now calculate the spatial averages as done for the full-field case in Sect. 3.1, we observe a similar response on the mean shear viscosity. Figure 14 shows substantial differences between the horizontal and vertical orientations in Fig. 14a and b, which correspond to the bridges aligned perpendicular and parallel to the wave front. We notice that the absolute values and the relative changes in viscosity are scaled to the size of this test-case configuration, and to the shape and proportion of ice floes in the domain. This configuration has a larger proportion of grease ice (see Fig. 5), which was chosen as a compromise to illustrate both the response of the grease ice and the role of ice floes' orientation. In the orientation perpendicular to the wave front (Fig. 14a), the domain-averaged sea ice viscosity is unaffected by the different wave periods because the difference between the curves is negligibly small. This can be attributed to the orientation, as both floes (including the bridge) behave more as a rigid body, making the system less sensitive to the wave period. In this configuration, the bridge primarily experiences compression and tension. In contrast, when the bridge is parallel to the wave direction (Fig. 14b), the floes can move more independently, creating a higher strain rate in the bridge and, consequently, a lower viscosity. As a result, we observe a greater separation between the curves with the largest wave period corresponding to the highest viscosity. This explains why different responses to the wave periods are observed across subdomains, as illustrated in Figs. A1 and A2. These test-case results also demonstrate that the shape and orientation of the bridge influence the oscillations in the viscosity curves shown in Fig. 7 and in Figs. A1 and A2.

Figure 14c and d presents the results for the inclined configuration. The straight inclined bridge is more similar to the results shown in Fig. 14a, while Fig. 14d exhibits more pronounced oscillations that are also observed in some subdomains of the full-field configuration (Figs. A1 and A2). These oscillations, which are of negligible magnitude for this test-case configuration, are more pronounced in the full-field subdomains. We attribute this to the shape of the zigzag connection, which adds complexity to the shear dynamics. While noise from numerical interface approximations cannot be excluded, the possible emergence of internal resonance – due to the domain periodicity and wavelength selection – as well as harmonics amplified in the full-field case also requires further investigation.

4 Discussion

The proposed model, WIce-FOAM 1.0, is designed to explore the dynamic response of the heterogeneous sea-ice cover, composed of consolidated ice floes and interstitial grease ice, to harmonic wave forcing. The rationale for this study originates from the need to describe the characteristics of the Antarctic marginal ice zone, which is a mosaic of pancake ice floes cemented together (e.g. Wadhams et al., 2018; Alberello et al., 2019; Vichi, 2022). Our results highlighted that the heterogeneity of the sea-ice surface, and not only the thickness, alters the overall mechanical response of sea ice to wave stresses. The spatially averaged viscosity field is sensitive to the wave direction and the period, and the intensity of the response depends on the percentage of solid floes. We found a lower shear viscosity in the regions where the floes aggregate in narrower formations. Through selected test cases, we found that the presence and orientation of the narrow junctions between ice floes relative to the wave direction play a key role in the mean viscosity of an area of heterogeneous sea ice. We also observed changes in the magnitude and oscillatory behaviour of the mean shear viscosity that are indicative of resonance at smaller scales that can be propagated to the kilometre scales.

Our approach is agnostic, and we deliberately did not consider in our simulations the likely existence of an underlying relationship between the patterns of heterogeneity and the direction of the waves. It is important to remember that these patterns are not sea-ice bands (e.g. Saiki et al., 2021), since this is not an open drift condition and the scales are smaller than the 10 km usually observed for the bands. Nevertheless, we observe a clear anisotropy in the response of the shear viscosity to the wave direction. We removed the wave patterns in Fig. 4, which indicated a propagation from west-northwest, and we do observe a series of bands oriented perpendicularly to this direction (see Fig. 5a). The bands may have been arranged this way by the wave motion, and we notice that the narrow bands show a lower viscosity with waves coming from the west (Fig. 6d and h).

The addition of thermodynamic processes alters the viscosity of sea ice in the dynamic model by 3 % at t=24 h. This small percentage increases slightly with an increase in the wave period. However, the coupling between dynamics and thermodynamics exhibits a nonlinear response that grows over time. Despite the detailed formulation of the fully coupled model, the dynamics requires a small time step of less than a second, resulting in significant computational costs. While differences observed over a 24 h period remain limited, they accumulate over time. This indicates that as computational capabilities improve, enabling extended simulations over a longer period, the impact of these differences will become increasingly apparent. However, given that environmental and waves' characteristics may change with the storm scales of a few days, this nonlinearity may not grow indefinitely. In addition, the lack of phase change between grease ice and solid floes limits the length of our simulations, since grease ice grows faster than the solid floes and would become thicker. Some degree of compaction or interaction with the rims of the ice floes is expected to occur, but there are no direct observations of this process.

One major outcome of our approach is an emerging linear relationship between the mean shear viscosity and the ice floe percentage within each subdomain (Fig. 9). Several test cases considering different shapes and sizes of the solid floes were implemented during the onset of this work, and led to responses of the mean viscosity field that were of difficult interpretation. This is demonstrated by the figures in the Appendix, in which different configurations lead to different responses of the mean viscosity to the waves characteristics. Only the use of the pseudo-realistic full-field configuration from an SAR image allowed us to discriminate emerging patterns. There is not enough knowledge to accurately determine the type of ice from radar intensity, and we used a threshold derived from the field itself to arbitrarily assign rheology and thickness to the ice field. We do not imply generalization to all SAR images, but we suggest that this model can already be used to determine the scales of heterogeneity and inform the design of parameterisations that include the effects of waves on sea-ice mechanics. Without explicitly resolving for the waves, it should be possible to derive a parameterisation of the heterogeneous rheology based on the percentage of young ice. Current sea ice models do not yet implement a deterministic calculation of young ice, but a few parameterisations are already included to resolve polynya conditions (Fang et al., 2024). Our scaling analysis indicates that the percentage of ice floes can be obtained from SAR images at the scales of current sea ice models, and used to derive a mean viscosity value that accounts for the wave action based on this emergent linear relationship. The wave direction and the period partly affect the magnitude, and we notice in Fig. 9 that there is a likely range between 30 %–70 % over which this happens. Our range of values is limited to 69 %, but we expect that at 100 % floe coverage, the viscosity would converge to the maximum extrapolated value of 2.6×10⁸ kg s⁻¹ from Eq. (13).

We acknowledge that this model contains major assumptions and compromises due to the limitations of our current knowledge on these processes in the Antarctic and the chosen computational framework. The model does not include leads because it focuses on regions where the cover is 100 % but still thin, and thus is affected by the penetration of waves. A multiphase approach involving more than two phases cannot be implemented in OpenFOAM without abandoning the Volume of Fluid (VOF) technique currently used, the IsoAdvector method. This would degrade the representation of the interface between interstitial grease ice and consolidated ice floes. For this preliminary study, which builds on the work of Marquart et al. (2021, 2023), we preferred to preserve the interfaces of the mosaic, but we do not exclude a further development to include a seawater phase and the phase transition between water and the different types of sea ice. Our results have highlighted the importance of the floes' percentage in determining the mechanical response to the waves, and therefore other more diffusive methods that ensure mass conservation at the expense of details can be considered in the future.

Another relevant assumption is that the rheology of both constituents is known, with viscous grease ice and viscous-plastic solid floes (see Sect. 2.1.2). This latter method is the standard parameterisation proposed by Hibler (1979) used in most sea-ice models, which do not account for the complex stratigraphy and the prevalence of granular structure observed in Antarctic floes (e.g. Lytle and Ackley, 2001; Skatulla et al., 2022; Tison et al., 2020; Johnson et al., 2023). There are no recent works on the compressive or shear response of granular ice (Paul et al., 2023), and the literature, as well as the numerical models, assume that sea ice is columnar in its structure. Finally, we only considered waves to be harmonic and unidirectional, and ignored the attenuation. This is justified by the duration of our simulations in terms of time; for this reason, we considered experiments with thin ice, which also allowed thermodynamic effects to become apparent within a one-day time scale.

5 Conclusions

We presented a modelling framework, WIce-FOAM 1.0, developed in OpenFOAM to explicitly resolve the dynamic and thermodynamic processes in sea ice composed of two distinct ice types at the sub-kilometre scale. In our model, cells identified as ice floes or grease ice may contain both ice types, though one type dominates. This model serves as an initial testing platform to explore the response of heterogeneous thin ice to the effects of waves. It is a useful tool for designing experimental in situ and laboratory setups and for deriving emerging relationships to inform the parameterisation of larger-scale numerical sea ice models. We demonstrated its functionality through both realistic and idealised simulations based on Antarctic sea ice examples, where the floes consist of agglomerated pancake ice and the interstitial component is grease ice. However, the model can be applied to any sea ice configuration where solid floes and interstitial ice coexist. Despite being constrained by limited computational time scales, our initial results show an emerging linear relationship between the fraction of solid floes across multiple spatial scales and the mean shear viscosity. The mean viscosity is affected by the direction and period of the incident waves, which, in principle, would require sea ice models to explicitly resolve wave features. However, our findings suggest the possibility to parameterise the mechanical response of heterogeneous sea ice to waves, even without their explicit representation, assuming the model can incorporate wave characteristics and simulate the presence of interstitial ice. Since our work assumes unattenuated waves, further research is required to include the feedback between the penetrating waves and sea ice, as well as the phase transition between interstitial grease ice and solid floes.

Appendix A

Figure A1 Spatially-averaged sea ice shear viscosity (in kg s⁻¹) over a 24 h period across 36 subdomains, for three different wave periods, with wave propagation from the south. The black, blue and red curves correspond to wave periods of T=8.8 s, T=12.4 s, and T=15.2 s, respectively. In most panels, shear viscosity is the highest for the longer wave periods and the lowest for the shorter. The shortest wave period exhibits a more dynamic response, which changes substantially between the subdomains and with respect to the domain-averaged results in Fig. 7. Panel number 3 consists solely of grease ice, resulting in shear viscosity values five orders of magnitude smaller than those of the adjacent panels. A clear trend is observed between grease ice viscosity and wave period, with the shear viscosity increasing as the wave period increases. Moreover, we observe that the spatially-averaged shear viscosity of grease ice is independent of the wave direction, as the results from panel 3 are the same as in Fig. A2.

Figure A2 As in Fig. A1 but for wave propagation from the west.