New constitutive equations for the teardrop and parabolic lens yield curves in viscous-plastic sea ice models

,

Sea ice dynamical models are an integral part of the CMIP6 models  munity, 2020). Following the increase of computation power, the models' resolution was 45 formation (N. Hutter et al., 2019;Hutchings et al., 2005). LKFs are narrow bands where most of the sea ice deformation takes place (Kwok, 2001 with open water or thin ice where the bulk of the heat and matter transfer between at-49 mosphere and ocean takes place (Maykut, 1978), so that it important to represent LKFs  The VP rheological model requires the definition of a yield curve and a flow rule.

62
First, the yield curve sets the stress limit at which sea ice deforms plastically. Several

79
For useful climate modeling, sea ice models need to be stable and efficient, while giving 80 an accurate prediction of sea ice motion. The stability of the sea ice model can be ex-81 pressed by system energy considerations (Dukowicz, 1997;Schulkes, 1996;Pritchard, 2005).

82
For example, a negative bulk viscosity would make the rheology act as a spurious en-83 ergy source that leads to model instabilities. Instabilities lead to poor numerical con-84 vergence and inefficiency as the numerical convergence is more difficult to obtain.

85
In this paper, we document three issues in the formulation of the Teardrop (TD) 86 and Parabolic Lens (PL) yield curves with normal flow rules, as they are described in  prevents the stresses from remaining on or within the yield curve. We then propose so-91 lutions to these issues that ensure numerical convergence for high-resolution sea ice mod-92 els. We test these solutions in an idealized experiment and show that they lead to im-93 proved numerical convergence and efficiency.

94
The paper is structured as follows: Section 2 reviews the original TD and PL viscous-95 plastic rheologies (Zhang & Rothrock, 2005). Section 3 discusses the problems and so- Following general practice, we simulate sea ice as a (vertically integrated) 2D viscousplastic material. The ice velocities are calculated from the sea ice momentum equations: where ρ is the ice density, h is the mean sea ice thickness, u is the ice drift velocity field, 102 f is the Coriolis parameter, k is the vertical unit vector, τ a is the surface air stress, τ o 103 is the ocean drag, ∇Φ s is the acceleration due to the gradient of geopotential (i.e., sea 104 -4-manuscript submitted to JGR: Oceans surface) height, and σ is the vertically integrated internal ice stress tensor defined by 105 the sea ice VP constitutive equations. 106 We use the constitutive equation for the elliptical yield curve (Hibler, 1979(Hibler, , 1977: where ζ and η are the bulk and shear viscosities, and˙ ij are the strain rates defined as  The local ice strength P is defined as a function of the mean ice thickness h and concentration A, as where P is the compressive strength of 1 m thick ice and C is a model parameter defin-

110
ing the ice strength dependence on ice concentration.

111
Equation (2) can also be written in stress invariant form as where the stress invariants are and the strain rates invariants arė With the factor 1 2 for˙ I and˙ II , we follow the definition of the strain rate invariants given where q = 1 2 for the TD and q = 1 for the PL yield curve. P is the local isotropic com- where Equations (9) or (10), together with the associated flow rule conditions, constitute systems of three equations and three unknowns (σ I , σ II , λ) or (u, y, λ). The solution of these systems can be written as (15) and where l =˙ İ II . Note that a second root for x TD is discarded because it leads to flow rules pointing inward of the yield curve (not shown). The bulk and shear viscosities ζ and η can be calculated from the above constitutive equations (Eqs. 4 -5) as, and for the PL. (21) In the limit where deformations tend to zero, the viscous coefficients η and ζ be- denominator have different signs, i.e., for: (the pressure term), we obtain, By comparing this formulation to Eq (4), new definitions for ζ TD and ζ PL emerge: and the pressure term in its original form P 2 becomes 2−kt 3 P and 1−kt 2 P . These where α ( 1, e.g., α = 0.95) ensures that η is never equal to zero; i.e., σ I = T or −P .   190 To overcome this issue, we impose limits of ζ and η as 191 ζ = min ζ, ζ max min 1, ζ η , where ζ max = η max are model parameters corresponding to ζ max for the elliptical yield