On the calculation of normalized viscous–plastic sea ice stresses

Calculating and plotting the normalized states of stress for viscous–plastic sea ice models is a common diagnostic for evaluating the numerical convergence and the physical consistency of a numerical solution. Researchers, however, usually do not explain how they calculate the normalized stresses. Here, we argue that care must be taken when calculating and plotting the normalized states of stress. A physically consistent and numerically converged solution should exhibit normalized stresses that are inside (viscous) or on (plastic) the normalized yield curve. To do so, two possible mistakes need to be avoided. First, when using an implicit solver, normalized stresses should be computed from viscous coefficients and replacement pressure calculated using the previous numerical iterate and the strain rates at the numerator calculated from the latest iterate. Calculating the stresses only from the latest iterate falsely indicates that the solution has numerically converged. Second, for both implicit and explicit (i.e., the EVP) solvers, the stresses should be normalized by the ice strength and not by the replacement pressure. Using the latter, normalized states of stress only lie on the yield curve (i.e., falsely indicating there are no viscous states of stress).

1 Introduction

Sea ice deformations, associated with the formation of leads, pressure ridges and shear lines, strongly influence the evolution of the sea ice cover in both polar oceans. As they affect the thickness distribution, sea ice deformations have an important impact on the exchange of heat, moisture and momentum between the atmosphere and the underlying ocean. To properly represent these processes in a model, it is essential that rheology (i.e., the relation between applied stresses, material properties and resulting deformations) is correctly formulated.

Although some authors have recently proposed new sea ice rheologies (e.g., Girard et al., 2011), most sea ice models are still based on the viscous–plastic (VP) formulation introduced by Hibler (1979). With the VP rheology, the ice is treated as a very viscous fluid (creep flow) when the internal stresses are small. However, once the stresses reach critical values defined by a yield curve, the ice flows as a plastic material and large deformations (i.e., large spatial gradients of the velocity field) can occur.

Calculating and plotting the normalized states of stress with respect to the normalized yield curve is a useful diagnostic for assessing the physical consistency and numerical convergence of a VP solution. Indeed, this method can confirm whether a sea ice rheology is properly implemented in a model. The method is also helpful for evaluating numerical convergence. This is especially true for the explicit elastic-VP (EVP) solver (e.g., Hunke, 2001), which does not include a measure of convergence such as a residual. Unfortunately, researchers usually do not explain how they calculate this diagnostic (e.g., Zhang and Hibler, 1997; Hunke, 2001; Lemieux and Tremblay, 2009; Wang and Wang, 2009; Kimmritz et al., 2015). As demonstrated here, care must be taken when calculating the normalized stresses as two potential mistakes could lead to a misinterpretation of modeling results. The purpose of this article is to provide a short guide on how to calculate and plot the normalized states of stress for assessing physical consistency and convergence of numerical solutions.

2 The viscous–plastic sea ice rheology

With the Hibler (1979) VP rheology, the components σ_ij of the stress tensor are given by

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{ϵ}}}_{ij}+[\mathit{\zeta }-\mathit{\eta }]{\dot{\mathit{ϵ}}}_{kk}{\mathit{\delta }}_{ij}-P_p{\mathit{\delta }}_{ij}/\mathrm{2},i,j=\mathrm{1},\mathrm{2},\end{array}$$

where δ_ij is the Kronecker delta, ${\dot{\mathit{ϵ}}}_{ij}$ are the strain rates defined by ${\dot{\mathit{ϵ}}}_{\mathrm{11}}=\frac{\partial u}{\partial x}$, ${\dot{\mathit{ϵ}}}_{\mathrm{22}}=\frac{\partial v}{\partial y}$ and ${\dot{\mathit{ϵ}}}_{\mathrm{12}}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$ with u and v the components of the horizontal sea ice velocity vector, ${\dot{\mathit{ϵ}}}_{kk}={\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}$, ζ is the bulk viscosity, η is the shear viscosity and P_p is the ice strength (we follow the notation of Kreyscher et al., 2000).

The formulation of the viscosities depends on the yield curve and the flow rule. In the following, ζ and η are based on the widely used elliptical yield curve with a normal flow rule (Hibler, 1979):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{P_p}{\mathrm{2}\mathrm{\Delta }}},\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\eta }=\mathit{\zeta }{e}^{-\mathrm{2}},\end{array}$$

where $\mathrm{\Delta }={\left[({\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}{)}^{\mathrm{2}}+e^{-\mathrm{2}}({\dot{\mathit{ϵ}}}_{\mathrm{11}}-{\dot{\mathit{ϵ}}}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{e}^{-\mathrm{2}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}$, and e is the aspect ratio of the ellipse, i.e., the ratio of the long and short axes of the elliptical yield curve.

When Δ tends toward zero, Eqs. (2) and (3) become singular. To avoid this problem, Hibler (1979) proposed to limit the maximum values of viscosities which is equivalent to limiting the minimum value of Δ. Hence, ζ is expressed as

$$\begin{array}{}\text{(4)}& \mathit{\zeta }={\displaystyle \frac{P_p}{\mathrm{2}{\mathrm{\Delta }}^{*}}},\end{array}$$

where ${\mathrm{\Delta }}^{*}=max(\mathrm{\Delta },{\mathrm{\Delta }}_{\min })$ with ${\mathrm{\Delta }}_{\min }=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ s⁻¹. Note that other approaches for limiting the viscous coefficients have been proposed (e.g., Kreyscher et al., 2000; Lemieux and Tremblay, 2009).

A drawback of the standard VP rheology is that the term $-P_p{\mathit{\delta }}_{ij}/\mathrm{2}$ in Eq. (1) can cause the ice to deform even in the absence of forcing. To remedy this problem, $-P_p{\mathit{\delta }}_{ij}/\mathrm{2}$ is replaced by $-P{\mathit{\delta }}_{ij}/\mathrm{2}$, where P is a function of the strain rates. The simplest formulation of P is

$$\begin{array}{}\text{(5)}& P=P_p{\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}},\end{array}$$

where P tends toward zero for small deformations while it tends toward P_p for large deformations.

P is sometimes referred to as the replacement pressure (e.g., Hunke and Lipscomb, 2010). The use of a replacement method such as the one described above is now widely used in VP sea ice models (e.g., Wang and Wang, 2009; Losch et al., 2010; Hunke and Lipscomb, 2010).

3 The normalized yield curve

Using Eqs. (1), (3), (4), (5) and the definition of Δ, one can obtain

$$\begin{array}{}\text{(6)}& P_p^{\mathrm{2}}{\left({\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}}\right)}^{\mathrm{2}}={\left[{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}+P\right]}^{\mathrm{2}}+e^{\mathrm{2}}\left[({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}\right].\end{array}$$

Introducing the principal stresses σ₁ and σ₂ given by

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\pm \sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}},\end{array}$$

Eq. (6) becomes

$$\begin{array}{}\text{(8)}& P_p^{\mathrm{2}}{\left({\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}}\right)}^{\mathrm{2}}={\left[{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+P\right]}^{\mathrm{2}}+e^{\mathrm{2}}\left[({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}{)}^{\mathrm{2}}\right].\end{array}$$

As demonstrated below, the correct way to normalize the stresses in Eq. (8) is to divide them by the ice strength P_p, which leads to

$$\begin{array}{}\text{(9)}& {\left({\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}}\right)}^{\mathrm{2}}={\left[{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+P}{P_p}}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}}{P_p}}\right]}^{\mathrm{2}}.\end{array}$$

Defining ${\mathit{\sigma }}_{\mathrm{1}}^n={\mathit{\sigma }}_{\mathrm{1}}/P_p$ and ${\mathit{\sigma }}_{\mathrm{2}}^n={\mathit{\sigma }}_{\mathrm{2}}/P_p$, we obtain

$$\begin{array}{}\text{(10)}& {\left({\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}}\right)}^{\mathrm{2}}={\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+{\displaystyle \frac{P}{P_p}}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}},\end{array}$$

which describes a family of ellipses that depend on the ratio $\mathrm{\Delta }/{\mathrm{\Delta }}^{*}$ for their size and on the ratio P∕P_p for their center. Equation (10) with $\mathrm{\Delta }/{\mathrm{\Delta }}^{*}=P/P_p=\mathrm{1}$ defines what we refer to as the normalized yield curve in principal stress space. Hence, according to our rheology, normalized plastic stresses should fall on the normalized yield curve, while normalized viscous stresses should lie on smaller ellipses inside the normalized yield curve (Geiger et al., 1998).

4 Experimental setup

The divergence of the stress tensor (described in Sect. 2), which is ∇⋅σ, is one of the terms of the sea ice momentum equation. The momentum equation is discretized in space and in time (see for example Lemieux et al., 2012, for details). It is either solved implicitly with a Picard solver (e.g., Zhang and Hibler, 1997; Losch et al., 2010) or with a Newton solver (e.g., Lemieux et al., 2012; Losch et al., 2014; Mehlmann and Richter, 2017) or it is solved explicitly with the EVP approach (Hunke, 2001) or using the modified EVP with pseudo-time stepping (e.g., Kimmritz et al., 2015).

The numerical simulations for this paper were conducted with the Picard solver of the McGill sea ice model (see Lemieux and Tremblay, 2009, for details). The spatial resolution is 10 km and the time step is 30 min. All the experiments with the elliptical yield curve were done with the ice strength parameter P^* set to 27.5×10³ Nm⁻² and e=2. The model was restarted on 1 January 2002 at 12:00 UTC from a long-term simulation. The states of stress were calculated from solutions obtained at the first time level (i.e., 12:30 UTC). We will discuss later how our conclusions apply to the other types of solvers.

With a Picard solver, one has to solve a nonlinear system of equations that can be concisely written as A(u)u=b(u), where u is a vector that contains all the u and v velocity components on the grid, A is a sparse matrix and b is a vector that contains terms such as the atmospheric stress. It is important to mention that the elements of the matrix A depend on the viscous coefficients ζ and η and that the vector b contains the replacement pressure P. Implicit solvers such as Picard solve a series of linearized systems of equations in order to find the solution u of the nonlinear system of equations. This algorithm can be expressed as
1. Start with an initial iterate u⁰
do  $k=\mathrm{1},k_{max}$
2. Solve $\mathbf{A}\left({\mathit{u}}^{k-\mathrm{1}}\right){\mathit{u}}^k=\mathit{b}\left({\mathit{u}}^{k-\mathrm{1}}\right)$ with a linear
solver
3. Stop if $\left|\right|\mathit{F}\left({\mathit{u}}^k\right)\left|\right|<{\mathit{\gamma }}_{\mathrm{nl}}\left|\right|\mathit{F}\left({\mathit{u}}^{\mathrm{0}}\right)\left|\right|$
enddo,
where $\mathit{F}\left({\mathit{u}}^k\right)=\mathbf{A}\left({\mathit{u}}^k\right){\mathit{u}}^k-\mathit{b}\left({\mathit{u}}^k\right)$ is the residual at iteration k, the symbol $\left|\right|\left|\right|$ denotes the Euclidean norm and γ_nl<1 is the nonlinear convergence parameter. The iterations of this loop are referred to as nonlinear iterations or as in Lemieux and Tremblay (2009) as outer loop iterations. The maximum number of iterations k_max is set here to a very high value (10 000). A “fully” converged solution for u is characterized by a very small residual (γ_nl needs to be set to a value ≪1). As the stresses are function of u, a “fully” converged velocity vector leads to states of stress that are either on (plastic) or inside (viscous) the yield curve.

In order to shorten the article, the presentation of the algorithm above has been simplified. For numerical stability, the water stress should be linearized with $({\mathit{u}}^{k-\mathrm{1}}+{\mathit{u}}^{k-\mathrm{2}})/\mathrm{2}$ (Hibler and Ackley, 1983). Lemieux and Tremblay (2009) also linearized the rheology term with $({\mathit{u}}^{k-\mathrm{1}}+{\mathit{u}}^{k-\mathrm{2}})/\mathrm{2}$. For faster convergence of the Picard solver, we recommend to use $({\mathit{u}}^{k-\mathrm{1}}+{\mathit{u}}^{k-\mathrm{2}})/\mathrm{2}$ only for the water stress and to linearize the rheology term with u^k−1. Hence, it is important to notice that when linearizing the system of equation (in step 2), ζ, η and P are expressed as a function of u^k−1.

5 The calculation of normalized states of stress

The steps for calculating and plotting the normalized stresses are given below.

1. Solve the nonlinear system of equations for   u^k∼u

2. Calculate
   ${\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }\left({\mathit{u}}^{k-\mathrm{1}}\right){\dot{\mathit{ϵ}}}_{ij}\left({\mathit{u}}^k\right)+\left[\mathit{\zeta }\right({\mathit{u}}^{k-\mathrm{1}})-\mathit{\eta }({\mathit{u}}^{k-\mathrm{1}}\left)\right]{\dot{\mathit{ϵ}}}_{kk}\left({\mathit{u}}^k\right){\mathit{\delta }}_{ij}-P\left({\mathit{u}}^{k-\mathrm{1}}\right){\mathit{\delta }}_{ij}/\mathrm{2}$,  i,j=1,2

3. Calculate  
   ${\mathit{\sigma }}_{\mathrm{1}}^n,{\mathit{\sigma }}_{\mathrm{2}}^n=\frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}{P}_p}\pm \frac{\mathrm{1}}{P_p}\sqrt{{\left(\frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}$

4. Plot the ${\mathit{\sigma }}_{\mathrm{1}}^n,{\mathit{\sigma }}_{\mathrm{2}}^n$ using symbols such as circles

5. Plot the normalized yield curve 
   ${\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+\mathrm{1}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}}=\mathrm{1}$
    as a reference,

where the calculations in steps 2 and 3 should be done for all the ice-covered grid cells (here grid cells with a concentration larger than 0.5 are considered). The σ_ij (step 2) and ${\mathit{\sigma }}_{\mathrm{1}}^n,{\mathit{\sigma }}_{\mathrm{2}}^n$ (step 3) are calculated so that they are collocated at the tracer point of our model C-grid. Step 2 should be omitted for the standard and modified EVP; the time-stepped stresses should be used directly for step 3. Note that normalized stresses can also be plotted using the stress invariants ${\mathit{\sigma }}_{\mathrm{I}}=\left({\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}\right)/\mathrm{2}$ and ${\mathit{\sigma }}_{\mathrm{II}}=\left({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}\right)/\mathrm{2}$.

Following this method allows one to assess the physical consistency and the numerical convergence of the solution. We mean by physical consistency and numerical convergence that the states of stress are at their final position inside (viscous) or on (plastic) the yield curve. Many authors (e.g., Zhang and Hibler, 1997; Lemieux and Tremblay, 2009) have indeed shown that an approximate solution that has not sufficiently converged exhibits unrealistic states of stress that are outside the yield curve. This is shown in Fig. 1. For 2 (Fig. 1a) or 10 nonlinear iterations (Fig. 1b), the approximate solution has not converged and shows unrealistic states of stress. The fully converged solution (Fig. 1c) demonstrates physical consistency and numerical convergence. The fully converged solution was obtained by setting γ_nl to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$. Note that, in general, the fact that states of stress are on or inside the yield curve does not imply full convergence; the final positions (on and inside the yield curve) are obtained once u^k is the fully converged solution (Lemieux and Tremblay, 2009).

Two mistakes need to be avoided in order to obtain similar results as in Fig. 1 and therefore to be able to evaluate the numerical convergence of the solution and physical consistency.

First, one has to consider the way the nonlinear system of equations is solved. It is crucial to note that the σ_ij in step 2 should be calculated from ζ, η and P, which are a function of the previous iterate u^k−1, and the strain rates at the numerator from the latest iterate u^k. Let us consider that one calculates the stresses only based on the latest iterate u^k, i.e., the viscous coefficients ζ and η, and the replacement pressure are functions of u^k instead of u^k−1. Figure 2 shows the normalized states of stress that are obtained in this case after only two nonlinear iterations. One might conclude from this figure that the solution has converged as all the states of stress appear to be VP while we know this is not the case from Fig. 1a. This is important because a “true” converged solution exhibits better defined sea ice leads (and deformations; Lemieux and Tremblay, 2009), where large moisture–energy–salt fluxes are present between the sea ice, the ocean and the atmosphere.

This apparent numerical convergence of the solution is a consequence of the use of a rate-independent plastic rheology. This can be easily understood by considering a 1D VP example. Assuming that sea ice does not have tensile strength and that it exhibits a large convergent deformation, the 1D relation between the stress (σ) and the deformation ($\dot{\mathit{ϵ}}=\frac{\partial u}{\partial x}$) is given by

$$\begin{array}{}\text{(11)}& \mathit{\sigma }=\mathit{\zeta }\dot{\mathit{ϵ}}-{\displaystyle \frac{P}{\mathrm{2}}},\end{array}$$

where $\mathit{\zeta }=\frac{P_p}{\mathrm{2}\left|\dot{\mathit{ϵ}}\right|}$ and P=P_p for a large plastic deformation.

Correctly expressing ζ as a function of u^k−1 and $\dot{\mathit{ϵ}}$ as a function of u^k (with ${\dot{\mathit{ϵ}}}^k=\dot{\mathit{ϵ}}\left({\mathit{u}}^k\right)$), we obtain

$$\begin{array}{}\text{(12)}& \mathit{\sigma }={\displaystyle \frac{P_p}{\mathrm{2}\left|{\dot{\mathit{ϵ}}}^{k-\mathrm{1}}\right|}}{\dot{\mathit{ϵ}}}^k-{\displaystyle \frac{P_p}{\mathrm{2}}},\end{array}$$

which is equal to −P_p only once the numerical solution has converged.

On the other hand, expressing both ζ and $\dot{\mathit{ϵ}}$ as a function of u^k leads to

$$\begin{array}{}\text{(13)}& \mathit{\sigma }={\displaystyle \frac{P_p}{\mathrm{2}\left|\dot{{\mathit{ϵ}}^k}\right|}}{\dot{\mathit{ϵ}}}^k-{\displaystyle \frac{P_p}{\mathrm{2}}},\end{array}$$

which is always equal to −P_p whatever the velocity field u^k used.

A second possible mistake would be to normalize the principal stresses in step 3 with the replacement pressure P instead of using P_p. Indeed, dividing Eq. (8) by P², we get

$$\begin{array}{}\text{(14)}& \mathrm{1}={\left[{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+P}{P}}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}}{P}}\right]}^{\mathrm{2}}.\end{array}$$

Defining ${\mathit{\sigma }}_{\mathrm{1}}^n={\mathit{\sigma }}_{\mathrm{1}}/P$ and ${\mathit{\sigma }}_{\mathrm{2}}^n={\mathit{\sigma }}_{\mathrm{2}}/P$, we obtain

$$\begin{array}{}\text{(15)}& \mathrm{1}={\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+\mathrm{1}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}},\end{array}$$

which is the equation of an ellipse with a size and a center that are fixed. Equation (15) is in fact the same equation as the one for the normalized yield curve (i.e., Eq. 10 with $\mathrm{\Delta }/{\mathrm{\Delta }}^{*}=P/P_p=\mathrm{1}$). Simulated stresses normalized by P indeed converge toward this fixed ellipse. This is shown in Fig. 3 for 2 (a), 10 (b) and the fully converged solution (c). The converged normalized states of stress do not exhibit a realistic solution as all the stresses appear to be plastic.

Figure 1 Principal stresses normalized by the ice strength P_p after 2 (a), 10 (b) nonlinear iterations and the fully converged solution (c). ${\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+\mathrm{1}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}}=\mathrm{1}$ with e=2 is the normalized yield curve (solid black line).

Figure 2 Principal stresses after two nonlinear iterations calculated only from u^k and normalized by the ice strength P_p. The solution appears to be numerically converged because the σ_ij are only a function of u^k. ${\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+\mathrm{1}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}}=\mathrm{1}$ with e=2 is the normalized yield curve (solid black line).

Figure 3 Principal stresses normalized by the replacement pressure P after 2 (a), 10 (b) nonlinear iterations and the fully converged solution (c). ${\left[{\mathit{\sigma }}_{\mathrm{1}}^n+{\mathit{\sigma }}_{\mathrm{2}}^n+\mathrm{1}\right]}^{\mathrm{2}}+e^{\mathrm{2}}{\left[{\mathit{\sigma }}_{\mathrm{1}}^n-{\mathit{\sigma }}_{\mathrm{2}}^n\right]}^{\mathrm{2}}=\mathrm{1}$ with e=2 is the normalized yield curve (solid black line).

6 Broader considerations

The recommendations given above remain the same if another approach is used for limiting the viscous coefficients (see Eq. 4). Numerical experiments with the approach of Kreyscher et al. (2000) or with the hyperbolic tangent of Lemieux and Tremblay (2009) allow one to draw the same conclusions (not shown).

While it is not recommended to linearize the rheology term with the previous two iterates (as done by Lemieux and Tremblay, 2009), the stresses in step 2 (see beginning of Sect. 5) should in this case be obtained from ${\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }\left({\mathit{u}}_l\right){\dot{\mathit{ϵ}}}_{ij}\left({\mathit{u}}^k\right)+\left[\mathit{\zeta }\right({\mathit{u}}_l)-\mathit{\eta }({\mathit{u}}_l\left)\right]{\dot{\mathit{ϵ}}}_{kk}\left({\mathit{u}}^k\right){\mathit{\delta }}_{ij}-P\left({\mathit{u}}_l\right){\mathit{\delta }}_{ij}/\mathrm{2}$ with ${\mathit{u}}_l=({\mathit{u}}^{k-\mathrm{1}}+{\mathit{u}}^{k-\mathrm{2}})/\mathrm{2}$.

If one does not use a replacement pressure, the stresses in step 2 should be calculated the same way with P=P_p. Instead of lying on ellipses defined by Eq. (10), the normalized viscous states of stress would lie on concentric ellipses centered at ${\mathit{\sigma }}_{\mathrm{1}}^n={\mathit{\sigma }}_{\mathrm{2}}^n=-\mathrm{0.5}$ (Geiger et al., 1998).

As demonstrated below, our recommendations also apply when using other yield curves in a VP framework. As an example, additional numerical experiments were conducted with a Mohr–Coulomb yield curve with compressive capping (Ip et al., 1991). This different constitutive law is obtained by expressing the viscous coefficients and the replacement pressure as

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{P_p}{\mathrm{2}\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{*}\right|}},\text{(17)}& {\displaystyle }& {\displaystyle }P=\mathrm{2}\mathit{\zeta }\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|,\text{(18)}& {\displaystyle }& {\displaystyle }\mathit{\eta }={\displaystyle \frac{\left(\frac{P}{\mathrm{2}}-\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right)\sin \mathit{\varphi }}{\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathrm{s}}^{*}}},\end{array}$$

where ${\dot{\mathit{ϵ}}}_{\mathrm{I}}={\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}$ is the divergence, $\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{*}\right|=max\left(\right|{\dot{\mathit{ϵ}}}_{\mathrm{I}}|,d_{\min })$ with d_min a small deformation similar to Δ_min, ϕ is the angle of friction, ${\dot{\mathit{ϵ}}}_{\mathrm{s}}^{*}=max({\dot{\mathit{ϵ}}}_{\mathrm{s}},s_{\min })$ with ${\dot{\mathit{ϵ}}}_{\mathrm{s}}={\left[{\left(\frac{{\dot{\mathit{ϵ}}}_{\mathrm{11}}-{\dot{\mathit{ϵ}}}_{\mathrm{22}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}$ the maximum shear strain rate and s_min another small deformation, here set equal to d_min. In terms of the stress invariants, the Mohr–Coulomb failure criterion is simply written as ${\mathit{\sigma }}_{\mathrm{II}}=-{\mathit{\sigma }}_{\mathrm{I}}\sin \mathit{\varphi }$. This Mohr–Coulomb implementation assumes a pure shear flow rule. Divergence (larger than d_min) can only occur at the tip of the triangle and convergence when ${\mathit{\sigma }}_{\mathrm{I}}=-P_p$.

It is observed that with this new rheology, the Picard solver really struggles to obtain a numerically converged solution. With $P^{*}=\mathrm{27.5}\times {\mathrm{10}}^{\mathrm{3}}$ Nm⁻², $d_{\min }=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ s⁻¹ and sin ϕ=0.5 (i.e., $\mathit{\varphi }=\mathrm{30}{}^{\circ }$), the solver does not converge. When calculating the normalized stresses the correct way (as in step 2 in Sect. 5), there are states of stress outside the yield curve (not shown). However, similar to the results obtained with the elliptical yield curve (see Fig. 2), the normalized stresses (shown in Fig. 4 in stress invariant space) after two nonlinear iterations appear to have converged if only u^k is used for calculating the σ_ij.

Figure 4 Stress invariants after two nonlinear iterations calculated only from u^k and normalized by the ice strength P_p. The solution appears to be numerically converged because the σ_ij are only a function of u^k. The yield curve (solid black line) is based on a Mohr–Coulomb failure criterion with compressive capping.

To obtain a fully converged solution (with ${\mathit{\gamma }}_{\mathrm{nl}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$), P^*, d_min and sin ϕ were respectively set to 5×10² Nm⁻², $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻¹ and 0.01. Consistent with the results obtained with the ellipse, the converged stresses normalized by P_p are either on or inside the yield curve (not shown). Again, normalizing the converged stresses by the replacement pressure falsely indicates there are no stresses in the viscous regime (Fig. 5). Strangely, there are no states of stress on the long side of the triangle; all the states of stress appear to be at the tip and the short side of the triangle. This can be easily understood by using Eqs. (16) and (17) to calculate the normalized first stress invariant

$$\begin{array}{}\text{(19)}& {\mathit{\sigma }}_{\mathrm{I}}^n={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{I}}}{P}}={\displaystyle \frac{\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{I}}}{P}}-{\displaystyle \frac{P}{\mathrm{2}P}}={\displaystyle \frac{{\dot{\mathit{ϵ}}}_{\mathrm{I}}}{\mathrm{2}\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

which is equal to zero (tip of the triangle) when ${\dot{\mathit{ϵ}}}_{\mathrm{I}}>\mathrm{0}$ and equal to −1 (short side of the triangle) when ${\dot{\mathit{ϵ}}}_{\mathrm{I}}<\mathrm{0}$.

Figure 5 Fully converged stress invariants normalized by the replacement pressure P. The yield curve (solid black line) is based on a Mohr–Coulomb failure criterion with compressive capping.

7 Conclusions

We have described how the normalized states of stress should be calculated and plotted in order to assess the numerical convergence and physical consistency of a VP solution. To do so, modelers should avoid two possible mistakes.

First, to evaluate the numerical convergence of an approximate solution, one should calculate stresses from viscous coefficients and replacement pressure that are a function of the previous iterate u^k−1 and the strain rates at the numerator from the latest iterate u^k. This conclusion applies to all implicit solvers. As the EVP and modified EVP approaches include time-stepping equations for the stresses, one simply needs to calculate the normalized stresses from the stress outputs. This issue of misinterpretation of numerical convergence with normalized stresses is therefore more prone to occur with Picard and Newton solvers.

Second, the stresses should be normalized by the ice strength – not by the replacement pressure. Using the latter causes all the normalized stresses to lie on the normalized yield curve, falsely indicating there are no stresses in the viscous regime. This issue can affect the implicit solvers as well as the EVP and modified EVP approaches.

This article should serve as a guide on how to calculate and plot normalized VP states of stress for assessing physical consistency and convergence of numerical solutions. It also complements and gives more details about one of the sea ice diagnostics suggested for the CMIP6 sea-ice intercomparison project (Notz et al., 2016).