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 (Reference Hibler III1979) or Gray & Morland (Reference Gray and Morland1994). 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. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005; Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006; Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). 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 (Reference Savage1984) 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. (Reference Jop, Forterre and Pouliquen2006) 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. (Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005). 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. (Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005) and Herman (Reference Herman2022) 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.

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 Reference Adams and Fournier2003). 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 (Reference Glowinski, Lions and Trémoliéres1981, § 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. Reference Glowinski, Lions and Trémoliéres1981), 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. Reference Ham2023). 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 Reference Wood1990).

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. Reference Ham2023) 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.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.1. Hibler's model

The most widely used model for sea ice is Hibler's model, which was first presented in Hibler (Reference Hibler III1979) 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 Reference Rothrock1975).

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 Reference Zhang and Rothrock2005; Ringeisen, Losch & Tremblay Reference Ringeisen, Losch and Tremblay2023), and the Mohr–Coulomb criterion (Ip et al. Reference Ip, Hibler and Flato1991). 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. (Reference Ip, Hibler and Flato1991) and Gutfraind & Savage (Reference Gutfraind and Savage1997). 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 (Reference Gutfraind and Savage1997). A significant difference between the two expressions is that in (6) from Gutfraind & Savage (Reference Gutfraind and Savage1997), 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. (Reference Shen, Hibler III and Leppäranta1987) 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 (Reference Herman2022) for a clear description of the collisional rheology). For this, the model proposed in Feltham (Reference Feltham2005) 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. (Reference Shen, Hibler III and Leppäranta1987).

To the knowledge of the authors, the only other study using the $\mu (I)$ rheology to model sea ice is Herman (Reference Herman2022). 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 (Reference Herman2022). 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 (Reference Herman2022) 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 (Reference Herman2022) over a large range of $I$ values.