Modelling sea ice in the marginal ice zone as a dense granular flow with rheology inferred from discrete element model data

Abstract

The marginal ice zone represents the periphery of the sea ice cover. In this region, the macroscale behaviour of the sea ice results from collisions and enduring contact between ice floes. This configuration closely resembles that of dense granular flows, which have been modelled successfully with the $\mu (I)$ rheology. Here, we present a continuum model based on the $\mu (I)$ rheology that treats sea ice as a compressible fluid, with the local sea ice concentration given by a dilatancy function $\varPhi (I)$. We infer expressions for $\mu (I)$ and $\varPhi (I)$ by nonlinear regression using data produced with a discrete element method (DEM) that considers polygon-shaped ice floes. We do this by driving the sea ice with a one-dimensional shearing ocean current. The resulting continuum model is a nonlinear system of equations with the sea ice velocity, local concentration and pressure as unknowns. The rheology is given by the sum of a plastic term and a viscous term. In the context of a periodic patch of ocean, which is effectively a one-dimensional problem, and under steady conditions, we prove this system to be well-posed, present a numerical algorithm for solving it, and compare its solutions to those of the DEM. These comparisons demonstrate the continuum model's ability to capture most of the DEM results accurately. The continuum model is particularly accurate for ocean currents faster than 0.25 m s$^{-1}$; however, for low concentrations and slow ocean currents, the continuum model is less effective in capturing the DEM results. In the latter case, the lack of accuracy of the continuum model is found to be accompanied by the breakdown of a balance between the average shear stress and the integrated ocean drag extracted from the DEM. Since this balance is expected to hold independently of our choice of rheology, this finding indicates that continuum models might not be able to describe sea ice dynamics for low concentrations and slow ocean currents.

1. Introduction

The periphery of the ice cover is known as the marginal ice zone (MIZ) and consists of relatively small, typically polygon-shaped ice floes. It is often defined as the region where ocean waves play an important role in shaping the morphological properties of the ice (Dumont 2022). On large scales, sea ice dynamics is typically described with Hibler's model (Hibler 1979), which treats ice as a viscoplastic fluid whose yield strength depends on the sea ice concentration and thickness. This model was developed for the central ice pack, where ice floes are closely interlocked, and deformation is due mostly to the opening of leads or the formation of ridges. Recently, elasto-brittle rheologies have also been used to model the evolution of the central ice pack. This class of rheological models, which was first proposed by Girard et al. (2011), appears to be superior to Hibler's model in capturing the ice deformation field (Rampal et al. 2019).

In the MIZ, however, it is the collisions and enduring contact between ice floes that give rise to the macroscale dynamical properties of the ice cover (Feltham 2005; Herman 2011, 2022). This configuration closely resembles that of dense granular flows, albeit at different spatial scales, since practically all studies for granular materials consider e.g. polystyrene beads, glass beads and sand, whose particles’ diameters are of the order of 0.1 and 1 mm (GDR MiDi 2004). Dense granular flows have been modelled successfully with the so-called $\mu (I)$ rheology (Da Cruz et al. 2005). The dense granular flow regime is understood as a transition between the quasi-static and dilute flow regimes. Whenever grain inertia is negligible, a quasi-static regime emerges that is often modelled as an elastoplastic solid (Nedderman 1992). The critical state at which plastic deformation occurs is characterized with a Coulomb-like criterion dependent on a so-called internal angle of friction (Wood 1990). Conversely, under great agitation and/or dilute concentrations of grains, particles interact only through binary, instantaneous, uncorrelated collisions. As a result, ideas from kinetic theory become applicable in this dilute regime (Jenkins & Savage 1983). However, in dense granular flows, grains interact via collisions and enduring contacts, such that inertial effects are important yet the collisions may no longer be assumed to be binary, instantaneous, or uncorrelated in general. This transitional regime is characterized in terms of the inertial number $I$ and an effective friction coefficient $\mu$ that is dependent on $I$ (Da Cruz et al. 2005).

Existing models for the MIZ recognize the importance of both collisions and plastic deformation, and derive rheological models based on first principles (Shen, Hibler & Leppäranta 1987; Gutfraind & Savage 1997; Feltham 2005). Recently, Herman (2022) suggests the use of the $\mu (I)$ rheology for modelling sea ice in the MIZ, and derives a $\mu (I)$ function from computations performed with a discrete element method (DEM). In these computations, disk-shaped ice floes are sheared by a moving wall in the classical manner of rheological studies. Unlike the previous models for the MIZ, Herman (2022) infers the rheological properties from data generated by a DEM. In particular, Herman (2022) fits a $\mu (I)$ function to the DEM data, although the resulting continuum model and its accuracy in replicating the DEM results is not examined.

This work represents an advance in the development of a continuum model for the MIZ that could improve the accuracy of Hibler's model, which is currently used in large-scale climate models over the MIZ; see e.g. Danabasoglu et al. (2020). A comparison of Hibler's model with the DEM data is presented in § 6.1, where it is demonstrated that it cannot capture the DEM results accurately. We extend the investigation initiated in Herman (2022), and explore the $\mu (I)$ rheology's accuracy in modelling sea ice dynamics in the MIZ. We infer a $\mu (I)$ function from data produced with the DEM implemented in SubZero (Manucharyan & Montemuro 2022). This DEM considers polygon-shaped ice floes that are driven by oceanic currents in an open patch of ocean, a set-up that we believe to be more natural for studying sea ice than the classical shearing test with a moving wall and disk-shaped ice floes. This inference results in a continuum viscous fluid model whose rheology is given by the sum of a viscous term and a plastic term. Moreover, for this system to be well-posed, the emerging model problem requires the continuum to be compressible and complemented with a constraint on the global sea ice concentration. Assuming the continuum to be compressible requires the inference of a dilatancy function $\varPhi (I)$ from the DEM computations that establishes a relationship between local sea ice concentration and the inertial number $I$.

The contributions of this paper can be summarized as follows. (1) Inference of the $\mu (I)$ and $\varPhi (I)$ constitutive functions for sea ice in the MIZ from data produced with the DEM. These computations are performed in an open ocean configuration where the sea ice is sheared by ocean currents. (2) Analysis of the resulting continuum model, establishing the existence and uniqueness of solutions. (3) Determination of the continuum model's range of validity by comparing its numerical solutions to those of the DEM.

We remark that the analysis of the continuum model and its comparisons with the DEM are restricted to a steady one-dimensional set-up. The model can be extended to unsteady two-dimensional problems as explained in § 2.1, although we expect that new complications will arise with these extensions. For example, Barker et al. (2015) demonstrate the emergence of time-dependent instabilities in $\mu (I)$ models, which Schaeffer et al. (2019) remedy with further modifications of the model. These potential complications should be studied carefully in future investigations.

This paper is structured as follows. In § 2, we formulate the continuum model, first in a general two-dimensional unsteady setting, then in the one-dimensional steady configuration considered in this paper. In this formulation, two functions, $\mu (I)$ and $\varPhi (I)$, are to be inferred from DEM data. This inference is presented in § 3. Section 4 contains a detailed analysis of the continuum model resulting from this inference. This analysis examines several properties of the momentum equation, the numerical solution of the continuum model, and its well-posedness. Then in § 5, we compare the continuum model and the DEM. This comparison allows us to establish the range of validity of the continuum model and its limitations. In § 6, we discuss the similarities and differences between our continuum model and other sea ice models, such as Hibler's model. We then end this paper with § 7, where we recommend potential extensions of this work to be explored in the future.

2. Mathematical formulation of the continuum problem

The dense flow regime represents a transition between the quasi-static and dilute flow regimes (Da Cruz et al. 2005). This transitional regime is characterized in terms of the inertial number $I$, an effective friction coefficient $\mu (I)$, and, whenever the continuum is assumed to be compressible, a dilatancy function $\varPhi (I)$. Below, we define these three terms and present a general formulation of the $\mu (I)$ rheology in two dimensions and in the one-dimensional steady configuration considered in the subsequent sections of this paper.

2.1. The two-dimensional setting

Although the problems presented in this paper are effectively one-dimensional, we first present the general form of the flow model in two dimensions for completeness. We denote the ice velocity, concentration and Cauchy stress tensor by $\boldsymbol {u}$, $A$ and $\boldsymbol {\sigma }$, respectively. We write the components of the Cauchy stress tensor and the velocity vector field as

(2.1a,b)\begin{equation} \boldsymbol{\sigma} = \left[\begin{array}{@{}cc@{}} \sigma_{xx} & \sigma_{xy} \\ \sigma_{xy} & \sigma_{yy} \end{array}\right] \quad \text{and} \quad \boldsymbol{u} = (u,v), \end{equation}

respectively. We assume that the morphology of the ice floes remains invariant by neglecting all thermomechanical effects, such as fracturing, melting or ridging, that can change the shape of a floe. For simplicity, we also neglect the Coriolis force, ocean tilting and the atmospheric drag (we assume low-wind conditions). Under these conditions, conservation of momentum and mass lead to the following system of equations:

(2.2a)$$\begin{gather} \rho H\,\frac{\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma} + \boldsymbol{t}_o, \end{gather}$$

(2.2b)$$\begin{gather}\frac{\mathrm{D} A}{\mathrm{D} t} ={-} A\,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}; \end{gather}$$

see Hibler (1979) or Gray & Morland (1994). For any scalar or vector-valued function $f$, the material derivative is given by $\mathrm {D} f/\mathrm {D} t = \partial f/\partial t + (\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla })f$. Here, we assume the ice thickness $H$ to be spatially uniform for simplicity, although in general we require an additional equation, analogous to (2.2b) but in terms of $H$, for mass to be conserved. Equation (2.2a) is a depth-averaged statement of conservation of momentum of the sea ice layer; here, $\boldsymbol {t}_o$ is the drag force exerted by the ocean on the sea ice. Given the surface ocean velocity field $\boldsymbol {u}_o$, this drag force is generally parametrized in terms of the drag coefficient $C_o$ and the ocean water density $\rho _o$ by

(2.3)\begin{equation} \boldsymbol{t}_o := \rho_o C_o\,\|\boldsymbol{u}_o - \boldsymbol{u}\|\,(\boldsymbol{u}_o - \boldsymbol{u}), \end{equation}

with $\|\cdot \|$ denoting the Euclidean norm of a vector.

The conservation laws (2.2) must be accompanied by constitutive relations. To write these, we first decompose the Cauchy stress tensor into a pressure term $p$ and its deviatoric component $\boldsymbol {\tau }$,

(2.4)\begin{equation} \boldsymbol{\sigma} = \boldsymbol{\tau} - p\boldsymbol{{I}}, \end{equation}

where $\boldsymbol {{I}}$ is the identity tensor, and define the strain rate tensor $\boldsymbol {{D}}$ and its deviatoric component $\boldsymbol {{S}}$ as

(2.5a,b)\begin{equation} \boldsymbol{{D}} := \tfrac{1}{2}(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{\rm T}) \quad \text{and}\quad \boldsymbol{{S}} := \boldsymbol{{D}} - \tfrac{1}{2}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}) \boldsymbol{{I}}. \end{equation}

In the following, for a given tensor $\boldsymbol {{T}}$, its second invariant is denoted by

(2.6)\begin{equation} \|\boldsymbol{{T}}\|=\sqrt{\tfrac{1}{2}\operatorname{tr}(\boldsymbol{{T}}^2)}. \end{equation}

The fundamental idea behind the $\mu (I)$ rheology is that the constitutive relation for a granular flow depends on the inertial number, a dimensionless quantity defined as

(2.7)\begin{equation} I := \bar{d}\,\sqrt{\frac{H\rho_i}{p}}\,\|\boldsymbol{{S}}\|, \end{equation}

where $\rho _i$ is the ice density, $\bar {d}$ denotes an average ice floe size, and $p$ is the pressure emerging in (2.4) (Da Cruz et al. 2005; Jop, Forterre & Pouliquen 2006; Pouliquen et al. 2006). Throughout this paper, we set $\bar {d}$ to be spatially constant over the whole domain, avoiding the need to consider the transport of this quantity. Savage (1984) interprets the quantity $I^2$ as the ratio between collisional (i.e. inertial) stresses and the total shear stress exerted on the material. Accordingly, for low values of $I$, the inertial effects of grains become negligible, and the flow approaches a quasi-static regime; conversely, as $I$ increases, collisional forces become increasingly important relative to the external forces exerted on the material. The functional relationship that establishes the material's rheology is written in terms of an effective friction $\mu (I)$, defined as

(2.8)\begin{equation} \mu(I) := \frac{\|\boldsymbol{\tau}\|}{p}. \end{equation}

We remark that the effective friction $\mu$ is defined in analogy with Coulomb's model of friction as the ratio between the shear (tangential) stress and the pressure (normal stress). Moreover, it is also helpful to think of the pressure $p$ as a quantification of the material's strength and its resistance to viscous and plastic deformation, as is made clear in § 4. It should be noted that the $\mu (I)$ model is a phenomenological model that has been found to work well with granular media, yet it is unclear if it represents some kind of limit for a large particle system.

To obtain a relationship between stress and strain, we need an additional constitutive law. Jop et al. (2006) propose the following equality that aligns $\boldsymbol {{S}}$ with $\boldsymbol {\tau }$:

(2.9)\begin{equation} \frac{\boldsymbol{{S}}}{\|\boldsymbol{{S}}\|} = \frac{\boldsymbol{\tau}}{\|\boldsymbol{\tau}\|}. \end{equation}

Combining (2.8) and (2.9), the relationship between deviatoric components of the stress tensor and the shear strain can be written as

(2.10)\begin{equation} \boldsymbol{\tau} = \mu(I)\,p\,\frac{\boldsymbol{{S}}}{\|\boldsymbol{{S}}\|}. \end{equation}

Compressible granular flows require a dilatancy law that relates the concentration $A$ to the inertial number $I$:

(2.11)\begin{equation} A = \varPhi(I); \end{equation}

see Da Cruz et al. (2005). In general, $\varPhi$ is found to be a strictly decreasing function of $I$, in such a way that the concentration $A$ decreases with the rate of shearing $\|\boldsymbol {{S}}\|$, a phenomenon know as dilatancy. Moreover, if $\varPhi$ is strictly decreasing, then it is invertible, and one can write an expression for the pressure $p$ analogous to an equation of state in thermodynamics:

(2.12)\begin{equation} p = \rho_i \bar{d}^2H\left(\frac{\|\boldsymbol{{S}}\|}{\varPhi^{{-}1}(A)}\right)^2, \end{equation}

where we have combined (2.7) and (2.11). In the problems considered in this paper, we find the spatial variations in sea ice concentration to be small. Although this would make the assumption of incompressibility reasonable, the periodic one-dimensional nature of these problems renders the dilatancy law (2.11) necessary for the model to be well-posed. This point is explained in § 2.2.

2.2. The steady one-dimensional periodic ocean problem

The model problem considered in this paper consists of a square patch of ocean of length $L$ with periodic boundaries in both the $x$- and $y$-directions. The ice floes floating on this patch are driven by an ocean velocity field that varies only in the $y$-direction, as depicted in figure 1. We neglect time-dependent effects, and consider only steady conditions in the forcing, i.e. $\boldsymbol {u}_o(x,y, t) = (u_o(y),0)$.

Figure 1. The periodic ocean set-up. The domain is a square patch of ocean, periodic in the horizontal and vertical directions. The ice floes are driven by the ocean velocity field (2.15) (in blue), which does not vary in the $x$-direction.

This configuration renders the continuum problem one-dimensional and independent of time, such that $\boldsymbol {u}(x,y, t) = (u(y), 0)$, $A(x,y, t) = A(y)$ and $p(x,y,t) = p(y)$. In this setting, the equations for conservation of momentum (2.2a), together with the constitutive equation (2.10), simplify to the following system on $(0,L)$:

(2.13a)$$\begin{gather} - \frac{\mathrm{d}\sigma_{xy}}{\mathrm{d} y}= \rho_o C_o\,|u_o - u|\,(u_o - u), \end{gather}$$

(2.13b)$$\begin{gather}\frac{\mathrm{d} p}{\mathrm{d} y} = 0, \end{gather}$$

(2.13c)$$\begin{gather}\sigma_{xy} = \mu(I)\,p\,\frac{\mathrm{d} u/\mathrm{d} y}{|\mathrm{d} u/\mathrm{d} y|}, \end{gather}$$

(2.13d)$$\begin{gather}I = \bar{d}\,\sqrt{\frac{H\rho_i}{p}}\,|\mathrm{d} u/\mathrm{d} y|. \end{gather}$$

Due to (2.13b), which represents the balance of momentum in the $y$-direction, the pressure is constant (but unknown) over the domain. Da Cruz et al. (2005) and Herman (2022) find $p$ by enforcing normal stress boundary conditions along a boundary of the domain, but we cannot do the same because the domain is periodic. In the DEM computations, which we introduce in § 3, we set a global ice concentration $A_0\in [0,1]$ that equals the domain-averaged value of the local concentration $A$, such that

(2.14)\begin{equation} \frac{1}{L}\int_0^L A\,\mathrm{d} y = A_0. \end{equation}

Therefore, if we assume the sea ice to behave like a compressible fluid, then condition (2.14) and the dilatancy law (2.11) close the system of equations. In §§ 3 and 5, we justify this modelling choice by demonstrating that dilatancy emerges in the DEM computations and that our model is capable of capturing it accurately.

In this paper, we consider only the following ocean velocity profile for simplicity:

(2.15)\begin{equation} u_o(y) = u_{o,max}(1 - |1 - 2y/L|), \end{equation}

for some maximum velocity $u_{o,max} > 0$. Figure 1 contains a plot of this velocity profile.

3. Inferring the constitutive equations of the system from the DEM

The system of equations presented in § 2.2 is incomplete because we need additional expressions for $\mu (I)$ and $\varPhi (I)$. In this paper, we infer these additional equations from data generated with SubZero, a DEM developed by Manucharyan & Montemuro (2022) and used for modelling sea ice dynamics with polygon-shaped ice floes. Following the set-up presented in § 2.2, we perform runs with $n = 2000$ floes over a square patch of ocean of length $L = 100\,{\rm km}$, driven by the ocean velocity field (2.15). This means that the $\mu (I)$ and $\varPhi (I)$ functions are inferred from DEM solutions with a constant number of floes (and therefore constant average floe size). In theory, the effects of floe size are included in the non-dimensional parameter $I$. In practice, the effects of $n$ are subtle and not too well captured with our continuum model; see figure 12 in § 5 below.

For a given number of floes and a global sea ice concentration $A_0$, the initial configuration of ice floes is generated with SubZero's packing algorithm, which is based on a Voronoi tessellation of the domain (Manucharyan & Montemuro 2022) (see figure 2(a) for an example of the outcome of this packing algorithm). Defining the floe size $d$ as the square root of the floe area, this generates a polydisperse floe size distribution whose histogram we can see in figure 3 for three values of $A_0$; we find that $d/L$ is approximately between 0.01 and 0.04. Although we do not study the effects of polydispersity on the rheology here, we remark that Herman (2022) finds the dilatancy law $\varPhi$ to vary visibly with the degree of polydispersity, while the effective friction $\mu$ presents much smaller variations.

Figure 2. Emergence of dilatancy in the DEM. Here, $A_0 = 0.8$ and $u_o = 1\,{\rm m}\,{\rm s}^{-1}$. Ten runs to steady state are performed from randomized initial conditions. (a,b) Examples of initial and final states, respectively. (c) Ocean and sea ice velocities. (d) Initial (red) and final (black) mean local sea ice concentration, with the standard deviation (bars). Clearly, $A$ decreases in the areas of maximum shear strain; compare (c) and (d).

Figure 3. Floe size distributions for different global sea ice concentrations $A_0$ and a total number of ice floes $n = 2000$. The floe size $d$ is defined as the square root of the floe area.

We run the DEM simulations for $2\times 10^4$ time steps, each of 5 s (the total running time is approximately 27.8 h). In general, we find that the velocity, stress and sea ice concentration, averaged over the last 25 % of the time steps, remain relatively unchanged when a longer computation is performed, hence we consider that a steady state has been reached. This temporal averaging is performed over data that at each time step has been averaged spatially over a grid as described in § A.2. For these simulations, the material parameters that determine the effects of the ocean drag and the collisions among floes are presented in table 1. Collisional forces and the resultant stresses, which determine the fields $\boldsymbol {\sigma }$ and $p$, are computed as explained in § A.1.

Table 1. Values for material parameters used throughout this paper. Here, $C_o$ is the drag coefficient for the ocean current, and $\rho _i$ and $\rho _o$ are the ice and ocean water densities, respectively. The Young's modulus $E$, Poisson's ratio $\nu$, and inter-floe friction coefficient $\mu ^\ast$ are used in the calculation of collisional forces, as described in § A.1.

In figure 4, we plot the values of the friction $\mu = |\sigma _{xy}|/p$ and local concentration $A$ against $I$ for different global concentrations $A_0$ between $0.7$ and $0.95$, and different maximum ocean velocities $u_{o,max}$ between $0.1$ and 1 m s$^{-1}$. The mean ocean velocities in the ocean patch are therefore between $0.05$ and 0.5 m s$^{-1}$, values that are consistent with real observations (Stewart, Klocker & Menemenlis 2019). In all of these computations, we set the ice thickness to $H = 2\,{\rm m}$. We find an increase in the friction $\mu$ and a decrease in the local concentration $A$ as $I$ increases. The decrease in the local concentration of sea ice $A$ with an increase in $I$ is due to dilatancy. Figure 2 presents an example of how dilatancy emerges in the DEM computations: given a random initial distribution of ice floes, when a steady state is reached, the concentration $A$ decreases in the areas where the largest shearing occurs ($y = 1/4$ and $y = 3/4$), and increases elsewhere. Since the global sea ice concentration $A_0$ is constant, in this context dilatancy represents a reorganization of the local concentration profile $A(y)$.

Figure 4. Plots of (a) friction $\mu = |\sigma _{xy}|/p$ and (b) concentration $A$, against $I$. Ten friction and concentration values are extracted from each DEM run by averaging along the grid cells plotted in figure 15. The black lines correspond to the functions (3.1a,b) fitted to the data by minimizing the least squares misfit.

The trends found in the data in figure 4 are well fitted with the following family of functions:

(3.1a,b)\begin{equation} \mu(I) = \mu_0 + \mu_1 I \quad \text{and} \quad \varPhi(I) = 1 - \phi_0 I^\alpha. \end{equation}

A linear behaviour is also found for $\mu$ in Da Cruz et al. (2005). The four parameters $\mu _0,\mu _1,\phi _0,\alpha$ are calculated by minimizing the least squares misfit problem between the points in figure 4 and the functions in (3.1a,b). The resulting values are shown in table 2. For the remainder of the paper, any numerical solution of the one-dimensional system (2.13) is obtained by setting the parameters in functions (3.1a,b) to the values given in table 2.

Table 2. Parameters for the functions $\mu (I)$ and $\varPhi (I)$ in (3.1a,b) obtained by minimizing the least squares misfit with the data plotted in figure 4. These are the numerical values used for computing solutions to the continuum model (2.13) throughout this paper.

Departures from the fitting curves are most visible when the ocean velocities and sea ice concentrations are low; see the case where $u_{o,max} = 0.1\,{\rm m}\,{\rm s}^{-1}$ and $A_0 = 0.7$. Unsurprisingly, in § 5, we also find the greatest misfit between the DEM and the continuum model precisely in this setting, when $u_{o,max} = 0.1\,{\rm m}\,{\rm s}^{-1}$ and $A_0 = 0.7$; see figure 9(m) below. In particular, in this setting, the fundamental balance between shear stress and ocean drag in the DEM is found to no longer hold; see § 5.

The constitutive equation in two dimensions resulting from functions (3.1a,b) is the sum of a plastic term and a viscous term:

(3.2)\begin{equation} \boldsymbol{\tau} = \underbrace{\mu_0 p\,\frac{\boldsymbol{{S}}}{\|\boldsymbol{{S}}\|}}_{plastic} + \underbrace{\mu_1 \bar{d}^2 \sqrt{\rho_i H p}\,\boldsymbol{{S}}}_{viscous}. \end{equation}

A consequence of this linear behaviour is that $\mu$ is approximately constant for small values of $I$, and this is precisely what we see for $I < 10^{-2}$ in figure 4. For $I\ll 1$, we have that $\mu (I) \approx \mu _0$, therefore it is the plastic term that dominates the rheology. This is essentially the quasi-static regime, where collisions are negligible. This plastic term follows from a Mohr–Coulomb yield criterion with internal angle of friction $\tan ^{-1}(\mu _0)$; examples of the Mohr–Coulomb yield criterion used for sea ice modelling can be found in Ip, Hibler & Flato (1991), Gutfraind & Savage (1997) and Ringeisen et al. (2019). The viscous term, which becomes increasingly important as the inertial number $I$ increases, can be associated with the collisional component of the rheology. A viscous rheology is derived in Shen et al. (1987) for modelling the rheological effects of collisions in sea ice, which, as we explain in § 6, is very similar to the viscous component in (3.2).

4. Analysis of the inferred continuum model

This section focuses on the one-dimensional system of equations presented in § 2.2, with the functions $\mu$ and $\varPhi$ taking the form (3.1a,b). First, in § 4.1 we non-dimensionalize the system of equations to understand the relative importance of the different terms involved.

The remaining two subsections explore the existence and uniqueness of solutions to this system, a regularization technique that facilitates its numerical solution, and the behaviour of solutions under different limits. In particular, in § 4.2, we focus only on the momentum equation (2.13b). In § 4.2.1, we show that the momentum equation can be rewritten as a minimization problem. This allows us to establish that solutions to this equation exist and are unique, and it allows us to make sense of the plastic component (see (3.2)) in a rigorous sense (as a variational inequality). Moreover, this equivalence also motivates a regularization of the plastic term that simplifies its numerical solution considerably, as described in § 4.2.2. We end § 4.2 with an analysis of how the sea ice velocity behaves under different limits in parameter values: in § 4.2.3, we explore the behaviour of the sea ice velocity for small and large pressures $p$; then in § 4.2.4, we derive an analytical solution for the sea ice velocity in the purely plastic limit. Understanding these limit solutions for the velocity is helpful for interpreting the DEM results that we present in § 5. It is also useful for the analysis that we present in § 4.3, which considers the complete system of equations (2.13). We begin by presenting a numerical method for solving the complete system in § 4.3.1. Then in § 4.3.2, we sketch out a demonstration of the existence and uniqueness of solutions to the complete system.

We remark that §§ 4.2.1, 4.2.2 and 4.3 are mostly concerned with questions of a mathematical and numerical nature. Although we believe these to be important topics in establishing the suitability of our model for modelling purposes, they are not required for understanding the remainder of the paper. We also note that, for simplicity, the solutions presented throughout this section result from driving the ice floes with the ocean velocity profile (2.15), although the analysis can be extended to more general ocean velocity profiles by following the same steps.

4.1. Non-dimensionalization of the system

For the non-dimensionalization, we set the characteristic magnitudes

(4.1a–d)\begin{equation} [y] = L,\quad [H] = H,\quad [u] = u_{o,max} \quad \text{and} \quad [\sigma] = \rho_i [u]^2 [H] \end{equation}

for the length, thickness, velocity and stress, respectively. We scale the velocities $u$ and $u_o$ with $[u]$, the spatial variables $y$ and $\bar {d}$ with $[y]$, the thickness $H$ with $[H]$, and $\sigma _{xy}$ and $p$ with $[\sigma ]$.

From this point onwards, all variables considered are non-dimensional unless made explicitly clear to the contrary, or units are specified. Keeping the same notation as used for dimensional variables, the following normalized system of equations is derived for $u$, $I$ and $A$, and for the constant $p > 0$:

(4.2a)$$\begin{gather} - \epsilon\,\frac{\mathrm{d}}{\mathrm{d} y} \left( \mu_0 p\,\frac{\mathrm{d} u/\mathrm{d} y}{|\mathrm{d} u/\mathrm{d} y|} + \mu_1 \sqrt{p\,\frac{A_0}{n}}\,\frac{\mathrm{d} u}{\mathrm{d} y} \right)= \beta_o\,|u_o - u|\,(u_o - u), \end{gather}$$

(4.2b)$$\begin{gather}I = \sqrt{\frac{A_0}{pn}}|\,\mathrm{d} u/\mathrm{d} y|, \end{gather}$$

(4.2c)$$\begin{gather}A = 1 - \phi_0 I^\alpha, \end{gather}$$

(4.2d)$$\begin{gather}\int_0^1A\,\mathrm{d} y = A_0. \end{gather}$$

Here, $\epsilon = H/L$ and $\beta _o = \rho _o/\rho _iC_o$, and the non-dimensional average floe size $\bar {d}$ is set to $\sqrt {A_0/n}$. The system (4.2) is closed by enforcing periodic boundary conditions for $u$, $I$ and $A$. Following our findings in table 2, we assume that the parameters $\mu _0$, $\phi _0$ and $\alpha$ are strictly greater than zero. We also assume that $\mu _1 > 0$ in all but § 4.2.4, where we study the case when $\mu _1 = 0$ with the intention of understanding the plastic component of the momentum equation (4.2a).

All numerical results computed in this section take (2.15) as the ocean velocity, which is written as

(4.3)\begin{equation} u_o(y) = 1 - |1 - 2y| \end{equation}

for $y\in (0,1)$ when non-dimensionalized.

4.2. The momentum equation

In order to understand the system of equations (4.2), we first focus on the momentum equation (4.2a). When considering the entire system (4.2), the pressure $p\in \mathbb {R}$ is one of the unknowns. However, it is useful to first assume it to be known, in which case we can solve the momentum equation  (4.2a) for $u$, and study the effect of $p$ on $u$. Here, we show that (4.2a) can be understood as a minimization problem. This reformulation of the momentum equation allows us to establish the existence and uniqueness of solutions. Moreover, the optimality conditions for the minimization problem result in a different formulation of the plasticity component of the rheology, which avoids the singularity, present in (4.2a), at $\mathrm {d} u/\mathrm {d} y = 0$. With this reformulation of the plastic term, we are able to find analytical solutions to the purely plastic problem that arises when $I \ll 1$, near the quasi-static regime.

4.2.1. Reformulation of (4.2a) as a minimization problem

Given a pressure $p > 0$, solutions $u$ to (4.2a) minimize the functional

(4.4)\begin{equation} \mathcal{J}(u) := \epsilon\mu_0p \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|\mathrm{d} y + \frac{\epsilon\mu_1}{2} \sqrt{p\,\frac{A_0}{n}} \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^2\mathrm{d} y + \frac{\beta_o}{3}\int_0^1|u_o - u|^3\,\mathrm{d} y \end{equation}

over the space of velocity profiles

(4.5)\begin{equation} V := \lbrace v\in H^1((0,1)): v\ \text{is a periodic function}\rbrace. \end{equation}

In the definition of $V$, the space $H^1((0,1))$ denotes the Sobolev space of weakly differentiable and periodic functions on the unit interval (Adams & Fournier 2003). As explained in Appendix B, the functional $\mathcal {J}$ is strictly convex and coercive over $V$, and therefore admits a unique minimizer. In this sense, one can state that the momentum equation (4.2a) also has a unique solution.

To derive (4.2a) for the minimizer $u$ of $\mathcal {J}$, we first note that if $u$ minimizes $\mathcal {J}$, then

(4.6)\begin{equation} \frac{1}{t}\,(\mathcal{J}(u + t(v-u)) - \mathcal{J}(u)) \geq 0 \quad \forall\ v\in V,\enspace \forall\ t\in (0,1). \end{equation}

The three terms on the right-hand side of (4.4) are convex, with the last two differentiable over all $V$. By exploiting the convexity of the first term (the $L^1$ norm) and the differentiability of the other two terms, we find that

(4.7)\begin{align} & \epsilon\mu_0p \left(\left|\frac{\mathrm{d} v}{\mathrm{d} y}\right| - \left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|\right) + \epsilon\mu_1\,\sqrt{p\,\frac{A_0}{n}} \int_0^1 \frac{\mathrm{d} u}{\mathrm{d} y}\,\frac{\mathrm{d} (v-u)}{\mathrm{d} y}\,\mathrm{d}\kern0.07em x \nonumber\\ &\quad - \beta_o \int_0^1 |u_o - u|\,(u_o - u)(v-u)\,\mathrm{d}\kern0.07em x \geq 0 \quad \forall v\in V. \end{align}

A variational statement as in (4.7) is known as a variational inequality. Under the assumption that the solution $u$ not only is in $V$ but is twice continuously differentiable, we may deduce that

(4.8a)$$\begin{gather} - \epsilon\,\frac{\mathrm{d} }{\mathrm{d} y}\left(\sigma_{xy}^P + \mu_1\, \sqrt{p\,\frac{A_0}{n}}\,\frac{\mathrm{d} u}{\mathrm{d} y}\right)= \beta_o\,|u_o - u|\,(u_o - u), \end{gather}$$

(4.8b)$$\begin{gather}|\sigma_{xy}^P| \leq \mu_0 p, \end{gather}$$

(4.8c)$$\begin{gather}\sigma_{xy}^P \left|\frac{\mathrm{d} u}{\mathrm{d} y}\right| = \mu_0 p\,\frac{\mathrm{d} u}{\mathrm{d} y} . \end{gather}$$

A derivation similar to that of (4.8) from (4.7) can be found in Glowinski, Lions & Trémoliéres (1981, § 1.3). In (4.8), we have introduced $\sigma _{xy}^P$, the purely plastic component of the shear stress. Introducing this new variable allows us to reformulate (4.2a) such that the singularity at $\mathrm {d} u/\mathrm {d} y = 0$ is removed. Indeed, if $\mathrm {d} u/\mathrm {d} y \neq 0$, then it is clear that (4.8) is equivalent to (4.2a). In this case, we have that $|\sigma _{xy}^P| = \mu _0 p$, and we say that the material has reached its plastic yield strength $\mu _0 p$. Conversely, when $\mathrm {d} u/\mathrm {d} y = 0$, (4.8a) remains well defined, unlike (4.2a). We remark that $\mathrm {d} u/\mathrm {d} y = 0$ must follow from (4.8c) whenever $|\sigma _{xy}^P| < \mu _0 p$ (the material has not reached its plastic yield strength). Below, in § 4.2.4, we provide further insight into the plastic component of the shear stress by computing purely plastic solutions to the momentum equation analytically.

4.2.2. Regularization of the plastic term to facilitate its numerical solution

The first-order optimality condition for the minimization of $\mathcal {J}$ is a variational inequality (rather than a variational equality) because the first term on the right-hand side of (4.4) (the $L^1$ norm) is non-differentiable when $\mathrm {d} u/\mathrm {d} y = 0$. We can make $\mathcal {J}$ differentiable by regularizing it as follows:

(4.9)\begin{align} \mathcal{J}_\varDelta(u) := \epsilon\mu_0p \int_0^1\sqrt{\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^2 + \varDelta^2}\,\mathrm{d} y + \frac{\epsilon\mu_1}{2}\,\sqrt{p\,\frac{A_0}{n}} \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^2\mathrm{d} y + \frac{\beta_o}{3}\int_0^1|u_o - u|^3\,\mathrm{d} y, \end{align}

where $\varDelta > 0$ is a small parameter. The first-order optimality conditions for the minimization of $\mathcal {J}_\varDelta$ over $V$ corresponds with the equation

(4.10)\begin{equation} {-}\epsilon\,\frac{\mathrm{d}}{\mathrm{d} y}\left(\mu_0 p\,\frac{\mathrm{d} u/\mathrm{d} y}{\sqrt{|\mathrm{d} u/\mathrm{d} y|^2 + \varDelta^2}} + \mu_1\,\sqrt{p\,\frac{A_0}{n}}\,\frac{\mathrm{d} u}{\mathrm{d} y} \right)= \beta_o\,|u_o - u|\,(u_o - u). \end{equation}

Although the system (4.8) can be solved numerically by e.g. introducing a Lagrange multiplier (Glowinski et al. 1981), it is easier to solve (4.10). This is the strategy that we use for solving the momentum equation, and, as we explain below in § 4.3.1, the complete system (4.2). To do so, we use the finite element method (FEM) implemented in Firedrake (Ham et al. 2023). In particular, we approximate the velocity profile $u$ with continuous piecewise linear functions. In figure 5, we plot solutions to (4.10) for two different values of $p$ and a range of $\varDelta > 0$. Convergence of the velocity profiles as $\varDelta \to 0$ is clearly visible in these plots; in fact, for $\varDelta \leq 10^{-2}$, the solutions become indistinguishable. We remark that if we remove the viscous component of the rheology in the regularized equation (4.10), then we essentially arrive at Hibler's model in one dimension, given below by (6.1a).

Figure 5. Solutions $u$ to the momentum equation (4.10) for different values of the ice pressure $p$ and the numerical regularization parameter $\varDelta$. We set $u_o$ equal to (2.15), normalized with $[u] = u_{o,max}$, and $\epsilon = 2\times 10^5$, $A_0 = 0.8$ and $n = 2000$. We observe convergence to a solution as $\varDelta \to 0$ and (a) for small $p$, $u \to u_o$, and (b) for large $p$, $u\to 0.5$, as expected from the theory.

4.2.3. Velocity profiles in the limits of small and large pressures

The pressure or ice strength $p$ is a fundamental variable in the continuum model; understanding its effect on $u$ is fundamental for making sense of our sea ice model. Figure 5 suggests that for small $p$, the velocity profile $u$ approaches the ocean's $u_o$, and for large $p$, $u$ flattens and comes close to a constant-valued velocity profile. We can deduce this behaviour from the functional $\mathcal {J}$. For small values of $p$,

(4.11)\begin{equation} \lim_{p\to 0} \mathcal{J}(u) = \frac{\beta_o}{3}\int_0^1|u_o - u|^3\,\mathrm{d} y, \end{equation}

therefore since $u$ minimizes $\mathcal {J}$, it must follow that $u \to u_o$. On the other hand, for large values of $p$, we see that

(4.12)\begin{equation} \mathcal{J}(u) \approx \epsilon\mu_0p \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|\mathrm{d} y \quad \text{as}\ p\to \infty, \end{equation}

and in principle, any constant velocity profile $u_c \in \mathbb {R}$ minimizes (4.12). However, this constant velocity field is constrained by the total force balance of the system. That is, due to the periodic boundary conditions, if we integrate (4.8a) along $(0,1)$, then we must have that

(4.13)\begin{equation} \beta_o \int_0^1 |u_o - u|\,(u_o - u)\,\mathrm{d} y = 0. \end{equation}

Therefore, the constant value to which $u$ tends as $p\to \infty$ will be a solution to (4.13) with $u = u_c \in \mathbb {R}$. In § 4.2.4, we show that a critical pressure $p_c$ can be found such that $u$ is constant whenever $p > p_c$ and $I \ll 1$.

4.2.4. Purely plastic solutions to the momentum equation

In figure 4, we can see that $\mu (I)$ approximately becomes constant for small inertial numbers $I \ll 1$, such that $\mu (I) \approx \mu _0$ and the flow rheology is plastic. This regime is closely related to the quasi-static regime for granular media, with the material behaving like a purely plastic flow characterized by a critical state at which plastic deformation occurs (Wood 1990).

The momentum equation for the purely plastic problem where $\mu (I) = \mu _0$ is given by

(4.14)\begin{equation} \epsilon\,\frac{\mathrm{d} \sigma_{xy}^P}{\mathrm{d} y} ={-} \beta_o\,|u_o - u|\,(u_o - u). \end{equation}

Equation (4.14) must be complemented with (4.8b) and (4.8c). Here, we present a method for calculating purely plastic solutions. Additionally, we find a critical pressure $p_c$ such that for $p > p_c$, the velocity profiles $u$ that solve (4.14) remain constant, and no shear strain occurs in the sea ice. In § 5, we show that this critical pressure approximates the pressure computed from the DEM when the global sea ice concentration is high. When following the derivation of purely plastic solutions, it is helpful to look at their plots in figure 6.

Figure 6. Solutions $u$ to the purely plastic momentum equation (4.14). We use the same parameter values as those used for figure 5, and $\varDelta = 10^{-2}$. The numerical solutions (black) to the purely plastic problem are indistinguishable from the analytical solution given by (4.20) (red dashed lines). As the ice pressure $p$ increases, the velocity profiles flatten at the critical pressure $p_c$ given by (4.19), and $p_c \approx 13.2$ in this case.

Conditions (4.8b) and (4.8c) for the plastic stress tensor indicate that there exist two distinct regions of the flow field: a region where the sea ice has yielded and $|\sigma _{xy}^P| = p\mu _0$, and another region where the ice has not yielded and $|\sigma _{xy}^P| \leq p\mu _0$ and $\mathrm {d} u/\mathrm {d} y = 0$. By working with this distinction, we can find a purely plastic solution to (4.14). Due to the symmetry of the problem, we assume that $\sigma _{xy}^P = 0$ at $y = 0$, $1/2$ and $1$. Then integrating (4.14) along the interval $(0,y)$ for some $y\in (0,1)$, we find that

(4.15)\begin{equation} \sigma_{xy}^P(y) ={-}\frac{\beta_o}{\epsilon}\int_0^y|u_o - u|\,(u_o - u)\,\mathrm{d} y. \end{equation}

Since $\sigma _{xy}^P(0) = 0$, we must necessarily have an interval $(0,y_1)$ where the ice has not yielded and the velocity equals a constant $u_1$. In the context of the ocean velocity profile (4.3), it makes sense to assume that $u_1 > 0$, and therefore $u_1 \geq u_o$ near $y = 0$, so that

(4.16)\begin{equation} \sigma_{xy}^P(y) = \frac{\beta_o}{6\epsilon}\,((2y - u_1)^3 + u_1). \end{equation}

Since the material has not yielded for $y\in (0,y_1)$, we have that $|\sigma _{xy}^P(y)| < \mu _0 p$. Equation (4.16) tells us that $\sigma _{xy}^P(y)$ increases with $y$ over this interval; this means that

(4.17)\begin{equation} \lim_{y\to y_1,\ y < y_1} \sigma_{xy}^P(y) = \mu_0 p, \end{equation}

and we find that

(4.18)\begin{equation} u_1 = \left(\frac{6\epsilon\mu_0p}{\beta_0}\right)^{1/3}. \end{equation}

Clearly, an upper bound is needed for $u_1$ in (4.18) because it grows indefinitely with $p$, yet it is senseless for the ice to circulate at speeds greater than the maximum ocean velocity when a steady state has been reached. We can make sense of this paradox by first assuming the existence of an interval $(y_1, y_2)$, where $y_1 < y_2 < 1/2$, in which the sea ice has yielded and $\sigma _{xy}^P = \mu _0 p$. In this interval, we must have that $u= u_o$ because $\sigma _{xy}^P$ is constant and therefore the ocean drag is zero. This means that $y_1 = u_1/2$. Moreover, repeating the same argument as that used for deriving (4.18), we assume that $u = u_2$ for some constant $u_2 < 1$ on $(y_2, 1/2)$, and find that $u_2 = 1 - u_1$ and $y_2 = u_2/2$. Now, the assumption that $y_1 < y_2 < 1/2$ will hold only for values of $p$ for which $u_1 \leq u_2$, i.e. $u_1 \leq 1/2$, and this upper bound on $u_1$ defines a critical pressure $p_c$ given by

(4.19)\begin{equation} p_c = \frac{\beta_o}{48\epsilon\mu_0}. \end{equation}

For $p \geq p_c$, the integral force balance along the domain must hold (see (4.13)); as a result, $u_1$ can be at most equal to $1/2$. Putting these results together, we may write the analytical solution to (4.14) as

(4.20a)$$\begin{gather} u(y) = \left\{\begin{array}{@{}ll@{}} u_1, & 0< y<\dfrac{u_1}{2}, \\ u_o(y), & y_1< y<\dfrac{1}{2} - y_1, \\ 1 - u_1, & \dfrac{1}{2} - y_1< y<\dfrac{1}{2}, \end{array}\right.\quad \text{if}\ p < p_c, \end{gather}$$

(4.20b)$$\begin{gather}u(y) = \frac{1}{2} \quad \text{if}\ p \geq p_c, \end{gather}$$

for $y\in [0,1]$. We test the validity of (4.20) by showing that it is indistinct from its numerical counterpart in figure 6. We compute this numerical solution by regularizing the shear stress $\sigma _{xy}^P$ as in (4.10), and setting $\varDelta = 10^{-3}$.

4.3. Solutions to the complete continuum model

In § 4.2, we have seen that, given a value for the pressure $p$, we can solve the momentum equation and find a velocity profile $u$ for the sea ice. We also prove that solutions to the momentum equation must exist and be unique. However, in general, the pressure $p$ is one of the unknowns in the system of equations (4.2), together with the sea ice concentration $A$ and the inertial number $I$. Here, we first present a numerical method for solving (4.2), and show that solutions to this system always exist and, under some circumstances, are unique.

4.3.1. A numerical method for the complete model (4.2)

In order to solve the system (4.2), we follow the approach discussed in § 4.2.2 for solving the momentum equation. There, the regularized equation (4.10) is solved numerically with the FEM. When solving the complete system (4.2), we find that also regularizing the inertial number $I$ improves the convergence properties of our numerical solver substantially. Therefore, we solve for

(4.21a,b)\begin{equation} I_\varDelta := \sqrt{\frac{A_0}{pn}\left(\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^2 + \varDelta^2\right)} \quad \text{and} \quad A = 1 - \phi_0I_\varDelta^\alpha \end{equation}

instead of (4.2b) and (4.2c). Then we use the FEM to solve the system of equations given by the regularized momentum equation (4.10), (4.21a,b), and the constraint for global concentration (4.2d). We approximate $u$ with continuous piecewise linear functions, and $I_\varDelta$ and $A$ with piecewise constant functions. Our solver is implemented in Firedrake (Ham et al. 2023) in such a way that the complete nonlinear system is solved with Newton's method. For small values of $\varDelta$, Newton's method tends to fail unless a very good initial guess for the solution $(u,I_\varDelta,A,p)$ is given. For this reason, we find that solving for a sequence of decreasing values of $\varDelta$, using the solution for the previous $\varDelta$ as the initial guess for the next $\varDelta$, yields a robust solution method.

We test the sensitivity of solutions $(u,I_\varDelta,A,p)$ to changes in $\varDelta$ by solving the system for values of $\varDelta$ between $10^{-3}$ and $10$, and for $A_0 = 0.75$ and $A_0 = 0.9$. The numerical results, which are plotted in figure 7, indicate that the solutions become more sensitive to $\varDelta$ as $A_0$ increases; this is natural, since $I$ decreases with $A_0$, and the plasticity term becomes more important. It is also clear from these plots that although the velocity profiles $u$ come very close to convergence for the smallest values of $\varDelta$, the local concentration profiles still experience visible changes around the symmetry points $y = 0$, $0.5$ and 1. In these points, the shear strain drops to 0, and by the definition of $I$, we expect $A = 1$ there. However, as we argue in § 5, when comparing the continuum model and the DEM, we consider such drastic changes in the local concentration to be artificial. This argument motivates the use of $\varDelta$ not just as a numerical parameter that helps us to solve the system numerically, but as a parameter that improves the model and may have a physical significance.

Figure 7. Solutions $u$ and $A$ to the continuum model (4.2) with (a,b) $A_0 = 0.75$, and (c,d) $A_0 = 0.9$, regularized with a parameter $\varDelta$ as explained in § 4.3.1. (e) Plot of the variation of the pressure $p$ with $\varDelta$. We use the same parameter values as those used for figure 5.

4.3.2. Existence and uniqueness of solutions

We conclude the analysis of the continuum model by showing that at least one solution to (2.13) must exist, and whenever $u_o$ is given by (2.15), the solution is unique. To do so, we first reformulate (2.13) solely in terms of $u$ and $p$ by eliminating $A$ and $I$ as follows. Substituting (4.2c) into (4.2d) yields

(4.22)\begin{equation} \int_0^1 I^\alpha\,\mathrm{d} y = \frac{1 - A_0}{\phi_0}. \end{equation}

Then by substituting the definition of $I$ in (4.2b) into the integrand in the expression above, we arrive at

(4.23)\begin{equation} \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^\alpha\mathrm{d} y = \frac{1 - A_0}{\phi_0}\left(p\,\frac{n}{A_0}\right)^{{\alpha}/{2}}. \end{equation}

Therefore, the system of equations (4.2) is equivalent to solving the momentum equation (4.2a) together with the constraint (4.23) for $u$ and $p$. We can interpret this problem as the minimization of the functional (4.4) over a set of velocity profiles subject to the constraint (4.23). Next, we define $\mathcal {F}:\mathbb {R}_+\to \mathbb {R}_+$ by

(4.24)\begin{equation} \mathcal{F}(p) := \int_0^1\left|\frac{\mathrm{d} u}{\mathrm{d} y}\right|^\alpha\mathrm{d} y, \end{equation}

with $u$ denoting the solution to (4.2a) with the pressure set to $p$; this operator is well-defined because, given a pressure $p > 0$, a unique velocity profile $u$ solves the momentum equation (4.2a). We also define $\mathcal {C}:\mathbb {R}_+\to \mathbb {R}_+$, given by the left-hand side of (4.23), that is,

(4.25)\begin{equation} \mathcal{C}(p) := \frac{1 - A_0}{\phi_0}\left(p\,\frac{n}{A_0}\right)^{{\alpha}/{2}}. \end{equation}

It is then clear that a pressure $p$ is part of the solution to the system of equations (4.2) if and only if

(4.26)\begin{equation} \mathcal{C}(p) =\mathcal{F}(p). \end{equation}

In § 4.2.3, we show that the velocity profile $u$ has two distinct limit behaviours. We find that $u$ approaches $u_o$ as $p \to 0$, and that $\mathrm {d}u/\mathrm {d} y$ tends to $0$ as $p\to \infty$. This means that

(4.27a,b)\begin{equation} \lim_{p\to 0} \mathcal{F}(p) = \int_0^1\left|\frac{\mathrm{d} u_o}{\mathrm{d} y}\right|^\alpha\mathrm{d} y \quad \text{and} \quad \lim_{p\to \infty} \mathcal{F}(p) = 0. \end{equation}

On the other hand, the function $\mathcal {C}$ is strictly increasing, with $\mathcal {C}(0) = 0$. Therefore, whenever $\int _0^1| \mathrm {d} u_o/\mathrm {d} y|^\alpha \,\mathrm {d} y > 0$ (i.e. $u_o$ is not a flat velocity profile) and $\mathcal {F}$ is a continuous function, we must have at least one solution $p$ to (4.26). Moreover, if $\mathcal {F}$ is a decreasing function, then this pressure must be unique.

In figure 8, we plot the functions $\mathcal {C}$ and $\mathcal {F}$ for several values of $A_0$ and an ocean velocity profile given by (4.3). In this case, we have that

(4.28)\begin{equation} \int_0^1\left| \frac{\mathrm{d} u_o}{\mathrm{d} y}\right|^\alpha\mathrm{d} y = 2^\alpha, \end{equation}

and the function $\mathcal {F}$ appears to be strictly decreasing. This is expected, because an increase in pressure is accompanied by a flattening of the sea ice velocity profile. This means that there is a unique $p$ for which (4.26) holds; since there is only one velocity profile $u$ that solves the momentum equation (4.2a), it follows that the solution to the continuum model (4.2) must be unique.

Figure 8. Functions $\mathcal {F}(p)$ and $\mathcal {C}(p)$, defined in (4.24) and (4.25), respectively, for different values of $A_0$. The lighter tones of blue and red are associated with lower values of $A_0$, although $\mathcal {F}$ hardly changes with $A_0$. A pressure $p$ solves (2.13) whenever $\mathcal {C}(p) =\mathcal {F}(p)$ (circles). To generate these curves numerically, the same problem setting as that of figure 5 is used.

5. Comparison of the DEM with the continuum model

The continuum model (4.2) is designed with the objective of capturing the averaged behaviour of the sea ice simulations carried out with the DEM. Here, we verify that in the one-dimensional setting of the steady ocean periodic problem, the continuum model is indeed capable of replicating most of the results of the DEM. We remind the reader that the DEM solutions can be considered steady and one-dimensional only in the sense that for spatially and temporally averaged data, we can expect variations in time and in the $x$-direction to be negligible. At a small scale, we certainly expect the data to present variations in time and in the $x$-direction, and the vertical velocity of the sea ice to be non-zero. Moreover, DEM computations are initialized using a random floe initialization; however, as displayed in figure 2, the (spatially and temporally averaged) steady states that the DEM computations reach from different floe initializations are almost indistinguishable.

Throughout this section, we use the parameters in tables 1 and 2 when solving the continuum model. The continuum model is solved with the FEM as explained in § 4.3.1, setting $\varDelta = 10^{-3}$ and using a uniform mesh of 300 cells.

5.1. Variation in global concentration $A_0$ and ocean speed $u_{o,max}$

We first evaluate the continuum model's accuracy in replicating the DEM results used for fitting the functions $\mu (I)$ and $\varPhi (I)$ in § 3. These results are computed for the ocean velocity profile (2.15) and ice thickness $H = 2\,{\rm m}$. We consider six global sea ice concentrations $A_0$ between $0.7$ and $0.95$, and four ocean velocities $u_{o,max}$ between $0.1$ and $1\,{\rm m}\,{\rm s}^{-1}$. For each of these cases, we run the DEM until a steady state is reached, then we extract 10 values of the sea ice velocity and concentration along the $y$-direction, as explained in § 3, and one value for the pressure $p$, averaged over the whole square patch of ocean.

Figure 9 displays the velocity and sea ice concentration profiles $u$ and $A$, and the pressure $p$, for both the continuum model and the DEM. The velocity profile $u$ and the pressure $p$ are normalized as explained in § 4.1 using $[u] = u_{o,max}$. A consequence of this normalization is that the non-dimensional solutions $u$, $A$ and $p$ to the continuum model (4.2) are indifferent to the value of $u_{o,max}$. For most of the DEM results, this is also the case. For each value of $A_0$, almost all of the normalized velocity profiles (figures 9a,c,e,g,i,k) and pressure points (figure 9m) appear to collapse onto a single curve or point, which is well approximated by the continuum model.

Figure 9. Solutions to the non-dimensional continuum model (4.2) (black lines) compared with results from the DEM (markers), with the ocean velocity in blue. For each pair of panels in the first three rows, we fix the global sea ice concentration $A_0$, vary the maximum ocean velocity $u_{o,max}$ for the DEM, and plot the velocities in (a,c,e,g,i,k) and the local concentrations in (b,d, f,h,j,l). In (m), we plot the pressure $p$ in terms of $A_0$. Due to the normalization in terms of $[u] = u_{o,max}$, the solutions to the continuum model are indifferent to a change in $u_{o,max}$.

Departures from the other normalized DEM results are most visible for slow ocean currents and low concentrations. This becomes particularly clear when looking at the pressure in figure 9(m) when $u_{o,max} = 0.1\,{\rm m}\,{\rm s}^{-1}$ and $A_0 \leq 0.8$, where the pressure values from the DEM (circles) depart substantially from the prediction of the continuum model (black solid line). These mismatches signal the continuum model's limitations to capture the DEM results for the range of regimes considered. When fitting the dilatancy function $\varPhi (I)$, the largest misfit is also found for points of slow ocean currents and low concentration ($u_{o,max} = 0.1\,{\rm m}\,{\rm s}^{-1}$ and $A_0 \leq 0.75$); see figure 4(b). If the DEM results are to approximate a continuum model as in (4.2), then we expect the equality

(5.1)\begin{equation} \sigma_{xy}(y) ={-} \int_0^y t_{ox}\,\mathrm{d} y \end{equation}

to hold approximately for the DEM results, with $t_{ox}$ denoting the horizontal component of the ocean drag. If we denote by $\sigma _{xy,i}$ and $t_{ox,i}$ the values extracted from the DEM at the grid cell of width $\Delta y$ located at $y_i$ (see § A.2 for an overview of how these quantities are obtained), then a discrete balance analogous to (5.1) is

(5.2)\begin{equation} \sigma_{xy,i} ={-} \Delta y \sum_{j = 1}^i t_{ox,j}. \end{equation}

We remark that the term $t_{ox,j}$ results from averaging the ocean drag on each ice floe, as opposed to introducing the averaged sea ice velocity into (2.3). We also note that (5.2) ignores any choice of rheology, and should hold independently of our choice of functions $\mu (I)$ and $\varPhi (I)$. A surprising result that we find when investigating the DEM data is that the terms on the left- and right-hand side of the equality in (5.2) differ in orders of magnitude whenever the sea ice concentration is low and the ocean currents are slow; see figure 10(a). Conversely, for faster ocean currents and denser concentrations, these terms become approximately equal (see figures 10c–f). Inertial effects are found to be negligible in all of the cases that we consider; moreover, since DEM quantities are averaged along the whole length of the domain in the $x$-direction, the term corresponding to $\partial \sigma _{xx}/\partial x$, which should be considered in a two-dimensional setting, becomes zero. Therefore, this finding raises the questions of whether a fundamentally different continuum model should be used for low sea ice concentrations and slow ocean currents, or if a continuum model is valid at all.

Figure 10. Shear stress $\sigma _{xy,i}$ and integrated ocean drag $\int t_{ox,i} := \Delta y \sum _{j = 1}^i t_{ox,j}$ extracted from the DEM for (a,b,c) $A_0 = 0.7$ and (d,e, f) $A_0=0.9$, and $u_{o,max}$ values (a,d) $0.1\,{\rm m}\,{\rm s}^{-1}$, (b,e) $0.5\,{\rm m}\,{\rm s}^{-1}$ and (c, f) $1\,{\rm m}\,{\rm s}^{-1}$.

The DEM local concentration $A$ is accurately captured by the continuum model for low sea ice concentrations; see figures 9(b,d, f,h). For $A_0 = 0.90$ and 0.95 (figures 9j,l), the continuum model overestimates the degree of dilatancy, although the general trends are visibly similar, with regions of higher concentration around $y=0$, $1/2$ and $1$, where the strain rates are lowest. Since an increase in $\varDelta$ diminishes the local variations in $A$ (see figure 7), this suggests that $\varDelta$ could be adjusted to improve the fit with the data. In this case, it would enter the model as a physical parameter whose interpretation should be explored further.

In figure 9(m), we also test two limit approximations of the pressure that we find in § 4 when $u_o$ is given by (4.3). When $p \ll 1$, which occurs for the smaller values of $A_0$, we expect $\mathcal {C}(p) \approx 2^\alpha$ because $u \approx u_o$. This implies that

(5.3)\begin{equation} p \approx 4\,\frac{A_0}{n}\left(\frac{\phi_0}{1 - A_0}\right)^{2/\alpha}\quad \text{as}\ A_0 \to 0. \end{equation}

On the other hand, high values of $A_0$ result in $I \ll 1$ (see figure 4), therefore $\mu (I) \approx \mu _0$, leading to the purely plastic regime studied in § 4.2.4. Figure 9(k) shows that the velocity profiles are mostly flat in this region, suggesting that the critical pressure $p_c$ calculated in (4.19) is a good approximation of $p$. By plotting (5.3) and $p_c$ in figure 9(m), we find that these are indeed good approximations in their respective limits.

5.2. Variation in ice thickness and number of floes

Next, we test the effectiveness of the continuum model in capturing the DEM results for different ice thicknesses and floe sizes. In figure 11, we plot results for the DEM and the continuum model for steady states with $H= 0.5$ and $4\,{\rm m}$, different ocean velocities $u_{o,max}$, and a global sea ice concentration $A_0 = 0.8$. The ice thickness $H$ enters the continuum model via the parameter $\epsilon = H/L$ in the momentum equation (4.2a). An increase in $H$ is accompanied by an increase in $\epsilon$, which effectively acts as a decrease in the normalized external forcing. We therefore expect an increase in $H$ to result in a decrease in the shear strain and pressure. As observed in figures 11(a,c), this is the case for the velocity profiles of both the continuum model and the DEM; an increase in $H$ is accompanied by a flattening of the normalized velocity profiles. The continuum model is once again accurate in capturing the averaged velocities of the DEM and the general trends in the variation of the sea ice concentration. The pressure resulting from the continuum model and the DEM also decreases as the ice thickness increases (see figure 11e); for the pressure, we find a decent fit between the DEM and the continuum model.

Figure 11. Solutions to the non-dimensional continuum model (4.2) (black lines) compared with results from the DEM (markers). For each pair of panels in the first row, we fix the sea ice thickness $H$ and vary the maximum ocean velocity $u_{o,max}$ (the ocean velocity is plotted in blue in (a,c)). In (e), we plot the pressure $p$ in terms of $H$. Due to the normalization in terms of $[u] = u_{o,max}$, the solutions to the continuum model are indifferent to a change in $u_{o,max}$.

Finally, figure 12 contains a comparison between the DEM and the continuum model for different numbers of floes $n$. In particular, we present results with $n = 500$, 2000 and 5000 for $A_0 = 0.8$ and $0.9$, and $u_{o,max} = 0.5\,{\rm m}\,{\rm s}^{-1}$. An increase in $n$ implies a decrease in the effective viscosity of the model, and we therefore expect a decrease in the ice strength or pressure $p$; in the limit where $n\to \infty$, the rheology becomes purely plastic, this can be seen by taking this limit in (4.2a). In figure 12, one can see that the pressure in the continuum model indeed decreases with $n$, but this is not the case with the DEM. In fact, the results of the DEM appear to change little with $n$, especially for $A_0 = 0.9$. Despite this, the results from the DEM do not depart from those of the continuum model excessively.

Figure 12. Solutions to the non-dimensional continuum model (4.2) (lines) compared with results from the DEM (markers). The solutions to the continuum model are plotted in red for $n = 500$, blue for $n = 2000$, and green for $n = 5000$; the tone is set to light if $A_0 = 0.8$, and dark if $A_0 = 0.9$ (see legend in (e)). For each pair of panels in the first row, we fix the sea ice thickness $H$ and vary the maximum ocean velocity $u_{o,max}$ (the ocean velocity is plotted in blue in (a,c)). In (e), we plot the pressure $p$ in terms of $n$.

6. Comparisons with existing models for sea ice

The continuum model studied in this paper shares features with existing models for sea ice. Here, we examine those similarities and also establish differences with our continuum model. Due to its ubiquity in sea ice modelling, we first consider Hibler's model in § 6.1, before considering other models in § 6.2.

6.1. Hibler's model

The most widely used model for sea ice is Hibler's model, which was first presented in Hibler (1979) and treats sea ice like a viscoplastic material. Under the one-dimensional conditions of the steady square ocean patch and the non-dimensionalization in § 4, Hibler's model reduces to the form

(6.1a)$$\begin{gather} - \epsilon\,\frac{\mathrm{d} }{\mathrm{d} y}\left(\frac{p}{2e}\, \frac{\mathrm{d} u/\mathrm{d} y}{\sqrt{|\mathrm{d} u/\mathrm{d} y|^2 + \delta^2}}\right)= \beta_o\,|u_o - u|\,(u_o - u), \end{gather}$$

(6.1b)$$\begin{gather}p = \frac{P^\ast}{\rho_iu_{o,max}^2}\exp{({-}20(1-A))}. \end{gather}$$

In (6.1), the parameter $P^\ast$ is an empirical constant, $\delta$ is a regularization parameter, and $e = 2$ represents the eccentricity of the elliptical yield curve from which Hibler's model is derived. That is, Hibler's model assumes a plastic rheology based on an elliptical yield curve that is then regularized. Inside the plane of principal stresses, the yield curve is set in such a way that the ice resists only compression, not pure extension. This geometrical configuration is motivated by the observation that in pack ice, deformation occurs through ridging (compression) and the opening of leads (extension); of these two mechanisms, only ridging requires a non-negligible amount of plastic work (Rothrock 1975).

In the steady one-dimensional setting, (6.1a) can be recovered from the regularized momentum equation (4.10) by setting $\mu _0 = 1/(2e)$ and $\mu _1 = 0$ (pure plasticity). In fact, $1/(2e) = 0.25$, which comes very close to the value $\mu _0 = 0.26$ that we infer from the DEM. This means that close to the quasi-static regime, when $I \ll 1$ and $\mu (I) \approx \mu _0$, we essentially work with the momentum equation from Hibler's model. This is a very reasonable coincidence, because Hibler's model was designed for the central ice pack where the sea ice concentration is very high.

The main departure between our model and Hibler's is the expression for the pressure (6.1b): as $u_{o,max}$ increases, $p$ decreases in Hibler's model, unlike our continuum model, where $p$ is independent of $u_{o,max}$. This makes it impossible for Hibler's model to capture the invariance of the non-dimensional sea ice velocity and pressure with $u_{o,max}$, which we find in most of the DEM results in § 5. This is made clear in figure 13, where we show solutions to Hibler's model and compare them with the DEM for different values of $A_0$ and $u_{o,max}$; an increase in $u_{o,max}$ is accompanied by excessively large changes in the velocity profile. We set $P^\ast = 5\times 10^4\,{\rm N}\,{\rm m}^{-1}$ and $\delta = 0.1$ in order to get a good fit with the DEM for some values of $u_{o,max}$. We remark that Hibler's model is designed for ice that fractures and ridges; in this setting, an increase in the ocean drag weakens the ice through this mechanical deformation. Below, in § 7, we state that future extensions of our continuum model should include the effects of fracturing and ridging, available in SubZero. These future investigations should study the validity of (6.1b) when fracturing and ridging are accounted for.

Figure 13. Solutions to Hibler's model for different global concentrations $A_0$ and ocean velocities $u_{o,max}$. We also plot the results from the DEM as a reference. Here, we set $H = 2\,{\rm m}$ and $n = 2000$.

The choice of an elliptical yield curve in Hibler's model is motivated by the simplicity of the resulting rheological formulation, yet other possibilities consistent with the requirement of null resistance to pure extension are available, such as the parabolic lens and the teardrop yield curves (Zhang & Rothrock 2005; Ringeisen, Losch & Tremblay 2023), and the Mohr–Coulomb criterion (Ip et al. 1991). In fact, as mentioned at the end of § 3, the purely plastic rheology in the $\mu (I)$ formulation considered here (given by (3.2) with $\mu _1 = 0$) can be written as the Mohr–Coulomb rheology presented in Ip et al. (1991) and Gutfraind & Savage (1997). This becomes clear if we note that $\|\boldsymbol {{S}}\| = {D}_1 - {D}_2$, where ${D}_1$ and ${D}_2$ denote the principal components (eigenvalues) of $\boldsymbol {{D}}$, and we compare the plastic component in (3.2) with expression (6) in Gutfraind & Savage (1997). A significant difference between the two expressions is that in (6) from Gutfraind & Savage (1997), the viscosity is capped to a maximum value in order to avoid the inevitable blow-up that occurs with a purely plastic rheology; this is another form of regularization. This difference of plastic yield curves between Hibler's model and the purely plastic version of our model is another point of departure between both models in a two dimensional setting.

6.2. Other models

In § 3, we explain that the viscous component in (3.2) can be interpreted as the rheological effects of collisions, which become increasingly important as $I$ increases. As mentioned in that section, a collisional rheology is derived in Shen et al. (1987) that is also of a linearly viscous nature. In fact, this collisional rheology establishes a viscosity that is linearly proportional to $\bar {d} \sqrt {\rho _i H p}$, as in (3.2) (see Herman (2022) for a clear description of the collisional rheology). For this, the model proposed in Feltham (2005) comes close to (3.2), since it takes the ice rheology to be the sum of Hibler's rheology and the collisional rheology from Shen et al. (1987).

To the knowledge of the authors, the only other study using the $\mu (I)$ rheology to model sea ice is Herman (2022). There, the floes are driven by a moving wall (as opposed to an ocean or atmospheric current, as in our case), and the DEM is based on disk-shaped floes, with a much more severe polydispersity. Two main differences can be found between the function $\mu (I)$ that we infer (figure 4) and the one found in Herman (2022). First, although both cases consider a very similar range of $I$ values, the magnitude of the friction $\mu$ differs considerably, although it is of the same order of magnitude. Second, the $\mu (I)$ curve in Herman (2022) plateaus for $I > 10^{-1}$ – something that we do not observe. It remains unclear what may cause these differences. Regarding the second point, we justify our use of a linear function for $\mu$, as in (3.1a,b), by remarking that it simplifies the resulting constitutive equation, enabling a more thorough analysis of the model, and resembles the $\mu$ function proposed in Herman (2022) over a large range of $I$ values.

7. Conclusions and future work

In this paper, we have presented a novel mathematical model for sea ice in the MIZ based on the $\mu (I)$ rheology. We have inferred the form of this rheology in § 3 using data produced with DEM computations. With the analysis in § 4, we prove that the steady one-dimensional formulation of this model, given by (4.2), is well-posed in the sense that it has a unique solution. The numerical results presented in § 5 demonstrate that this model is capable of replicating most of the results of the DEM in the context of steady one-dimensional problems. The most visible departure from the continuum model occurs for low ocean velocities and sea ice concentrations, with $u_{o,max} = 0.1\,{\rm m}\,{\rm s}^{-1}$ and $A_0 \leq 0.85$; in this case, the DEM results indicate that the shear stress no longer balances the integrated ocean drag, signalling a breakdown of the underlying grid-averaged momentum balance equation. That is, since inertial effects are found to be negligible in the steady states that we consider, and the extensional stresses cancel out when integrating across the $x$-direction of the domain to perform the averaging of DEM quantities, our basic continuum model, prior to any choice of rheology, establishes a balance between shear stresses and ocean drag. In figure 10, this balance is found to hold approximately for all cases except those of slow ocean currents and low sea ice concentrations.

The lack of validity of the continuum model for slow ocean currents and low concentrations found in § 5 should be explored further, since this is a regime that we expect to find in parts of the MIZ (Stewart et al. 2019). As explained in § 5, this lack of accuracy is accompanied by the breakdown of the balance between the averaged shear stress and the integrated ocean drag extracted from the DEM. This balance precedes the choice of rheology and therefore indicates that either a fundamentally different continuum model should be used or some assumption necessary for a continuum model to even hold is no longer valid.

Mechanical processes such as fracturing, ridging and welding are fundamental processes in sea ice dynamics (Feltham 2008). Therefore, future improvements of the model presented here should consider the inclusion of these effects. Since ridging and fracturing are fundamental processes in the central ice pack (Rothrock 1975), our continuum model could also be valid in this area when these effects are included, yielding a unified sea ice model. Since SubZero includes these mechanical interactions between ice floes (Manucharyan & Montemuro 2022), it is a promising tool for exploring their macroscopic effects on the rheology. A first step could be to explore the consequences of including floe fracturing. We expect that this would set an upper bound on the pressure and shear stress that the ice cover can sustain. This bound will probably be closely related to the fracture criterion used at the ice floe level. Moreover, floe fracturing will create smaller floes that lead to higher degrees of polydispersity.

The comparisons between the continuum model and the DEM in § 5 are confined to the context of the steady periodic ocean problem, which is essentially one-dimensional. Any future studies must examine the accuracy of the $\mu (I)$ rheology in modelling sea ice dynamics under unsteady conditions and in a two-dimensional configuration. Barker et al. (2015) discovered the $\mu (I)$ formulation presented in § 2.1 to be mathematically ill-posed in time-dependent problems, and Schaeffer et al. (2019) has proposed modifications that avoid these instabilities while leaving the steady equations unchanged. A future investigation in the context of sea ice should take these studies into account.

The rheology that we propose in (3.2) is local in the sense that the viscosity and the pressure at a certain point of the domain depend only on other quantities and their derivatives at that same point. Yet granular materials create complex contact networks that enable the interaction of grains set far apart. This leads to non-local effects, and extensions of the $\mu (I)$ rheology that include these effects have been proposed in the context of granular flows (Kamrin & Koval 2012; Bouzid et al. 2013). Such effects will probably also arise in sea ice modelling and should be considered in future extensions of the model that we propose here.