# Analytical Study of an Isotropic Viscoplastic Sea Ice Model in Idealized Configurations

Jérôme Sirven, Bruno Tremblay

► To cite this version:

Jérôme Sirven, Bruno Tremblay. Analytical Study of an Isotropic Viscoplastic Sea Ice Model in Idealized Configurations. Journal of Physical Oceanography, 2015, 45 (2), pp.331-354. 10.1175/JPO-D-13-0109.1. hal-01135447

HAL Id: hal-01135447

https://hal.science/hal-01135447

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical Study of an Isotropic Viscoplastic Sea Ice Model in Idealized Configurations

JÉRÔME SIRVEN

UMR Locean Laboratory (Sorbonne Universités-Université Pierre et Marie Curie, CNRS, IRD, MNHN),
Paris, France

BRUNO TREMBLAY

Department of Atmospheric and Oceanic Sciences, McGill University, Montreal, Quebec, Canada

(Manuscript received 17 May 2013, in final form 27 February 2014)

## ABSTRACT

Analytic solutions of a mechanical sea ice model are computed in idealized configurations. They are then used to study the properties of this model. It classically assumes that the ice behaves at large scale as an isotropic viscoplastic medium. The plastic regime is characterized by a Mohr–Coulomb yield curve. The flow rule corresponds to the one used in granular mediums and depends on a parameter $\delta$ that characterizes the expansion properties of the medium. Using simple model configurations, this study first shows that a sliding of the ice along the coast must be permitted; otherwise, the model generally has no solution when the plastic regime is active. This study then shows that the viscous regime is reached only if the stress remains nearly uniform over a large area. For a stress having no particular properties, the plastic regime acts everywhere. In this case, the compressive stress may reach the maximum value allowed by the model close to the coastline. The extension of the domain where the compressive stress is at its maximum depends on $\delta$ and the direction of the forcing field. Over this domain, the ice behaves as a fluid material with a small negative viscosity. Last, the authors found that neither the existence of the solution nor its unicity are guaranteed in this stationary model. This result does not imply that the unicity is lost in the transient problem; it suggests that the evolution of sea ice depends not only on the forcing, but also on the initial conditions or history of the system.

## 1. Introduction

Sea ice is formed by the freezing of seawater. However, this process alone is not sufficient to explain how the ice fields are transformed during their life. Mechanical processes are important to understand how the combined action of temperature changes at the surface of the ice, wind, and ocean succeeds in modifying the young ice to create the pack ice that covers the polar oceans.

Sea ice is not a simple elastic solid whose properties have been well known for more than a century. At large spatial scale (larger than 100 km), the Arctic Ice Dynamics Joint Experiment campaign suggested that sea ice could be treated as an isotropic continuous medium (Coon et al. 1974). These assumptions have been recently

questioned (see, e.g., Coon et al. 2007; Rampal et al. 2008), and anisotropic (Hibler and Schulson 2000) or discontinuous (Schreyer et al. 2006) models have been developed. However, most of the models used for climate studies—a wide range of such models are described and their performances discussed in the Arctic Ocean Model Intercomparison Project (AOMIP; see Proshutinsky et al. 2011)—still assume continuity and isotropy. Among these models, those based on a viscous–plastic rheology are the most frequently used; this includes the so-called elastic–viscous–plastic models (Hunke and Dukowicz 1997), which have the same physics as the viscous–plastic models but differ by their numerical implementation. Models that represent the physical elastic stress are more rarely used because of the complexity of the numerical resolution. For instance, see Pritchard (2001) who explains how the elastic waves, which contaminate the elastic–plastic numerical simulations, can be eliminated.

Sea ice models generally assume that the ice behaves like a very viscous fluid when the strain rate is small

<footnote>*Corresponding author address:* Jérôme Sirven, Sorbonne Universités (UPMC, Univ. Paris 06)-CNRS-IRD-MNHN, Locean Laboratory, Université Pierre et Marie Curie, F-75005, Paris, France.
E-mail: jerome.sirven@locean-ipsl.upmc.fr</footnote>
(under a threshold value arbitrarily fixed). Beyond this value, the ice is assumed to behave like a plastic medium: the two principal stresses satisfy a relation that defines the “yield curve” (also known as the plasticity criterion) and the strain rate tensor is linked to the stress tensor by a “flow rule,” the most common one being the “normal flow rule.”

Such assumptions are exceptional—even unique to our knowledge—in the field of continuous medium mechanics. Indeed, when the stress remains small, the materials commonly studied in mechanics (often metals) behave as elastic media, of which the theory is well established. The theoretical investigations are therefore advanced for the static behavior of elastoplastic materials [see Hill (1960) for a classic book or Salençon (2002) for a more recent one], whereas they remain inexistent for viscoplastic materials. Moreover, laboratory experiments have been performed for a long time (more than a century) in order to determine the validity of the assumptions that are made to study the plastic regime of metals [Hill (1960) provides many historical details]. Such studies are just beginning for sea ice (Schulson and Nickolayev 1995; Schulson 2001) and suggest that sea ice can show a brittle behavior when the strain rates are high (see also Marsan et al. 2004).

Despite these difficulties or discrepancies, the viscoplastic models are widely used and turn out to be efficient to predict the thickness, concentration, and drift of sea ice. Obviously, improvements are needed, and a large body of work has been devoted to this effort [see Proshutinsky et al. (2011) and the publications listed in their article and associated with the AOMIP project]. The current search around the viscoplastic rheology of sea ice is mainly devoted to determine the best yield curve and flow rule (e.g., Tremblay and Mysak 1997; Hibler and Schulson 2000; Ukita and Moritz 1995; Wilchinsky and Feltham 2004; Taylor et al. 2006) or to look for more stable and accurate numerical models (e.g., Hunke and Dukowicz 1997; Heil and Hibler 2002; Hunke and Dukowicz 2002; Lemieux et al. 2008; Lemieux et al. 2010). However, little attention has been paid to the mathematical properties of the models themselves and their solutions. Only Gray (1999) investigated the well-posedness of the equations of the viscous–plastic sea ice model of Hibler (1979) and showed that, for a uniaxial divergent flow, the equations could be unstable and ill posed because of the loss of hyperbolicity.

This paper investigates the mathematical properties and solutions of viscous–plastic sea ice models, but adopts a more basic approach than that followed by Gray (1999) and considers models different from the one of Hibler (1979). Indeed, stationary exact solutions are searched in very simple configurations for isotropic models with a Mohr–Coulomb yield curve. This choice is not the most

usual in sea ice modeling but is now considered in a few papers (see, e.g., Heil and Hibler 2002; Sedlacek et al. 2007) and is simple enough to allow analytic computations. Moreover, recent studies (Weiss et al. 2007) based on the observations of Richter-Menge et al. (2002) give it an experimental support.

Section 2 presents the model we study [similar to that of Tremblay and Mysak (1997)]. In section 3, the idealized configuration we consider is introduced and the corresponding solutions are computed and studied when the system behaves as a plastic medium. The following section considers the simple and exceptional case of a uniform free drift, where a viscous behavior may be observed in a part of the domain. A discussion of the results concludes the paper.

## 2. Isotropy and Coulomb friction law in sea ice models

The sea ice models usually are one of many components of climate models, which include at least atmosphere, ocean, and land models. Here, these components are not considered, and the physical variables that they allow to compute (for example the wind stress or the current velocity in the ocean) will be prescribed. The focus is on the mechanical processes that play a role in the evolution of sea ice. This approach is somewhat artificial, since they interact with thermodynamical processes; however, we make this assumption to isolate the internal stresses and strains in the sea ice from the other variables, and, by doing so, we hope to clarify the properties of sea ice models.

### a. Free drift of sea ice

The acceleration and advection terms can be neglected in comparison with other terms in the momentum equation when long time scales (longer than about 1 day) are considered (see, e.g., Thorndike 1986; Leppäranta 2005). The free drift of sea ice $\mathbf{u}_{fd}$ is then determined by a balance between the Coriolis force, the friction with the atmosphere and the ocean, and the gradient of dynamic height at the ocean surface. This leads to the equation

$$ \rho_i h f \mathbf{k} \times \mathbf{u}_{fd} = -\rho_i h g \mathbf{V}(H_d) + \boldsymbol{\tau}_a - \boldsymbol{\tau}_w, $$

where $\boldsymbol{\tau}_a$ and $\boldsymbol{\tau}_w$ are the stresses exerted respectively by the wind and the ocean on the ice, $H_d$ is the dynamic height of the ocean, $h$ is the ice thickness, $\rho_i$ is the ice density, $f$ is the Coriolis parameter assumed constant, $\mathbf{k}$ is a unitary vector normal to the ice surface, and $g$ is the gravitational acceleration. The ice concentration $A$ is not taken into account in this equation for simplicity, but could be introduced by multiplying $\boldsymbol{\tau}_a - \boldsymbol{\tau}_w$ by $A$. In the
<figure><img src="image_2.png" /><figcaption>FIG. 1. Variation of (top) $\mathcal{N}$ and (bottom) $\theta_0$ as functions of the ice thickness $h$.</figcaption></figure>

free drift regime, the ice concentration is very low, and it is important to include it correctly in the external forcing. Here, we will only consider regimes where the internal stress is high; in this case, its effect is very weak according to Tremblay and Mysak (1997).

The stress $\tau_a$ exerted by the wind on the ice surface can be assumed independent of the ice drift because the wind velocity is usually much larger than the ice velocity. On the contrary, the ocean drag on the ice flow $\tau_w$ depends on the ocean currents as well as on the ice drift. Following the formulation proposed by Hutchings et al. (2005), we set $\tau_w = \rho_w c_w \mathbf{R}_{\theta_w}(\mathbf{u}_{fd} - \mathbf{u}_g)$, where $\mathbf{u}_g$ is the geostrophic ocean current, $\theta_w$ is a fixed positive rotation factor, and $\mathbf{R}_{\theta_w}(\cdot)$ is the rotation matrix of angle $\theta_w$. Quadratic drag laws (McPhee 1975) are more correct and more frequently used; the linear law is used in the following for simplicity. Rewriting $\rho_i h g \mathbf{V}(H_d)$ in term of geostrophic velocity, we obtain

$$ \rho_i h f \mathbf{k} \times (\mathbf{u}_{fd} - \mathbf{u}_g) + \rho_w c_w \mathbf{R}_{\theta_w} (\mathbf{u}_{fd} - \mathbf{u}_g) = \boldsymbol{\tau}_a. $$

Introducing the matrix

$$ \mathbf{M} = \begin{bmatrix} \rho_w C_w \cos\theta_w & -\rho_i f h - \rho_w C_w \sin\theta_w \\ \rho_i f h + \rho_w C_w \sin\theta_w & \rho_w C_w \cos\theta_w \end{bmatrix}, $$

the previous equation becomes $\mathbf{u}_{fd} = \mathbf{u}_g + \mathbf{M}^{-1}\boldsymbol{\tau}_a$. To facilitate the interpretation of this relation, we write $\mathbf{M} = \mathcal{N}\mathbf{R}_0$, where $\mathcal{N} = \sqrt{(\rho_i h f)^2 + 2\rho_i \rho_w h f C_w \sin\theta_w + C_w^2 \rho_w^2}$ is a positive number and $\mathbf{R}_0$ a rotation matrix of angle $\theta_0$:

$$ \mathbf{R}_0 = \begin{bmatrix} \cos\theta_0 & -\sin\theta_0 \\ \sin\theta_0 & \cos\theta_0 \end{bmatrix}, $$

with $\cos\theta_0 = \rho_w C_w \cos\theta_w/\mathcal{N}$ and $\sin\theta_0 = (\rho_i h f + \rho_w C_w \sin\theta_w)/\mathcal{N}$. Hence, the free drift velocity can be written as

$$ \mathbf{u}_{fd} = \mathbf{u}_g + \mathcal{N}^{-1} \mathbf{R}_0^{-1} \boldsymbol{\tau}_a. \quad (1) $$

The variations of the coefficient $\mathcal{N}$ and angle $\theta_0$ as a function of the ice thickness $h$ are shown in Fig. 1 (the values of the constants are given in Table 1). The coefficient $\mathcal{N}$ increases with $h$, since the influence of the wind becomes less important when the ice becomes thicker. For the same reason, the turning angle $\theta_0$ increases with $h$.

b. *The assumptions about rheology: Isotropy and Coulomb friction law*

The material behavior of sea ice in viscous-plastic models is characterized by the stress tensor integrated over the thickness of the ice $\boldsymbol{\sigma}$. Its component $\sigma_{ij}$ acts on a plane perpendicular to the $i$ axis in the $j$ direction. The stress tensor is a function of the only strain rate tensor $\dot{\epsilon}$, which suffices to characterize the deformation of sea ice (the strain tensor is ignored). Both tensors being symmetric, they can be diagonalized after a rotation of the coordinate axes $x$ and $y$. For a common material, two distinct rotations are needed; but for an isotropic material, the same rotation can be used. In the rotated coordinate axes these tensors take the form
<table><caption>TABLE 1. Constants used in the numerical computations.</caption><thead><tr><th>Variable</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Water density</td><td>ρ<sub>w</sub></td><td>1025 kg m<sup>-3</sup></td></tr><tr><td>Ice density</td><td>ρ<sub>i</sub></td><td>900 kg m<sup>-3</sup></td></tr><tr><td>Air drag coefficient</td><td>C<sub>a</sub></td><td>2.3 × 10<sup>-3</sup> m s<sup>-1</sup></td></tr><tr><td>Water drag coefficient</td><td>C<sub>d</sub></td><td>0.66 × 10<sup>-3</sup> m s<sup>-1</sup></td></tr><tr><td>Water turning angle</td><td>θ<sub>w</sub></td><td>25°</td></tr><tr><td>Coulombic friction angle</td><td>φ</td><td>30°</td></tr><tr><td>Ice thickness</td><td>h</td><td>1.8 m</td></tr><tr><td>Maximal compressive stress</td><td>P<sub>max</sub></td><td>7. × 10<sup>3</sup> kg s<sup>-2</sup></td></tr><tr><td>Maximal viscosity</td><td>η<sub>max</sub></td><td>1. × 10<sup>12</sup> kg s<sup>-1</sup></td></tr></tbody></table>

$$ \boldsymbol{\sigma} = \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix} \quad \text{and} \quad \dot{\boldsymbol{\epsilon}} = \begin{bmatrix} \dot{\epsilon}_1 & 0 \\ 0 & \dot{\epsilon}_2 \end{bmatrix}, \qquad (2) $$

where $\sigma_1$, $\sigma_2$, $\dot{\epsilon}_1$, and $\dot{\epsilon}_2$ are the principal stresses and strain rates. As shown in appendix A, the isotropy hypothesis thus implies

$$ \left\{ \begin{aligned} \sigma_{xx} &= -P + \eta(\partial_x u - \partial_y v) \\ \sigma_{xy} &= \sigma_{yx} = \eta(\partial_x v + \partial_y u), \\ \sigma_{yy} &= -P - \eta(\partial_x u - \partial_y v) \end{aligned} \right. \qquad (3) $$

with

$$ P = -(\sigma_1 + \sigma_2)/2, $$

$$ \eta = S/D, $$

$$ S = (\sigma_1 - \sigma_2)/2, \quad \text{and} $$

$$ D = \dot{\epsilon}_1 - \dot{\epsilon}_2 = +\sqrt{(\partial_x u - \partial_y v)^2 + (\partial_x v + \partial_y u)^2}. $$

A more usual way of writing Eq. (3) is $\sigma_{ij} = -P\delta_{ij} - \eta\text{Tr}(\dot{\boldsymbol{\epsilon}})\delta_{ij} + 2\eta\dot{\epsilon}_{ij}$, where $\delta_{ij}$ is the Kronecker symbol and Tr($\mathbf{A}$) is the trace of a matrix $\mathbf{A}$. This form makes explicit a difference between the model described here and the standard (Hibler) viscous-plastic model. In the latter, the relation $\sigma_{ij} = -[P + (\zeta - \eta)\text{Tr}(\dot{\boldsymbol{\epsilon}})]\delta_{ij} + 2\eta\dot{\epsilon}_{ij}$ is used. The second viscosity $\zeta$ does not appear in our model; it is included in the unknown $P$.

Two supplementary relations must be specified to determine $P$ and $S$. When the deformation of sea ice is plastic, the results of Weiss et al. (2007)—the second figure of their article is reproduced here in Fig. 2— coming from data of the Surface Heat Budget of the Arctic Ocean (SHEBA) and Sea Ice Mechanics Initiative (SIMI) experiments (Richter-Menge et al. 2002) lead to adopt a coulombic rheology; it inherits its name from the Coulomb’s friction law that characterizes the sliding of an object on a plane. It assumes an affine

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2. Figure reproduced from Weiss et al. (2007). Stress states recorded during the Arctic winter at the sensor Baltimore from mid-October 1997 to end of June 1998 (1 measure per hour). The top panel is plotted in a principal stress space $\sigma_1, \sigma_2$, and the bottom one is rotated 45°.</figcaption></figure>

relation between $\sigma_1$ and $\sigma_2$: $\sigma_1 = a\sigma_2 + b$. This relation defines a wedgelike yield curve; the constant $b$ allows us to take into account a possible tensile strength. We rewrite this relation as follows:

$$ S = \frac{\sigma_1 - \sigma_2}{2} = \left(-\frac{\sigma_1 + \sigma_2}{2} - P_0\right) \sin\phi = (P - P_0) \sin\phi. \qquad (4) $$

The constant $P_0$ represents the maximum tensile strength, and the angle $\phi$ is the effective angle of friction that links the shear and normal stresses acting on an object sliding on a plane.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Schema describing the mechanical states of the sea ice in the model. When the plastic-compressive regime is operating, the flow rule $\dot{\epsilon}_1 + \dot{\epsilon}_2 = (\dot{\epsilon}_1 - \dot{\epsilon}_2) \tan\delta$ no longer applies. The flow rule with $\tan\delta$ and the normal flow rule are represented by the dashed and continuous bold lines crossing the axes $\dot{\epsilon}_1$ and $\dot{\epsilon}_2$.</figcaption></figure>

A second relation must still be specified. It adds to the isotropy hypothesis and introduces a supplementary relation between the stress tensor $\sigma$ and the strain rate tensor $\dot{\epsilon}$. When the yield curve is elliptic, the most usual one is given by a normal flow rule: it assumes that the strain rate is normal to the stress on the yield curve. Here, this would lead to the relation $\dot{\epsilon}_1 + \dot{\epsilon}_2 = (\dot{\epsilon}_1 - \dot{\epsilon}_2) \sin\phi$. Normal flow rules result in attractive mathematical properties when the yield curve is convex—such as a guarantee of positive dissipation—and have been validated by experimental studies in metallurgy (see, e.g., Hill 1960). However, they are not often used with Mohr–Coulomb yield curves and the data of Stern and Moritz (2002) suggest that they are inappropriate in this case to represent sea ice rheology (see below); indeed they lead to too much dilatation (divergence associated with shear).

Here, we will assume in agreement with Tremblay and Mysak (1997) that the strain rate makes a constant angle (distinct from $\pi/2$) with the stress. This leads to replacing $\sin\phi$ by a coefficient, $\tan\delta$, so that $\dot{\epsilon}_2 + \dot{\epsilon}_1 = (\dot{\epsilon}_1 - \dot{\epsilon}_2) \tan\delta$. In the x, y coordinates, the divergence and the shear deformation are thus linked by the relation

$$ \partial_x u + \partial_y v = \tan\delta \sqrt{(\partial_x u - \partial_y v)^2 + (\partial_x v + \partial_y u)^2}. \quad (5) $$

The system is therefore under divergence as soon as $\tan\delta > 0$. Note that all the relations are given in this section for $\sigma_1 - \sigma_2 > 0$ and $\dot{\epsilon}_1 - \dot{\epsilon}_2 > 0$ (part of Figs. 2 and 3 located under the axis of symmetry $\sigma_1 = \sigma_2$).

Relation (5) was obtained by Tremblay and Mysak (1997), who were able to give it a physical basis. They considered a sea ice model based on a granular material rheology. Assuming that the granules (ice floes here) are small, moving planar objects (similar to hockey pucks of various diameters), the small-scale sliding of the individual granules may occur in a direction that differs from the large-scale sliding of the ice pack. Consequently, to model this process, they introduced an angle between the mean small-scale sliding planes and the large-scale sliding plane; they called it the dilatancy angle $\delta$. We adopt this formulation below.

Stern and Moritz (2002) studied the sea ice kinematics from the Canadian RADARSAT satellite over the SHEBA site and showed that the divergence is usually much smaller than the shear deformation: for example, in autumn and early winter, the ratio between these two values is smaller than 0.14 at the 200-km scale. Consequently, to avoid unrealistic divergence or convergence of the ice velocity, it is reasonable to assume that the angle $\delta$ is bounded by $\pm 10^\circ$ ($\tan\delta = \pm 0.176$). Considering
the results of Stern and Moritz (2002), the normal flow rule is clearly wrong. However, since it is frequently used when the yield curve is elliptic, we briefly investigate this unrealistic case. It corresponds to $tanδ = sinφ$ (for $φ = 30°$, $sinφ = 0.5$, and $δ = 24°$).

With these hypotheses, the plastic regime is characterized by Eqs. (3), (4), and (5). The open circles on the thin continuous lines limiting the triangle in Fig. 3 correspond to this regime.

Other states of the system may be characterized by points that are inside the triangle (symbolized by a + sign in Fig. 3). These points verify $S < (P - P_0)sinφ$ and $S > -(P - P_0)sinφ$ ($P > P_0$). The model assumes that ice flow is viscous in this case, following the choice made in the vast majority of publications. We thus set $η = η_m$, where $η_m$ is the maximum authorized value of the coefficient $η$. Since Eq. (3) may be based on a physical ground linked to the geometry of the ice floes (Tremblay and Mysak 1997), we assume that it is still verified in this case.

Last, the compressive stress cannot indefinitely increase in sea ice. To take this into account, various solutions have been proposed. For example, Hibler and Schulson (2000) closed the two straight lines associated with the coulombic rheology by an elliptic curve when $P$ exceeded a critical value. The resulting yield curve had the shape of an ice cream cone (see also Coon et al. 1974). Others' solutions have been suggested leading to yield curves similar to a triangle (Tremblay and Mysak 1997), a lemon (Taylor et al. 2006), and so on.

Here, we try to define a model as close as possible to the observations. The experimental data from Richter-Menge et al. (2002) (see Fig. 2) show that the set of stress states is not convex; consequently, the yield curve that enfolds them is also not convex (this contrasts with the usual hypothesis of convexity made in plasticity theory). Considering these experimental results, we abandon the idea of closing the yield curve. We assume that $P - P_0$ can reach the maximum value $P_{max}$ only at the extremity of the segments defining the plastic behavior (stars in Fig. 3), following a path that becomes tangent to these segments close to the extreme points. As a consequence, the Coulomb friction law (4) is still valid when $P - P_0 = P_{max}$, but Eq. (5) is not now applicable. It will be seen below that the flow rule is determined by the system itself; in particular the relation $\dot{\epsilon}_1 > \dot{\epsilon}_2$ is no longer verified when $\sigma_1 > \sigma_2$. This point is discussed in more detail in section 3.

### c. The sea ice momentum and continuity equations

Our main aim consists in understanding the consequences of the assumptions made about rheology in

section 2b. We thus assume that the sea ice thickness $h$ is known and constant.<fnref n="1" /> Figure 1 allows us to assess the impact of this hypothesis on the matrix $M$. It shows that the hypothesis remains reasonable when $h$ varies in a range smaller than about 2 m, that is, over a large surface of the Arctic Ocean. However, it is poor when $h$ sustains strong variations; indeed, $N$ varies from about 0.7 to 1.4 and $θ$ from 25° to 65° when $h$ increases from 0 to 8 m. In this case, the quantitative results are no longer valid, but the nature of the differential equations we study [see below, Eq. (6)] is not modified. Consequently, the qualitative results are not affected.

Everywhere the internal stress is negligible the velocity of the ice $u$ is equal to the free drift velocity $u_{fd}$. When contact forces appear, this equality is no longer verified and the Newton second law specifies how it is modified:

$$ \mathcal{N}\mathbf{R}_0(\mathbf{u}_i - \mathbf{u}_{fd}) = \mathbf{\nabla} \cdot \boldsymbol{\sigma}. \quad (6) $$

As a consequence of isotropy, the divergence of the stress tensor is given by

$$ \mathbf{\nabla} \cdot \boldsymbol{\sigma} = \left\{ \begin{array}{l} -\partial_x P + \partial_x [\eta(\partial_x u - \partial_y v)] + \partial_y [\eta(\partial_y u + \partial_x v)] \\ -\partial_y P + \partial_y [\eta(\partial_y u + \partial_x v)] - \partial_x [\eta(\partial_x u - \partial_y v)] \end{array} \right\} \quad (7) $$

[see Eqs. (3)]. The supplementary hypotheses that have been done (a coulombic rheology, a “granular flow rule,” and maximum authorized values for $η$ and $P$) lead to distinguish three states in the system:

* In the first state (circles in Fig. 3), sea ice behaves as a plastic material. We have

$$ \eta = \frac{(P - P_0) \sin\phi}{\sqrt{(\partial_x u - \partial_y v)^2 + (\partial_y u + \partial_x v)^2}} \quad (8) $$

associated with the dilatancy equation [Eq. (5)]. We refer to this case as the “plastic-shear regime” below.

* In the second state, the Coulomb friction law is abandoned as soon as $η$ becomes larger than $η_m$ (symbol + in Fig. 3). It is replaced by the relation $η = η_m$, and the dilatancy equation is kept unchanged. We refer to this case as the “viscous regime” below.

* Last, if $P - P_0$ reaches the maximum value $P_{max}$, the dilatancy equation is abandoned and $η$ is equal to

<footnote n="1">¹ A dynamic-thermodynamic model include a mass conservation equation; our hypothesis means that the unknown sources and sinks of mass would adjust in such a model, in such a way that $h$ remains constant.</footnote>
$\pm P_{\max} \sin\phi/\sqrt{(\partial_x u - \partial_y v)^2 + (\partial_y u + \partial_x v)^2}$. This state is characterized by the two stars in Fig. 3. We refer to this case as the “plastic-compressive regime” below. Note that this case is the consequence of a plastic behavior. However, we distinguish it from the first case and give it another name for clarity. Indeed, it differs from the plastic-shear regime on two points: the compressive stress becomes uniform and the flow rule, which led to Eq. (5), is no longer verified.

In the first two cases (plastic shear and viscous behavior), a system of three differential equations is obtained, which allows us to compute the three variables $P, u$, and $v$. In the third case, the pressure $P$ is known and two equations suffice to compute $u$ and $v$. The apparent simplicity of the problem (only three unknowns in a two-dimensional space) hides a considerable technical difficulty: the limit of the domain where each regime applies is unknown and thus must be computed. In the following section, we consider a simplified setup, which will make this computation possible.

### 3. Translation invariant solution of the model

To study the model described in the previous section, the geometry of the domain is first simplified; we first assume that the coast is a straight line that defines the y axis. Then, we suppose that the external forcings $u_g$ and $\tau_a$ yield a free drift $u_{fd} = (u_{fd}, v_{fd})$ that depends only on the distance $x$ to the coast [see Eq. (1)]. Rather than prescribing $u_g$ and $\tau_a$, we choose to prescribe $u_{fd}$. This is coherent with the hypothesis that the ice thickness is known and uniform ($\mathcal{N}$ and $\mathbf{R}_0$ are known).

With these simplified geometry and forcing terms, the model admits a solution $\{P(x), \mathbf{u} = [u(x), v(x)]\}$ that depends only on $x$. This solution is determined by the vector equation

$$ \mathcal{N}\mathbf{R}_0(\mathbf{u} - \mathbf{u}_{fd}) = \partial_x(-P + \eta\partial_x u)\mathbf{e}_x + \partial_x(\eta\partial_x v)\mathbf{e}_y \quad (9) $$

$(e_x$ and $e_y$ are unitary vectors respectively perpendicular and parallel to the coast), and supplementary equations that depend on the state of the system (see section 2c and Fig. 3). These equations are summarized below.

* Plastic-shear regime: $\eta = (P - P_0) \sin\phi/\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$ and $\partial_x u = \tan\delta\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$. The name “dilatancy equation” will hereinafter refer to this last equation, which is a particular case of the general dilatancy Eq. (5).

* Viscous regime: $\eta = \eta_m$ and $\partial_x u = \tan\delta\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$.

* Plastic-compressive regime: $P - P_0 = P_{\max}$ and $\eta = (\sigma_1 - \sigma_2)/[2(\dot{\epsilon}_1 - \dot{\epsilon}_2)] = \pm(P - P_0) \sin\phi/\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$. A different sign for $\sigma_1 - \sigma_2$ and $\dot{\epsilon}_1 - \dot{\epsilon}_2$ is here allowed (symbol $\pm$).

#### a. Plastic-shear regime

The dilatancy equation $\partial_x u = \tan\delta\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$ is easily solved. We rewrite it as $(\partial_x u)^2(1 - \tan^2\delta) = (\partial_x v)^2 \tan^2\delta$ that leads to

$$ v = \epsilon \frac{\sqrt{1 - \tan^2\delta}}{\tan\delta} u + v_i, \quad (10) $$

where $v_i$ is an arbitrary constant and $\epsilon = \pm 1$; $\partial_x u$ and $\tan\delta$ have the same sign, which adds a constraint on the sign of $\epsilon$. Note that the case $\delta = 0$ leads to the solution $u = u_0$ and therefore is meaningless here (see the discussion about $\delta$ in the conclusion).

Using the dilatancy equation, the expression of $\eta$ (first bullet in section 3), and Eq. (10), we have $\eta\partial_x u = (P - P_0)\sin\phi\tan\delta$ and $\eta\partial_x v = \epsilon(P - P_0)\sin\phi\sqrt{1 - \tan^2\delta}$, where the dilatancy angle $\delta$, which models the microscopic properties of the medium, is a specified constant. Consequently, the momentum equation takes the form

$$ \mathcal{N}\mathbf{R}_0(\mathbf{u} - \mathbf{u}_{fd}) = \frac{\partial P}{\partial x}(-\mathbf{e}_x + \sin\phi\mathbf{t}), \quad (11) $$

where $\mathbf{t} = \tan\delta\mathbf{e}_x + \epsilon\sqrt{1 - \tan^2\delta}\mathbf{e}_y$ is a unitary vector that depends only on $\delta$, and $P_0$ no longer appears. Equation (11) is satisfied only if $\mathbf{u} - \mathbf{u}_{fd}$ is colinear to the vector $\mathbf{K} = (K_x, K_y) = \mathbf{R}_0^{-1}(-\mathbf{e}_x + \sin\phi\mathbf{t})$, which is obtained from the vector $-\mathbf{e}_x + \sin\phi\mathbf{t}$ after a rotation of angle $-\theta_0$. Noting $\mathbf{K}^{\perp} = (-K_y, K_x)$, this condition becomes $\langle\mathbf{u} - \mathbf{u}_{fd}, \mathbf{K}^{\perp}\rangle = 0$, where $\langle\mathbf{a}, \mathbf{b}\rangle$ represents the scalar product of $\mathbf{a}$ and $\mathbf{b}$. Substituting the expression of $\mathbf{u} = u\mathbf{e}_x + v\mathbf{e}_y$ and using Eq. (10) yield

$$ u\langle\mathbf{t}, \mathbf{K}^{\perp}\rangle + v_i \tan\delta\langle\mathbf{e}_y, \mathbf{K}^{\perp}\rangle = \tan\delta\langle\mathbf{u}_{fd}, \mathbf{K}^{\perp}\rangle. $$

A supplementary condition should be added to determine the value of the constant $v_i$. If the domain where the ice has a plastic behavior extends up to the coast, the component of the velocity perpendicular to the coast $\mathbf{u}$ must vanish. A more restrictive condition cannot be imposed. This means that the ice may slip along the coast. If this domain ends along the unknown line $x = x_0 > 0$, the regime changes beyond this line. As sea ice on both sides of the line $x = x_0$, it seems natural to assume that $\mathbf{u}$ is continuous at $x = x_0$.

A fastidious computation (see appendix A) leads to the following expression for $\mathbf{u}$:

$$ \mathbf{u} = \frac{\langle \mathbf{u}_{fd} - \mathbf{u}_{fd}^{0}, \mathbf{K}^{\perp} \rangle}{\langle \mathbf{t}, \mathbf{K}^{\perp} \rangle} \mathbf{t} + \frac{\langle \mathbf{u}_{fd}^{0}, \mathbf{K}^{\perp} \rangle}{\langle \mathbf{e}_y, \mathbf{K}^{\perp} \rangle} \mathbf{e}_y + u^0 \frac{\mathbf{K}}{\langle \mathbf{K}, \mathbf{e}_x \rangle}. \quad (12) $$

Here, $\mathbf{u}_{fd}^0$ and $u^0$ represent the free drift velocity and the component along the $x$ axis of the drift velocity $\mathbf{u}$ at the
<figure><img src="image_6.png" /><figcaption>FIG. 4. (top to bottom) Free drift velocity, corresponding sea ice velocity, compressive stress (kg s⁻²) (the allowed maximum is indicated by the dashed line), and viscosity coefficient η (kg s⁻¹). All the fields are represented as a function of the distance from the coast (km). The solid black rectangle represents the land.</figcaption></figure>

unknown abscissa $x = x_0$. This expression is used just below with $x_0 = 0$—the plastic-shear regime occupies the whole domain in this case—and will be used in the next subsection to connect the plastic-shear and plastic-compressive regimes.

Multiplying Eq. (11) by $\mathbf{t}^{\perp}$, an explicit expression of $\partial_x P$ is obtained:

$$ \mathcal{N} \frac{\langle \mathbf{t}^{\perp}, \mathbf{u} - \mathbf{u}_{\text{fd}} \rangle}{\langle \mathbf{t}^{\perp}, \mathbf{K} \rangle} = \partial_x P. $$

Using Eq. (12), and the fact that the compressive stress vanishes at $x = L$, abscissa of the ice edge, $P$ is given by the relation

$$ \frac{P(x)}{\mathcal{N}} = \int_{L}^{x} \left( \frac{\langle \mathbf{u}_{\text{fd}}^{0} - \mathbf{u}_{\text{fd}}, \mathbf{t}^{\perp} \rangle}{\langle \mathbf{t}^{\perp}, \mathbf{K} \rangle} + \frac{u^{0} - \langle \mathbf{u}_{\text{fd}}^{0}, \mathbf{e}_{x} \rangle}{\langle \mathbf{e}_{x}, \mathbf{K} \rangle} \right) dx \quad (13) $$

(see appendix A for more details about the computations). Relations (12) and (13) are valid only if $\eta < \eta_{\text{max}}$
and $P - P_0 < P_{max}$. Since $\eta = (P - P_0) \sin\phi \tan\delta/\partial_x u$, this leads with Eq. (12) to the constraint

$$ \eta = (P - P_0) \sin\phi \frac{\langle \mathbf{t}, \mathbf{K}^{\perp} \rangle}{\langle \partial_x \mathbf{u}_{\text{fd}}, \mathbf{K}^{\perp} \rangle} < \eta_{\text{max}}. $$

Consequently, the plastic-shear regime is active as soon as the gradient of the free drift velocity is large enough; more precisely it suffices that the inequality $||\partial_x \mathbf{u}_{\text{fd}}|| \ge (P - P_0)/\eta_{\text{max}}$ is verified. For a compressive stress close to $P_{max}$, this means that spatial variations of free drift velocity of about 1 cm s⁻¹ over a distance not exceeding 400 km lead to a plastic behavior of sea ice. This qualitative analysis thus suggests that the viscous regime is only possible in areas where the free drift velocity is nearly uniform.

Two numerical examples are given in order to illustrate these results. We once more insist on the fact that the free drift is prescribed and used to force the model (see the beginning of the section). Consequently, the fields shown in Figs. 4 and 5 are computed from Eqs. (12) and (13). As the plastic-shear regime applies over the whole domain, the boundary condition at $x = 0$ is $u = 0$. This means that $u^0$ is set to 0 in Eqs. (12) and (13). Last, the dilatancy angle $\delta$ is set to +10°.

In Fig. 4 (first panel), the free drift points toward the coast and decreases linearly from 1.8 cm s⁻¹ at 1000 km from the coast to 0.4 cm s⁻¹ at the coast. Over the whole domain, $P$ remains smaller than $P_{max}$ (we took $P_0 = 0$) and $\eta$ is smaller than $\eta_{max}$ (Fig. 4, third and fourth panel). This shows that the regime is everywhere “plastic shear.” As expected, the component of the ice velocity parallel to the coast does not vanish at $x = 0$, inducing a sliding of sea ice along the coast.

When the free drift increases, the response of the model is dramatically modified. To illustrate this change, Fig. 5 shows the compressive stress when the free drift velocity is multiplied by a factor 5 in comparison with the previous experiment. The stress $P$ exceeds the threshold value $P_{max}$ as soon as the distance to the coast becomes smaller than 870 km. The plastic-compressive regime is therefore reached in an area whose size must be determined (see just below). As the coefficient $\eta$ remains everywhere smaller than $\eta_{max}$ (not shown), the viscous regime is not observed.

b. Plastic-compressive regime: $P - P_0 = P_{max}$

First of all, we recall that the maximum compressive stress $P_{max}$ does not depend on $x$; consequently, $\partial_x P_{max} = 0$. This agrees with the hypothesis that the ice thickness $h$ is constant.

In the plastic-compressive regime, the “viscosity coefficient” $\eta = (\sigma_1 - \sigma_2)/[2(\dot{\epsilon}_1 - \dot{\epsilon}_2)]$ is equal to $\pm P_{max} \sin\phi/\sqrt{(\partial_x u)^2 + (\partial_x v)^2}$. Differing from the previous

<figure><img src="image_7.png" /><figcaption>FIG. 5. Compressive stress (kg s⁻²) when the free drift is multiplied by 5 in comparison with the previous figure. The allowed maximum is indicated by the dashed line; this solution is thus no longer valid and must be modified by taking into account the plastic-compressive regime.</figcaption></figure>

case, negative values for $\eta$ are now authorized; the reason for this paradox will be explained below.

Since the compressive stress $P$ is known, the evolution of the system is given by the unique vector equation

$$ \mathcal{N}\mathbf{R}_0(\mathbf{u} - \mathbf{u}_{\text{fd}}) = \partial_x[(\eta\partial_x u)\mathbf{e}_x + (\eta\partial_x v)\mathbf{e}_y]. \quad (14) $$

Equation (14) cannot be solved by only analytical methods. However, a bit of mathematical calculus helps understand its properties and greatly facilitates its numerical resolution. Moreover, it also provides exact solutions when the free drift is constant.

The complex velocity $Z = u + iv$ is introduced ($i^2 = -1$), and Eq. (14) becomes

$$ \mathcal{N}e^{i\theta_0}(Z - Z_{\text{fd}}) = \pm P_{\text{max}} \sin\phi \partial_x (\partial_x Z / |\partial_x Z|). \quad (15) $$

Defining the complex transports $T = \int_a^x Z(u) du$ and $T_{\text{fd}} = \int_a^x Z_{\text{fd}}(u) du$, Eq. (15) can be integrated. This yields

$$ e^{i\theta_0}(T - T_{\text{fd}} + C) = K \partial_{xx} T / |\partial_{xx} T|, \quad (16) $$

where $C$ is a constant and $K = \pm P_{max} \sin\phi/\mathcal{N}$. Consequently, $|T - T_{\text{fd}} + C| = |K|$ and a function $\psi(x)$ such that $T - T_{\text{fd}} - C = K \exp[i(\psi - \theta_0)]$ exists. Equation (16) simplifies to

$$ e^{i\psi} = \frac{\partial_{xx}\{T_{\text{fd}} + K \exp[i(\psi - \theta_0)]\}}{|\partial_{xx}\{T_{\text{fd}} + K \exp[i(\psi - \theta_0)]\}|}. \quad (17) $$

The amplitude $\rho_{fd} > 0$ and the direction $\alpha_{fd}$ of the derivative of the free drift are introduced; $\partial_{xx}(T_{\text{fd}}) = \partial_x u_{\text{fd}} + i\partial_x v_{\text{fd}} = \rho_{\text{fd}} \exp(i\alpha_{\text{fd}})$, and finally Eq. (17) is rewritten
$\rho_{fd} \exp[-i(\psi - \alpha_{fd})] + K \exp(i\theta_0)(-\psi'^2 + i\psi'') = A$, where $A$ is an unknown positive real and the prime notes the derivative. The imaginary part of the first member of the equality must therefore be equal to 0, whereas the real part must be positive. This leads to the differential equation

$$K \cos\theta_0\psi'' + K \sin\theta_0\psi'^2 - \rho_{fd} \sin(\psi - \alpha_{fd}) = 0 \quad (18)$$

associated with the inequality

$$K \sin\theta_0\psi'' - K \cos\theta_0\psi'^2 + \rho_{fd} \cos(\psi - \alpha_{fd}) > 0, \quad (19)$$

which adds a constraint on the coefficient $K$. This constraint is more easily understood by eliminating $\psi''$ from Eq. (19). One obtains $\rho_{fd} \cos(\psi - \alpha_{fd} - \theta_0) > K\psi'^2$. If $K$ is positive, this inequality generally is not verified. In particular, when $\rho$ is equal to 0, a case that will be investigated in section 4, solutions exist only if $K$ is negative.

The inequality $K < 0$ means that the viscosity coefficient $\eta$ becomes negative. At first glance, this property appears as a flaw in the model; it means that energy is injected in a bounded area (where $P$ reaches its maximum value) and entropy decreases, which seems in contradiction with the second principle of thermodynamics. Moreover, it will require changes in order to guarantee the stability of the numerical schemes that could be used to solve a sea ice code.

However, this negative assertion can be mitigated. In fluid mechanics, the viscosity coefficient finds its origin in the collisions between the atoms or molecules that constitute the liquid or gaz. Here, this coefficient has not a molecular origin but depends on the large-scale characteristics of the flow. Such “negative viscosity coefficients”—whose origin depends on the large-scale macroscopic properties of a system—are known to exist and have been used for a long time in geophysics and even in fluid mechanics to parameterize the turbulence (see Starr 1968).

For the “real sea ice,” a maximum compressive stress induces the formation of compressive ridges; (potential) energy is accumulated in these ridges, and more energy is available to produce work. In thermodynamics, this property is often rephrased as “the order increases.” A fully disorderly state would then correspond to a state where the ice thickness is nearly uniform, and a more orderly state would be obtained when anomalies of thickness appear. This process is not represented here since we assumed that $h$ remained constant, but the negative value of $\eta$ may convey its existence.

Last, the fact that the coefficient $\eta$ is negative means that the ratio between the macroscopic fields $\dot{\epsilon}_1 - \dot{\epsilon}_2$ and $\sigma_1 - \sigma_2$ is negative. When the compressive stress reaches

$P_{\max} + P_0$, the mechanical state of the system is defined by the two points indicated by stars at the extremities of the yield curve in Fig. 3. The flow rule is no longer specified but is determined by the system itself. In the space of the rates of deformation, a flow in any direction then becomes possible. On the contrary, when the regime is plastic shear, the flow rule specifies a deformation rate (in the half plane $\dot{\epsilon}_1 - \dot{\epsilon}_2 > 0$ when $\sigma_1 - \sigma_2$ is positive). More generally, if the yield curve is closed and convex, the flow rule assumes a deformation rate in a direction exterior to the domain bounded by the yield curve to satisfy the second law of thermodynamics (Hill 1960; Ukita and Moritz 1995).

The yield curve studied here is not convex; it is even useless to assume that it is closed. The existence of two points where $P - P_0 = P_{\max}$ is simply assumed, but we do not specify how the cone defined by the Coulomb friction law is closed (this absence is symbolized by the dotted line in Fig. 3). These assumptions are in agreement with the results of Weiss et al. (2007) based on the data of Richter-Menge et al. (2002). But they may result in a deformation rate pointing in a direction inside the half plane $\dot{\epsilon}_1 - \dot{\epsilon}_2 < 0$ when $\sigma_1 - \sigma_2$ is positive and hence in a negative coefficient $\eta$. As illustrated by the numerical results below, the derivative of the velocity is discontinuous, which induces a discontinuity of $\eta$. Its absolute value may decrease by an order of magnitude when the plastic-compressive regime is reached.

The function $\psi$ characterizes the amplitude and direction of the velocity difference $\mathbf{u} - \mathbf{u}_{fd}$. Indeed, we have $Z = Z_{fd} + iK\psi' \exp[i(\psi - \theta_0)]$, that is,

$$\mathbf{u} = \mathbf{u}_{fd} + K \frac{d\mathbf{n}}{dx} \quad \text{with} \quad K = \pm \frac{P_{\max} \sin\phi}{N} \quad \text{and} \\ \mathbf{n} = \begin{bmatrix} \cos(\psi - \theta_0) \\ \sin(\psi - \theta_0) \end{bmatrix}. \quad (20)$$

More precisely, $|\psi'| P_{\max} \sin\phi/N$ defines the norm of the vector $\mathbf{u} - \mathbf{u}_{fd}$, and $\psi - \theta_0 + \pi/2$ defines its direction.

The differential Eq. (18) is too complex to be analytically solved but standard numerical methods provide solutions with a very good accuracy. As it is a second-order equation, two arbitrary constants must be determined with the help of boundary conditions.

The domain where $P = P_{\max}$ begins far from the coast, at a distance $x = x_0$, and ends at the coast; the transition to the open sea is ensured by the plastic-shear regime previously studied—or, when the free drift is nearly uniform, by the viscous regime. In this section, only the conditions that must be applied along the unknown line $x = x_0$ separating the “plastic-shear domain” from the “plastic-compressive domain” are specified. As the
<figure><img src="image_8.png" /><figcaption>FIG. 6. (top to bottom) The free drift $\mathbf{u}_{\text{fd}}$, the corresponding sea ice velocity $\mathbf{u}$, the compressive stress $P$, and the viscosity coefficient $\eta$. The angle $\delta$ is equal to $-10^{\circ}$.</figcaption></figure>

plastic-compressive regime is still a plastic regime and the medium does not change (sea ice exist on both sides of $x = x_0$), it is reasonable to assume that the sea ice velocity and the compressive stress are continuous at $x = x_0$. This leads to three relations:

* Continuity of $P$:

$$P(x_0^+) = P_{\max} + P_0.$$

This condition, combined with the expression of $P$ in section 3a, gives a relation between $x_0$ and the

component of the sea ice velocity along the $x$ axis at $x = x_0$: $u^0$.

* Continuity of $\mathbf{u}$: A simple calculation from Eq. (12) gives $(\mathbf{u} - \mathbf{u}_{\text{fd}})(x_0^+) = \{[u^0 - \mathbf{u}_{\text{fd}}(x_0^+)]/K_x\}\mathbf{K}$, where $K_x = \langle \mathbf{K}, \mathbf{e}_x \rangle$. On the other hand, Eq. (20) yields $(\mathbf{u} - \mathbf{u}_{\text{fd}})(x_0^-)$. Combining these two expressions, one obtains

$$\begin{cases} \psi(x_0^-) - \theta_0 = \hat{K} \\ \psi'(x_0^-)P_{\max} \sin\phi = \mathcal{N}[u^0 - \mathbf{u}_{\text{fd}}(x_0^+)]\|\mathbf{K}\|/K_x \end{cases}$$

where $\hat{K}$ is the angle that defines the direction of the vector $\mathbf{K}$.
<figure><img src="image_9.png" /><figcaption>FIG. 7. The free drift is unchanged in comparison with the previous figure, but the angle δ is now equal to +10°. (top) Corresponding sea ice velocity **u**, (middle) the compressive stress **P**, and (bottom) the viscosity coefficient η.</figcaption></figure>

Last, the component of the velocity perpendicular to the coast must vanish, leading to a fourth relation:

$$u_{\text{fd}}(0) = -\psi'(0)P_{\text{max}} \sin\phi \cos[\psi(0) - \theta_0]/\mathcal{N}.$$

These four relations, associated with the differential Eq. (18) and the inequality Eq. (19), allow us to determine the function $\psi$ and the two unknowns $u^0, x_0$.

Some of the numerical solutions presented below have a remarkable supplementary property. They verify $(\mathbf{u} - \mathbf{u}_{\text{fd}})(x_0) = 0$ and $\psi'(x_0) = 0$. This means that not only the ice drift is continuous, but also the derivative of the compressive stress at $x = x_0$; $\partial_x P(x_0) = 0$.

The results shown in Figs. 6 to 10 have been obtained for $P_0 = 0$. They result from a numerical resolution based on a predictor-corrector scheme of order two and correspond to five distinct situations.

In the first three cases, the free drift increases with the distance to the coast and is oriented “eastwards” ($y < 0$). When the dilatancy angle $\delta$ is equal to $-10°$ (Fig. 6), the

plastic-shear regime acts beyond $x = 330$ km; within this bound, the plastic-compressive regime (characterized by $P$ constant) is observed as shown in Fig. 6, third panel. Note that the derivative of the compressive stress is not continuous in this case. The direction of the ice velocity is reversed close to the coast; though the free drift is eastward, the ice flows westward when $x$ ranges between 0 and 100 km. The ice slides along the coast with a velocity of about $-2 \text{ cm s}^{-1}$. The viscosity coefficient is negative in the plastic-compressive regime, as expected. Everywhere, it remains very small, 20 times smaller than the threshold value $\eta_{\text{max}}$. These characteristics remain qualitatively unchanged when the dilatancy angle is set to +10° (Fig. 7). The domain where the plastic-shear regime acts is slightly reduced (it starts from $x = 400$ km). The derivative of the compressive stress is now continuous as it can be seen in the middle panel. When $\delta = 26°$ (Fig. 8), which corresponds to a normal flow rule, the domain where the plastic-shear regime acts is
<figure><img src="image_10.png" /><figcaption>FIG. 8. As in Fig. 7, but where $\delta = 26^{\circ}$.</figcaption></figure>

considerably reduced; it only extends over 430 km. Moreover, the flow becomes very small close to the coast and remains oriented eastward.

In the last two cases, the component of the free drift parallel to the coast points “westwards” (in order to mimic a situation similar to that occurring north of the Canadian Archipelago). When the angle $\delta$ takes the value $-10^{\circ}$, a smooth ice drift is obtained; the ice velocity along the coast is about $+2 \text{ cm s}^{-1}$. When $\delta$ is equal to $+10^{\circ}$, noticeable changes can be seen. The size of the compressive domain increases from 230 to 460 km, the viscosity coefficient is multiplied by about 6 and the direction of the ice velocity is reversed along the coast. The model dependence on the dilatancy angle thus strongly depends on the direction of the free drift.

The solution ceases to exist when the components of the free drift become too large, in particular the component normal to the coast. However, in the real world, the probability that such large values exist is null; the values we have used are realistic, since their maximum is

comparable with the maxima observed in oceanic general circulation models (see Proshutinsky et al. 2011). They allow the existence of a solution.

In any case, the model behaves as a nonlinear model; the dependence of the solutions on the parameters cannot be easily guessed, as exemplified by the previous figures. The existence of a solution is not guaranteed, and, as explained in section 4, the unicity of the solution is problematic.

c. Viscous regime: $0 < P - P_0 < P_{\text{max}}$ and $\eta = \eta_{\text{max}}$

Equation (6) is still verified, and consequently the ice velocity has the same expression as in the plastic-shear regime: $\mathbf{u} = (u/\tan\delta)\mathbf{t} + v_i\mathbf{e}_y$. On the contrary, Eq. (11) becomes

$$ \mathcal{N}\mathbf{R}_0(\mathbf{u} - \mathbf{u}_{\text{fd}}) = -\partial_x P \mathbf{e}_x + \eta_{\text{max}} \partial_{xx} \mathbf{u}. \quad (21) $$

Multiplying both side of the previous equation by $\mathbf{e}_y$, and using the expression of $\mathbf{u}$, we obtain the equation
<figure><img src="image_11.png" /><figcaption>FIG. 9. (top to bottom) The free drift $\mathbf{u}_{\text{fd}}$, corresponding sea ice velocity $\mathbf{u}$, the compressive stress $P$, and the viscosity coefficient $\eta$. The angle $\delta$ is equal to $-10^{\circ}$.</figcaption></figure>

$$ \partial_{xx} u - \lambda^{-2} u = \frac{\lambda^{-2} \tan \delta}{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle} (\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{e}_y \rangle v_i - \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{\text{fd}} \rangle), \quad (22) $$

where $\lambda$ represents a characteristic length depending on the parameters on the problem:

$$ \lambda^{-2} = \frac{\mathcal{N}\langle\mathbf{e}_y, \mathbf{R}_0\mathbf{t}\rangle}{\eta_{\text{max}}\langle\mathbf{t}, \mathbf{e}_y\rangle}. $$

This length is plotted in Fig. 11 as a function of the ice thickness for $\eta_{\text{max}} = 10^{12} \text{ kg s}^{-1}$. It ranges from 1150 to 1450 km.

The solution of Eq. (22) may be decomposed into an exponential part $u_e = A \cosh(x/\lambda) + B \sinh(x/\lambda)$ and a part that depends on the free drift. The length scale $\lambda$, which is proportional to the square root of $\eta_{\text{max}}$, gives the typical scale of variations of $u_e$. If the free drift varies over scales much smaller than this typical scale (in the example of section 3a, the free drift velocity was divided by 2 over a distance of about 500 km), the spatial
<figure><img src="image_12.png" /><figcaption>FIG. 10. The free drift is unchanged in comparison with the previous figure but the angle δ is now equal to +10°. (top) The corresponding sea ice velocity **u**, (middle) the compressive stress **P**, and (bottom) the viscosity coefficient η.</figcaption></figure>

variations of **u** are driven by those of **u**<sub>fd</sub>, and the co-efficient η, estimated from [(P - P<sub>0</sub>) sinφ tanδ]/|∂<sub>x</sub>u<sub>fd</sub>|, remains smaller than η<sub>max</sub>. Consequently, the viscous regime is never reached.

This analysis confirms the result shown in section 3a, namely, that the viscous regime is possible only in areas where the free drift is nearly uniform. Numerical tests have corroborated this result. Consequently, and for simplicity, the connection between the viscous regime and the plastic-compressive regime will be considered and studied only for a uniform free drift.

## 4. Study for a uniform free drift

The case of a uniform free drift **u**<sub>fd</sub> over a large area is certainly not observed but presents two theoretical advantages: the analytical computations can be further pursued—in particular, when the plastic-compressive regime operates, they become exact—and the viscous

regime, which is difficult to reach otherwise, can be investigated.

As the plastic-shear regime is excluded when the free drift is uniform, it suffices to consider solutions corresponding to the plastic-compressive and viscous regimes. As before, the solution corresponding to the plastic-compressive regime is valid in the domain defined by 0 ≤ x ≤ x<sub>0</sub>, where x<sub>0</sub> is an unknown abscissa. It is obtained from the resolution of Eq. (18) with ρ<sub>fd</sub> = 0. It yields ψ = ψ<sub>0</sub> + [ln(x + c)/tanθ<sub>0</sub>], where ψ<sub>0</sub> and c are two arbitrary constants. The sea ice velocity is thus equal to

$$ \mathbf{u} = \mathbf{u}_{\text{fd}} + \frac{P_{\text{max}} \sin\phi}{\mathcal{N} \tan\theta_0} \frac{\mathbf{n}}{x+c} \quad \text{with} \quad \mathbf{n} = -\sin\psi \mathbf{e}_x + \cos\psi \mathbf{e}_y, \qquad (23) $$

and the compressive stress is P = P<sub>max</sub> + P<sub>0</sub>. Note that the sign — and the angle θ<sub>0</sub> in Eq. (20) have been absorbed in the arbitrary constant ψ<sub>0</sub>. With this
<figure><img src="image_13.png" /><figcaption>FIG. 11. Variations (km) of $\lambda$ as a function of the ice thickness $h$. The top panel is for $\epsilon = 1$, and the bottom panel for $\epsilon = -1$.</figcaption></figure>

solution, the “viscosity” $\eta$ is equal to $-P_{\max} \sin\phi/\sqrt{(\partial_x u)^2 + (\partial_x v)^2} = -\mathcal{N} \tan\theta_0 \sin\theta_0 (x+c)^2$ and is always negative (see section 3b).

The solution corresponding to the viscous regime is valid in the domain defined by $x_0 \le x \le L$. Using the results of section 3c and the fact that a particular solution of Eq. (22) is a constant, one obtains

$$ \mathbf{u} = \left[ A \cosh\left(\frac{x-L}{\lambda}\right) + B \sinh\left(\frac{x-L}{\lambda}\right) + \frac{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{fd} \rangle - v_i \cos\theta_0}{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle} \right] \mathbf{t} + v_i \mathbf{e}_y. \quad (24) $$

The difference $x - L$ has been introduced to facilitate the ulterior computations. At $x = L$, the compressive stress vanishes; in order to keep the model consistent, the expression of the coefficient $\eta$ obtained on the yield curve ($\eta = P \sin\phi \tan\delta/\partial_x u$) must remain larger than $\eta_{\max}$. This implies $\partial_x u(x = L) = 0$ and consequently $B = 0$. Using this result, fastidious algebraic transformations (see appendix A) eventually lead to the expression

$$ \mathbf{u} - \mathbf{u}_{fd} = A \cosh\left(\frac{x-L}{\lambda}\right) \mathbf{t} + w \mathbf{R}_0^{-1} \mathbf{e}_x. \quad (25) $$

Here, $w = \langle \mathbf{t}^{\perp}, \mathbf{u}_{fd} - v_i \mathbf{e}_y \rangle / \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle$ depends on the arbitrary constant $v_i$ and will be considered as an arbitrary

new constant. The previous expression allows us to compute the pressure from Eq. (21). We first obtain

$$ -\partial_x P = \mathcal{N} w - \frac{\mathcal{N} \sin\theta_0}{\langle \mathbf{t}, \mathbf{e}_y \rangle} \left[ A \cosh\left(\frac{x-L}{\lambda}\right) \right], $$

and after integration, considering the fact that $P(L) = 0$, we obtain

$$ P(x) = \mathcal{N}\lambda \left[ w \frac{L-x}{\lambda} + \frac{A \sin\theta_0}{\langle \mathbf{t}, \mathbf{e}_y \rangle} \sinh\left(\frac{x-L}{\lambda}\right) \right]. \quad (26) $$

Equations (25) and (26) are valid as long as the compressive stress $P$ does not exceed the threshold value $P_{\max}$.

The solutions given by Eqs. (23) and (25) join at $x = x_0$. The conditions given in section 3b (continuity of $P$ and $\mathbf{u}$), which rest upon physical principles, are still valid as well as the one given at $x = 0$ (vanishing of $u$). However, they do not permit to determine the five arbitrary constants that have been introduced, namely, $A$, $w$, $x_0$, $\psi_0$, and $c$. A supplementary condition must be added.

When the system is close to the extrema point indicated by a star in Fig. 4, the viscosity on the yield curve is given by $\eta = P_{\max} \sin\phi \tan\delta/\partial_x u$; it must tend to $\eta_{\max}$.
when $x \to x_0$. The derivative of the velocity $\partial_x u$ must therefore verify

$$ \partial_x u(x = x_0^+) = P_{\max} \sin\phi \tan\delta/\eta_{\max}, $$

which provides the searched supplementary condition.

Appendix B gives technical details about the analytical method that allows us to compute the coefficients $A$, $w$, $x_0$, $\psi_0$, and $c$ from this last condition, the continuity of $P$ and $u$ at $x = x_0$, and the vanishing of $u$ on the coast $[u(0) = 0]$.

In the next two subsections, numerical solutions obtained from Eqs. (23), (25), and (26), plus the boundary and matching conditions, are shown for two different geometries (note that $P_0$ has been set to 0). In the first one, we consider a uniform free drift in a geometry similar to that of section 2c, except that the size of the domain covered with ice has been reduced to 500 km. In the second one, the geometry is changed; the sea ice is compressed between two walls, a situation that might occur if it occupies the whole width of a strait.

### a. Case of an open domain

An eastward uniform-free drift of $7.2 \text{ cm s}^{-1}$ is applied over a domain that extends 500 km from the coast (top panel in Fig. 12). The solution predicted by the model is first shown for a dilatancy angle equal to $-10^\circ$ (the coefficient $N$ and matrix $R_0$ remain unchanged).

It strongly differs from the one obtained when the plastic-shear and plastic-compressive regimes acted. First, the plastic-compressive regime is never reached (bottom panel). Second, the ice drift is much smaller than the free drift, never exceeding $1 \text{ cm s}^{-1}$ (middle panel). Last, the viscosity is equal (by hypothesis) to $\eta_{\max}$ and thus is considerably larger than the viscosity previously obtained for an eastward free drift (Figs. 6 to 8).

The response of the model significantly changes when $\delta$ is equal to $+10^\circ$ (Fig. 13). The plastic-compressive regime is active between $x = 0$ and $x = 170 \text{ km}$, and the solution is not realistic. Indeed, the maximum of the velocity is reached along the coast where it takes an unrealistic value of about $53 \text{ cm s}^{-1}$. Interestingly, other solutions may be built, which present the same characteristics, with a still larger velocity along the coast. The nonunicity of the solution is made explicit in the next subsection.

Even though high velocities of sea ice have been reported close to the coasts, it seems difficult to imagine that such a solution might be observed in the Arctic. The high velocity results from the combination of an unrealistic forcing field (the velocity is uniform over 500 km) and a strongly nonlinear behavior of the model when the plastic-compressive regime is reached. If the

domain is reduced to 230 km, the plastic-compressive regime is no longer observed; only the viscous regime acts, and the solution becomes ordinary again.

### b. Ice flow through a strait

To clarify the meaning of the solutions found for the plastic-compressive regime, a case where the latter acts alone is considered now. We suppose that sea ice have formed in a strait of width $L = 200 \text{ km}$ and occupies it completely. Moreover, we assume that the compressive stress has everywhere reached the maximum value $P_{\max}$. It is not clear if such a case could occur in practice—uniform, thick ice should occupy the strait in its entirety—but it has a great theoretical interest; indeed, the system is thus governed only by Eq. (20).

For a uniform drift ($u_{fd}$ independent of $x$), an analytical solution is given by Eq. (23) as explained at the beginning of this section. The boundary conditions at $x = 0$ and $x = L$ add two constraints, which allow us to determine the coefficients $c$ and $\psi_0$:

$$ \begin{aligned} \mathbf{u}_{\text{fd}} \cdot \mathbf{e}_x &= -(K/c) \sin\{\psi_0 + [\ln(c)/\tan\theta_0]\} \\ &= -[K/(c+L)] \sin\{\psi_0 + [\ln(c+L)/\tan\theta_0]\}, \end{aligned} \quad (27) $$

with $K = P_{\max} \sin\phi/(N \tan\theta_0)$.

Figure 14 shows the prescribed free drift inside the strait (top panel). A solution given by Eq. (23) and that satisfies Eq. (27) is shown in the bottom panel. The ice velocity peaks at the western coast (on the left) where it reaches $21.18 \text{ cm s}^{-1}$. Velocities larger than $10 \text{ cm s}^{-1}$ appear close to the western coast, at a distance shorter than 30 km. They become negligible beyond 100 km. The viscosity coefficient $\eta$ is everywhere negative, but its absolute value remains very small; the maximum is about 100 times smaller than $\eta_{\max}$ and the mean value 1000 times smaller.

This solution is not unique; other solutions given by Eq. (23) and that satisfy Eq. (27) may be found. This strange behavior is caused by the sine function in Eq. (27). It makes possible the determination of different couples ($\psi_0$, $c$); they correspond to successive oscillations of the sine function and satisfy the identities of Eq. (27). Figure 15 illustrates this property; it shows a second solution corresponding to another couple of values $\psi_0$ and $c$ (the free drift is the same as in Fig. 14). Other solutions have still been found; they showed a higher velocity along the coast and thus were less realistic (not shown).

This nonunicity property is made explicit in the case of a vanishing free drift along $e_x$. Equation (27) admits as solutions $\psi_0 + [\ln(c)/\tan\theta_0] = p\pi$ and $\psi_0 + [\ln(c+L)/\tan\theta_0] = q\pi$, where $p$ and $q$ are integers.
<figure><img src="image_14.png" /><figcaption>FIG. 12. (top) The free drift is now uniform, and the angle $\delta$ is equal to -10°. The (middle) sea ice velocity is given by the viscous regime since the (bottom) compressive stress remains smaller than $P_{max}$. The viscosity coefficient $\eta$ is equal to $\eta_{max}$.</figcaption></figure>

Thus, $\ln(1 + L/c) = (q - p)\pi \tan\theta_0$ and $c = L/\{-1 + \exp[(q - p)\pi \tan\theta_0]\}$. An infinite family of solutions is therefore obtained, indexed by the numbers $q - p$.

**5. Conclusions**

The model studied in this paper is a viscous-plastic model, close to the models used for climatological studies. It is based on a classical assumption of isotropy and assumes that the plastic regime may be characterized by a Mohr–Coulomb yield curve. The flow rule, which prescribes how plastic deformation takes place, has a physical ground, since it may be obtained from the study of granular materials (Tremblay and Mysak 1997). It contains a parameter $\delta$, the angle of dilatancy, which can be easily modified. In particular, a choice of this parameter corresponding to the famous normal flow rule may be used, even though it is not realistic in such a model.

To study the properties of this model, an analytical approach has been chosen. Indeed, because of the complexity of sea ice models, the efficiency and

precision of the numerical schemes do not seem to be precisely known. It is therefore difficult to discriminate between the inaccuracies due to the numerical method and those due to the physical parameterizations. Consequently, to complement the numerical studies, it seems interesting to look for and study exact solutions. A simple configuration has thus been defined in which the physical variables depend only on the distance to the coast. It allows us to pursue analytical computations far beyond what is usually done.

In this model, three distinct regimes of sea ice drift have been defined: the plastic-shear regime, the viscous regime, and a particular case of the plastic regime that we named the plastic-compressive regime:

* In the plastic-shear regime, the compressive stress $P$ remains smaller than the maximum value $P_{max} + P_0$ ($P_0$ is a tensile stress), and the ratio between the stress $\sigma_1 - \sigma_2$ and the strain rate $\dot{\epsilon}_1 - \dot{\epsilon}_2$, $\eta$, remains smaller than the maximum value $\eta_{max}$. As shown by the analytical study, this regime occurs as soon as the
<figure><img src="image_15.png" /><figcaption>FIG. 13. The free drift is unchanged in comparison to the previous figure, but the angle $δ$ is now equal to +10°. (top) The corresponding sea ice velocity **u**, (middle) the compressive stress **P**, and (bottom) the viscosity coefficient $η$.</figcaption></figure>

gradient of the forcing over the considered domain is sufficiently high; in a realistic situation, this condition is almost always fulfilled.

* In the viscous regime, $η$ is equal to $η_{max}$. This regime can occur only if the gradient of the forcing field is very weak (to give an order of magnitude, variations of the free drift not exceeding 0.2 cm s⁻¹ over a distance of 100 km).

* In the plastic-compressive regime, the compressive stress $P$ is equal to $P = P_{max} + P_0$. This regime is mainly active close to the coast, up to a distance depending on the forcing field and dilatancy angle.

This model behaves in such a way that it is useless to define a curve that links the extremities of the straight segments that define the plastic behavior. Indeed, when the compressive stress increases, the state predicted by the model gets closer to the plastic state, and it reaches the maximum value $P = P_{max} + P_0$ tangentially

to the segments defining the plastic behavior. The yield curve thus is not closed, and the whole pattern fits the analysis by Weiss et al. (2007) of the experimental data of Richter-Menge et al. (2002).

The flow rule is defined by the dilatancy angle $δ$. In the model proposed by Tremblay and Mysak (1997), $δ$ emerges naturally from considering the small-scale properties of the sea ice; it represents the angle made by the planes associated with small- and large-scale slidings. It thus depends on the properties of the medium, for example, on the size of the ice floes, on the ice concentration, and so on. It is thus highly probable that a relation is missing between $δ$ and some of the variables that characterize sea ice. This relation would constitute the state equation of the system. Note that, in the traditional Hibler model, a state equation has been introduced; it usually takes the form $P = P^*h \exp[C(1 - A)]$ and allows us to determine the size of the yield curve
<figure><img src="image_16.png" /><figcaption>FIG. 14. (top) The free drift is uniform with an amplitude of 2.4 cm s⁻¹. (bottom) The corresponding ice drift.</figcaption></figure>

as a function of ice thickness and ice concentration. An equation that would play a similar role is absent here.

The fact that we have not introduced a state equation defining $\delta$ is probably acceptable, since we assumed that $h$ is constant. The experiments have been done for $\delta = -10°$ (convergent ice drift), $\delta = 10°$ (divergent ice drift), and $\delta = 26°$, this last and unrealistic value corresponding to a normal flow rule. The response of the model depends on $\delta$ and the forcing field (here the free drift). For example, we observed that the size of the domain where the plastic-shear (plastic compressive) regime is active decreases (increases) when the angle $\delta$ increases. These variations may be moderate (cf. Figs. 6 and 7) or large (cf. Figs. 9 and 10). In the latter case, the direction of the ice drift even reverses close to the coast. These results suggest that the determination of $\delta$ (e.g., as a function of $h$ and $A$), or more generally of the flow rule, may have a significant impact on the ice dynamics predicted by a sea ice model.

Our results also suggest that the problem of the boundary conditions in the numerical models should be considered in a more detailed way. This study proves that, in a model with a Mohr–Coulomb yield curve, the classical no-slip condition too strongly constrains

the sea ice dynamics and consequently may prevent the existence of a solution in the plastic regime. On the contrary, if a sliding of the ice along the coasts is allowed, a solution can be computed.

It is not clear if this result is still valid when a more complex yield curve or a different model is used. However, it is known that sea ice models are not always numerically stable (see, e.g., Lemieux et al. 2010). This could be the consequence of a too constraining boundary condition. The velocity is assumed to vanish on the coasts in most of the sea ice models, but there is no clear physical reason for such an assumption. In solid mechanics, the sliding of a solid over another one may be allowed, and, adding to the vanishing of the component of the velocity perpendicular to the interface between the two solids, supplementary boundary conditions may apply on the components of the stress tensor itself (continuity of stress at the interface).

As reminded above, a specificity of this model is the existence of two extremal points on the curve of plasticity, where the compressive stress is maximal. These points, which define the plastic-compressive regime, lead to a dynamics of sea ice that dramatically differs

<figure><img src="image_17.png" /><figcaption>FIG. 15. Another unrealistic solution for the same free drift as in the previous figure. There is no unicity of the solutions.</figcaption></figure>
from that obtained in the other regimes. Indeed, when the plastic-compressive regime is reached, the ratio between the stress $\sigma_1 - \sigma_2$ and the strain rate $\dot{\epsilon}_1 - \dot{\epsilon}_2$ may become negative. This formally leads us to introduce a negative viscosity coefficient $\eta$.

Associated with the plastic-compressive regime, a second difficulty may happen; there may exist several (an infinity of) solutions that verify the model equations and the boundary conditions for a given forcing. This result is compatible with the nonlinear character of the equations defining this regime. It means that two different states could be reached for the same final stationary forcing, depending on the initial conditions that are prescribed or the history of the system.

These oddities have been previously discussed and reveal not only the richness of the mathematics and physics hidden in the model, but also the difficulties it could generate if a numerical integration is performed. Indeed, numerical schemes have not been developed to tackle these unusual (perhaps even unphysical) behaviors. They are probably caused by the existence of two salient points at the extremity of the segments defining the yield curve. Thus, a yield curve that would be convex in the vicinity of these points would certainly allow the elimination of such a phenomenon. The problem of the boundary conditions is probably more difficult to solve. A solution is to add dissipative terms (a Laplacian is sufficient as illustrated by the viscous case). A more sophisticated solution would be to modify the definition of $\eta$ as suggested by Gray (1999); the expression $\eta = S/(D + D_{\min})$ instead of $\eta = S/D$ would prevent the simplification that occurs in section 3a and would allow us to keep more derivatives in order to satisfy more constraining boundary conditions. It is, however, not clear whether such modifications have a physical basis.

Last, we want to emphasize that the results have been derived from a realistic model of the sea ice rheology. The setup is very simple (the domain has a straight boundary, the forcing field is uniform in the direction parallel to the boundary, and the ice thickness is constant) and limits the range of applicability of the results, but it allowed us to solve the model equations without doing any approximations that could simplify the physical mechanisms driving the rheology of sea ice. The obtained solutions are computed thanks to analytical methods, the use of numerical methods being restricted to the resolution of only one differential equation. Consequently, these solutions are not altered by an inappropriate numerical scheme or numerical errors, which are always difficult to find and analyze. They are just the consequence of the physical principles that defined the rheology of the model.

**Acknowledgments.** We thank the anonymous reviewers for their helpful comments. This research was supported in part by a grant from the French program LEFE/EVE, project TIAMMO. It has also received support from the European Union 6th Framework Programme under Grant Agreement 308290 NACLIM.

# APPENDIX A

## Isotropy and System of Eqs. (3)

The stress tensor $\sigma$ is symmetric and thus can be diagonalized using a rotation matrix; we have

$$ \begin{bmatrix} \sigma_{xx} & \sigma_{xy} \\ \sigma_{xy} & \sigma_{yy} \end{bmatrix} = \begin{bmatrix} \cos\psi & \sin\psi \\ -\sin\psi & \cos\psi \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix} \times \begin{bmatrix} \cos\psi & -\sin\psi \\ \sin\psi & \cos\psi \end{bmatrix}. $$

This matrix relation leads to the three equations:

$$ \begin{cases} \sigma_{xx} = (\sigma_1 + \sigma_2)/2 + \cos(2\psi)(\sigma_1 - \sigma_2)/2 \\ \sigma_{xy} = \sigma_{yx} = \sin(2\psi)(\sigma_1 - \sigma_2)/2 \\ \sigma_{yy} = (\sigma_1 + \sigma_2)/2 - \cos(2\psi)(\sigma_1 - \sigma_2)/2 \end{cases}. $$

The same relations are verified by the strain rate tensor because of the isotropy hypothesis. We use them to find the expression of $\cos(2\psi)$ and $\sin(2\psi)$. Since $\dot{\epsilon}_1$ and $\dot{\epsilon}_2$ are the eigenvalues of $\dot{\epsilon}$, we first have $\dot{\epsilon}_1 = [\text{Tr}(\dot{\epsilon}) + D]/2$ and $\dot{\epsilon}_2 = [\text{Tr}(\dot{\epsilon}) - D]/2$, where $D = \sqrt{[\text{Tr}(\dot{\epsilon})]^2 - 4(\dot{\epsilon}_{xx}\dot{\epsilon}_{yy} - \dot{\epsilon}_{xy}\dot{\epsilon}_{xy})}$. Using the expressions defining the component of the strain rate tensor, we obtain $D = \dot{\epsilon}_1 - \dot{\epsilon}_2 = +\sqrt{(\partial_x u - \partial_y v)^2 + (\partial_x v + \partial_y u)^2}$. This leads to $\sin(2\psi) = 2\dot{\epsilon}_{xy}/D$ and consequently to

$$ \sigma_{xy} = \eta(\partial_x u + \partial_y v), $$

with $\eta = S/D$ and $S = (\sigma_1 - \sigma_2)/2$.

Similarly, by subtracting the equation giving $\dot{\epsilon}_{yy}$ from that giving $\dot{\epsilon}_{xx}$, we have $\cos(2\psi) = (\dot{\epsilon}_{xx} - \dot{\epsilon}_{yy})/D$. This leads to

$$ \sigma_{xx} = -P + \eta(\partial_x u - \partial_y v) $$

and

$$ \sigma_{yy} = -P - \eta(\partial_x u - \partial_y v), $$

where $P = -(\sigma_1 + \sigma_2)/2$. The system of Eqs. (3) is thus obtained.
## APPENDIX B

### Details about Some Computations

This appendix gives the details of the calculations that were omitted in the text for the sake of simplicity. They are now presented for the reader who is interested in these technical points. Each computation is characterized by the labels of the equations to which it is related. The comments are reduced to the minimum:

* From Eq. (12) to Eq. (13). With Eq. (11), we have

$$ \mathbf{u} = u\mathbf{e}_x + \left[ u \frac{\epsilon\sqrt{1 - \tan^2\delta}}{\tan\delta} + v_i \right] \mathbf{e}_y. $$

Using the definition of $\mathbf{t}$ and the fact that $u$ is equal to $\tan\delta(-v_i\langle\mathbf{e}_y, \mathbf{K}^\perp\rangle + \langle\mathbf{u}_{fd}, \mathbf{K}^\perp\rangle)/\langle\mathbf{t}, \mathbf{K}^\perp\rangle$ (see section 3a), we also have

$$ \mathbf{u} = \frac{u}{\tan\delta} \mathbf{t} + v_i \mathbf{e}_y = \frac{-v_i \langle \mathbf{e}_y, \mathbf{K}^\perp \rangle + \langle \mathbf{u}_{fd}, \mathbf{K}^\perp \rangle}{\langle \mathbf{t}, \mathbf{K}^\perp \rangle} \mathbf{t} + v_i \mathbf{e}_y. $$

The unknown $u^0$ has been introduced in the text; it is the component of the velocity $\mathbf{u}$ along $\mathbf{e}_x$ at the unknown abscissa $x = x_0$. Using this definition, we have

$$ u^0 = \langle \mathbf{u}(x=x_0), \mathbf{e}_x \rangle = \frac{-v_i \langle \mathbf{e}_y, \mathbf{K}^\perp \rangle + \langle \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle}{\langle \mathbf{t}, \mathbf{K}^\perp \rangle} \langle \mathbf{t}, \mathbf{e}_x \rangle, $$

and consequently

$$ v_i \langle \mathbf{e}_y, \mathbf{K}^\perp \rangle = - \frac{\langle \mathbf{t}, \mathbf{K}^\perp \rangle}{\langle \mathbf{t}, \mathbf{e}_x \rangle} u^0 + \langle \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle. $$

Thus,

$$ \mathbf{u} = \frac{\langle \mathbf{u}_{fd} - \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle}{\langle \mathbf{t}, \mathbf{K}^\perp \rangle} \mathbf{t} + \frac{\langle \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle}{\langle \mathbf{e}_y, \mathbf{K}^\perp \rangle} \mathbf{e}_y + u^0 \mathbf{A}, $$

with

$$ \mathbf{A} = \frac{1}{\langle \mathbf{t}, \mathbf{e}_x \rangle} \left( \mathbf{t} - \frac{\langle \mathbf{t}, \mathbf{K}^\perp \rangle}{\langle \mathbf{e}_y, \mathbf{K}^\perp \rangle} \mathbf{e}_y \right). $$

The vector $\mathbf{A}$ is easy to compute by introducing the component of $\mathbf{t}$ and $\mathbf{K}^\perp$ in the orthonormal basis ($\mathbf{e}_x, \mathbf{e}_y$). It is equal to $\mathbf{K}/\langle\mathbf{e}_y, \mathbf{K}^\perp\rangle = \mathbf{K}/\langle\mathbf{e}_x, \mathbf{K}\rangle$. Equation (13) is thus shown.

* Equation (14) needs the computation of $\langle \mathbf{t}^\perp, \mathbf{u} - \mathbf{u}_{fd} \rangle$. First, using Eq. (13),

$$ \langle \mathbf{t}^\perp, \mathbf{u} \rangle = \frac{\langle \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle \langle \mathbf{t}^\perp, \mathbf{e}_y \rangle}{\langle \mathbf{e}_x, \mathbf{K} \rangle} + u^0 \frac{\langle \mathbf{t}^\perp, \mathbf{K} \rangle}{\langle \mathbf{e}_x, \mathbf{K} \rangle} $$

(we recall that $\langle \mathbf{e}_y, \mathbf{K}^\perp \rangle = \langle \mathbf{e}_x, \mathbf{K} \rangle$). We remark that

$$ \langle \mathbf{u}_{fd}^0, \mathbf{K}^\perp \rangle \langle \mathbf{t}^\perp, \mathbf{e}_y \rangle = \langle \mathbf{u}_{fd}^0, \mathbf{t}^\perp \rangle \langle \mathbf{e}_x, \mathbf{K} \rangle - \langle \mathbf{K}, \mathbf{t}^\perp \rangle \langle \mathbf{u}_{fd}^0, \mathbf{e}_x \rangle. $$

Consequently,

$$ \partial_x P = \mathcal{N} \left( \frac{\langle \mathbf{u}_{fd}^0 - \mathbf{u}_{fd}, \mathbf{t}^\perp \rangle}{\langle \mathbf{K}, \mathbf{t}^\perp \rangle} + \frac{u^0 - \langle \mathbf{u}_{fd}^0, \mathbf{e}_x \rangle}{\langle \mathbf{e}_x, \mathbf{K} \rangle} \right). $$

* From Eq. (25) to Eq. (26). We have to compute

$$ \frac{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{fd} \rangle - v_i \cos\theta_0}{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle} \mathbf{t} + v_i \mathbf{e}_y - \mathbf{u}_{fd}. $$

On the one hand,

$$ \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{fd} \rangle \mathbf{t} - \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle \mathbf{u}_{fd} = (u^x \sin\theta_0 + u^y \cos\theta_0)\mathbf{t} - (t^x \sin\theta_0 + t^y \cos\theta_0)\mathbf{u}_{fd}, $$

where $(u^x, u^y)$ and $(t^x, t^y)$ note the components of $\mathbf{u}_{fd}$ and $\mathbf{t}$ in the basis $(\mathbf{e}_x, \mathbf{e}_y)$. The previous equality simplifies and yields

$$ \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{fd} \rangle \mathbf{t} - \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle \mathbf{u}_{fd} = \langle \mathbf{u}_{fd}, \mathbf{t}^\perp \rangle \mathbf{R}_0^{-1} \mathbf{e}_x $$

(note that the relation $\cos\theta_0\mathbf{e}_x - \sin\theta_0\mathbf{e}_y = \mathbf{R}_0^{-1}\mathbf{e}_x$ has been used).

On the other hand,

$$ \cos\theta_0 \mathbf{t} - \langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle \mathbf{e}_y = t^x \mathbf{R}_0^{-1} \mathbf{e}_x = \langle \mathbf{t}^\perp, \mathbf{e}_y \rangle \mathbf{R}_0^{-1} \mathbf{e}_x. $$

Finally, we obtain

$$ \frac{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{u}_{fd} \rangle - v_i \cos\theta_0}{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle} \mathbf{t} + v_i \mathbf{e}_y - \mathbf{u}_{fd} = \frac{\langle \mathbf{u}_{fd} - v_i \mathbf{e}_y, \mathbf{t}^\perp \rangle}{\langle \mathbf{e}_y, \mathbf{R}_0 \mathbf{t} \rangle} \mathbf{R}_0^{-1} \mathbf{e}_x $$

and relation (26) follows.

## APPENDIX C

### Computation of A, w, $\psi_0$, c, and $x_0$ (Section 4)

The condition on the viscosity coefficient $\lim_{x \to x_0} \eta_{\max} \partial_x u = P_{\max} \sin\phi \tan\delta$ leads to the relation

$$ A \eta_{\max} \sinh[(x_0 - L)/\lambda] = \lambda P_{\max} \sin\phi, $$
which defines A as a function of $x_0$.

The condition $P(x = x_0) = P_{max} + P_0$ yields the relation

$$w(x_0 - L) = \frac{P_{max}}{N} \left( \frac{\sin\theta_0 \sin\phi}{\langle \mathbf{R}_0 \cdot \mathbf{t}, \mathbf{e}_y \rangle} - 1 \right) - \frac{P_0}{N},$$

where the expressions of $\lambda^2$ and A have been used to simplify it. This defines w as a function of $x_0$.

To remark, we have

$$\partial_x \mathbf{u} = \frac{P_{max} \sin\phi}{\eta_{max}} \frac{\sinh[(x - L)/\lambda]}{\sinh[(x_0 - L)/\lambda]} \mathbf{t}.$$

The consistency of the model is thus ensured since the viscosity given by the plastic-shear regime $\eta = P_{max} \sin\phi \tan\delta/\partial_x u$ is equal to $\eta_{max}\{\sinh[(x_0 - L)/\lambda]/\sinh[(x - L)/\lambda]\}$ and is everywhere larger than $\eta_{max}$.

The conditions on the ice velocity at $x = 0$ and $x = x_0$ add three new relations that allow us to determine the remaining constants. In the interval $0 \le x \le x_0$, the ice drift is given by the relation

$$\mathbf{u} - \mathbf{u}_{\text{fd}} = \frac{P_{\text{max}} \sin\phi}{N \tan\theta_0} \frac{1}{x+c} (-\sin\psi \mathbf{e}_x + \cos\psi \mathbf{e}_y),$$

with $\psi = \psi_0 + [\log(x + c)]/\tan\theta_0$. At $x = 0$, the component of the velocity along $\mathbf{e}_x$ must vanish. Consequently,

$$\frac{1}{c} \sin \left( \psi_0 + \frac{\log c}{\tan \theta_0} \right) = \frac{N \tan \theta_0 u_{\text{fd}}(0)}{P_{\text{max}} \sin \phi}. \quad (28)$$

We define the vector $\mathbf{N} = N[\cos(\hat{N})\mathbf{e}_x + \sin(\hat{N})\mathbf{e}_y]$ by the following relation:

$$\mathbf{u}(x_0^+) - \mathbf{u}_{\text{fd}}(x_0^+) = \mathbf{N}.$$

The notation $x_0^+$ means that we consider a limit in the interval $x_0 \le x \le L$. Consequently, the expression of $\mathbf{u}$ is given by Eq. (26) with $x = x_0$. As A and w depend only on $x_0$, the two functions N and $\hat{N}$ also depends only on $x_0$.

At $x = x_0$, the components of the velocity must be continuous. The two following relations are thus obtained:

$$\psi_0 + \frac{\log(x_0 + c)}{\tan\theta_0} + \frac{\pi}{2} = \hat{N}(x_0), \quad (29)$$

and

$$\frac{P_m \sin\phi}{N \tan\theta_0(x_0 + c)} = N(x_0), \quad (30)$$

where we have recalled that N and $\hat{N}$ depends on $x_0$. Combining Eqs. (28) and (29), one obtains

$$\frac{1}{c} \cos \left[ \hat{N}(x_0) - \frac{1}{\tan\theta_0} \log \left( 1 + \frac{x_0}{c} \right) \right] = - \frac{N \tan\theta_0 u_{\text{fd}}(0)}{P_m \sin\phi}. \quad (31)$$

The unknown c is easily expressed as a function of $x_0$ from Eq. (30). This expression is introduced into Eq. (31) and the resulting equation allows us to find the unknown $x_0$.

## REFERENCES

Coon, M. D., G. A. Maykut, R. S. Pritchard, D. A. Rothrock, and A. S. Thorndike, 1974: Modeling the ice pack as an elastic-plastic material. *AIDJEX Bull.*, **24**, 1–106. [Available online at http://psc.apl.washington.edu/nonwpProjects/aidjex/files/AIDJEX-24.pdf.]

—, R. Kwok, G. Levy, M. Pruis, H. Schreyer, and D. Sulsky, 2007: Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate. *J. Geophys. Res.*, **112**, C11S90, doi:10.1029/2005JC003393.

Gray, J. M. N. T., 1999: Loss of hyperbolicity and ill-posedness of the viscous-plastic sea ice rheology in uniaxial divergent flow. *J. Phys. Oceanogr.*, **29**, 2920–2929, doi:10.1175/1520-0485(1999)029<2920:LOHAIP>2.0.CO;2.

Heil, P., and W. D. Hibler III, 2002: Modeling the high-frequency component of Arctic sea ice drift and deformation. *J. Phys. Oceanogr.*, **32**, 3039–3057, doi:10.1175/1520-0485(2002)032<3039:MTHFCO>2.0.CO;2.

Hibler, W. D., III, 1979: A dynamic thermodynamic sea ice model. *J. Phys. Oceanogr.*, **9**, 815–846, doi:10.1175/1520-0485(1979)009<0815:ADTSIM>2.0.CO;2.

—, and E. M. Schulson, 2000: On modeling the anisotropic failure and flow of flawed sea ice. *J. Geophys. Res.*, **105**, 17105–17120, doi:10.1029/2000JC900045.

Hill, R., 1960: *The Mathematical Theory of Plasticity*. Oxford University Press, 355 pp.

Hunke, E. C., and J. K. Dukowicz, 1997: An elastic-viscous-plastic model for sea ice dynamics. *J. Phys. Oceanogr.*, **27**, 1849–1867, doi:10.1175/1520-0485(1997)027<1849:AEVPMF>2.0.CO;2.

—, and —, 2002: The elastic-viscous-plastic sea ice dynamics model in general orthogonal curvilinear coordinates on a sphere— Incorporation of metric terms. *Mon. Wea. Rev.*, **130**, 1848–1865, doi:10.1175/1520-0493(2002)130<1848:TEVPSI>2.0.CO;2.

Hutchings, J. K., P. Heil, and W. D. Hibler III, 2005: Modeling linear kinematic features in sea ice. *Mon. Wea. Rev.*, **133**, 3481–3497, doi:10.1175/MWR3045.1.

Lemieux, J. F., B. Tremblay, S. Thomas, J. Sedláček, and L. A. Mysak, 2008: Using the preconditioned generalized minimum residual (GMRES) method to solve the sea ice momentum equation. *J. Geophys. Res.*, **113**, C10004, doi:10.1029/2007JC004680.

—, —, J. Sedláček, P. Tupper, S. Thomas, D. Huard, and J.-P. Auclair, 2010: Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton-Krylov method. *J. Comput. Phys.*, **229**, 2840–2852, doi:10.1016/j.jcp.2009.12.011.

Leppäranta, M., 2005: *The Drift of Sea Ice*. Springer, 280 pp.
Marsan, D., H. Stern, R. Lindsay, and J. Weiss, 2004: Scale dependence and localization of the deformation of Arctic sea ice. *Phys. Rev. Lett.*, **93**, 178501, doi:10.1103/PhysRevLett.93.178501.

McPhee, M., 1975: Ice-ocean momentum transfer for the AIDJEX ice model. *AIDJEX Bull.*, **29**, 93–112. [Available online at http://psc apl.washington.edu/nonwpprojects/aidjex/files/AIDJEX-29.pdf.]

Pritchard, R. S., 2001: Long-term sea ice dynamics simulations using an elastic-plastic constitutive law. *J. Geophys. Res.*, **106**, 31333–31343, doi:10.1029/2000JC000638.

Proshutinsky, A., and Coauthors, 2011: Recent advances in Arctic Ocean studies employing models from the Arctic Ocean Model Intercomparison Project. *Oceanography*, **24**, 102–113, doi:10.5670/oceanog.2011.61.

Rampal, P., J. Weiss, D. Marsan, R. Lindsay, and H. Stern, 2008: Scaling properties of sea ice deformation from buoy dispersion analysis. *J. Geophys. Res.*, **113**, C03002, doi:10.1029/2007JC004143.

Richter-Menge, J., S. L. McNutt, J. E. Overland, and R. Kwok, 2002: Relating Arctic pack ice stress and deformation under winter conditions. *J. Geophys. Res.*, **107**, 8040, doi:10.1029/2000JC000477.

Salençon, J., 2002: De l'Élasto-plasticité au calcul de la rupture. Éditions de l'École polytechnique, Palaiseau, 262 pp.

Schreyer, H. L., D. L. Sulsky, L. B. Munday, M. D. Coon, and R. Kwok, 2006: Elastic-decohesive constitutive model for sea ice. *J. Geophys. Res.*, **111**, C11S26, doi:10.1029/2005JC003334.

Schulson, E. M., 2001: Brittle failure of ice. *Eng. Fract. Mech.*, **68**, 1839–1887, doi:10.1016/S0013-7944(01)00037-6.

—, and O. Nickolayev, 1995: The failure of columnar saline ice under biaxial compression: Failure envelopes and the brittle-to-ductile transition. *J. Geophys. Res.*, **100**, 22383–22400, doi:10.1029/95JB02513.

Sedlacek, J., J.-F. Lemieux, L. A. Mysak, L.-B. Tremblay, and D. M. Holland, 2007: The granular sea ice model in spherical coordinates and its application to a global climate model. *J. Climate*, **20**, 5946–5962, doi:10.1175/2007JCLI1664.1.

Starr, V. P., 1968: *Physics of Negative Viscosity Phenomena*. McGraw-Hill Book Company, 256 pp.

Stern, H. L., and R. E. Moritz, 2002: Sea ice kinematics and surface properties from RADARSAT synthetic aperture radar during the SHEBA drift. *J. Geophys. Res.*, **107**, 8028, doi:10.1029/2000JC000472.

Taylor, P. D., D. L. Feltham, P. R. Sammonds, and D. Hatton, 2006: Continuum sea ice rheology determined from subcontinuum mechanics. *J. Geophys. Res.*, **111**, C11015, doi:10.1029/2005JC002996.

Thorndike, A. S., 1986: Kinematics of sea ice. *The Geophysics of Sea Ice*, N. Untersteiner, Ed., Nato ASI Series, Vol. 146, Pleumum Press, 489–549.

Tremblay, L.-B., and L. A. Mysak, 1997: Modelling sea ice as a granular material, including the dilatancy effect. *J. Phys. Oceanogr.*, **27**, 2342–2360, doi:10.1175/1520-0485(1997)027<2342:MSIAAG>2.0.CO;2.

Ukita, J., and R. E. Moritz, 1995: Yield curves and flow rules of pack ice. *J. Geophys. Res.*, **100**, 4545–4557, doi:10.1029/94JC02202.

Weiss, J., E. M. Schulson, and H. L. Stern, 2007: Sea ice rheology form in-situ, satellite and laboratory observations: Fracture and friction. *Earth Planet. Sci. Lett.*, **255**, 1–8, doi:10.1016/j.epsl.2006.11.033.

Wilchinsky, A. V., and D. L. Feltham, 2004: A continuum anisotropic model of sea-ice dynamics. *Proc. Roy. Soc. London, A460*, 2105–2140, doi:10.1098/rspa.2004.1282.