Geosci. Model Dev., 13, 1763–1769, 2020
https://doi.org/10.5194/gmd-13-1763-2020
Geosci. Model Dev., 13, 1763–1769, 2020
https://doi.org/10.5194/gmd-13-1763-2020

Methods for assessment of models 02 Apr 2020

Methods for assessment of models | 02 Apr 2020 # On the calculation of normalized viscous–plastic sea ice stresses

On the calculation of normalized viscous–plastic sea ice stresses
Jean-François Lemieux1 and Frédéric Dupont2 Jean-François Lemieux and Frédéric Dupont

Abstract

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).

Share
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., ), most sea ice models are still based on the viscous–plastic (VP) formulation introduced by . 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., ), 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., ). 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 VP rheology, the components σij of the stress tensor are given by

$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{ij}+\left[\mathit{\zeta }-\mathit{\eta }\right]{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{kk}{\mathit{\delta }}_{ij}-{P}_{p}{\mathit{\delta }}_{ij}/\mathrm{2},\phantom{\rule{1em}{0ex}}i,j=\mathrm{1},\mathrm{2},\end{array}$

where δij is the Kronecker delta, ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{ij}$ are the strain rates defined by ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}=\frac{\partial u}{\partial x}$, ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}=\frac{\partial v}{\partial y}$ and ${\stackrel{\mathrm{˙}}{\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, ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{kk}={\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}+{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}$, ζ is the bulk viscosity, η is the shear viscosity and Pp is the ice strength (we follow the notation of ).

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 (Hibler1979):

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

where $\mathrm{\Delta }={\left[\left({\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}+{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}{\right)}^{\mathrm{2}}+{e}^{-\mathrm{2}}\left({\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}-{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}{\right)}^{\mathrm{2}}+\mathrm{4}{e}^{-\mathrm{2}}{\stackrel{\mathrm{˙}}{\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, 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 }=\frac{{P}_{p}}{\mathrm{2}{\mathrm{\Delta }}^{*}},\end{array}$

where ${\mathrm{\Delta }}^{*}=max\left(\mathrm{\Delta },{\mathrm{\Delta }}_{\mathrm{min}}\right)$ with ${\mathrm{\Delta }}_{\mathrm{min}}=\mathrm{2}×{\mathrm{10}}^{-\mathrm{9}}$ s−1. Note that other approaches for limiting the viscous coefficients have been proposed (e.g., ).

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}\frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}},\end{array}$

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

P is sometimes referred to as the replacement pressure (e.g., ). The use of a replacement method such as the one described above is now widely used in VP sea ice models (e.g., ).

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(\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[\left({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}{\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}\right].\end{array}$

Introducing the principal stresses σ1 and σ2 given by

$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}}=\frac{{\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}±\sqrt{{\left(\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(\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[\left({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}{\right)}^{\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 Pp, which leads to

$\begin{array}{}\text{(9)}& {\left(\frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}\right)}^{\mathrm{2}}={\left[\frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+P}{{P}_{p}}\right]}^{\mathrm{2}}+{e}^{\mathrm{2}}{\left[\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(\frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{*}}\right)}^{\mathrm{2}}={\left[{\mathit{\sigma }}_{\mathrm{1}}^{n}+{\mathit{\sigma }}_{\mathrm{2}}^{n}+\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 PPp 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 .

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 , for details). It is either solved implicitly with a Picard solver (e.g., ) or with a Newton solver (e.g., ) or it is solved explicitly with the EVP approach (Hunke2001) or using the modified EVP with pseudo-time stepping (e.g., ).

The numerical simulations for this paper were conducted with the Picard solver of the McGill sea ice model (see , 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×103 Nm−2 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 u0
do  $k=\mathrm{1},\phantom{\rule{0.33em}{0ex}}{k}_{max}$
2. Solve $\mathbf{A}\left({\mathbit{u}}^{k-\mathrm{1}}\right){\mathbit{u}}^{k}=\mathbit{b}\left({\mathbit{u}}^{k-\mathrm{1}}\right)$ with a linear
solver
3. Stop if $||\mathbit{F}\left({\mathbit{u}}^{k}\right)||<{\mathit{\gamma }}_{\mathrm{nl}}||\mathbit{F}\left({\mathbit{u}}^{\mathrm{0}}\right)||$
enddo,
where $\mathbit{F}\left({\mathbit{u}}^{k}\right)=\mathbf{A}\left({\mathbit{u}}^{k}\right){\mathbit{u}}^{k}-\mathbit{b}\left({\mathbit{u}}^{k}\right)$ is the residual at iteration k, the symbol $||\phantom{\rule{0.25em}{0ex}}||$ 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 as outer loop iterations. The maximum number of iterations kmax 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 $\left({\mathbit{u}}^{k-\mathrm{1}}+{\mathbit{u}}^{k-\mathrm{2}}\right)/\mathrm{2}$ . also linearized the rheology term with $\left({\mathbit{u}}^{k-\mathrm{1}}+{\mathbit{u}}^{k-\mathrm{2}}\right)/\mathrm{2}$. For faster convergence of the Picard solver, we recommend to use $\left({\mathbit{u}}^{k-\mathrm{1}}+{\mathbit{u}}^{k-\mathrm{2}}\right)/\mathrm{2}$ only for the water stress and to linearize the rheology term with uk−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 uk−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   uku

2. Calculate
${\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }\left({\mathbit{u}}^{k-\mathrm{1}}\right){\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{ij}\left({\mathbit{u}}^{k}\right)+\left[\mathit{\zeta }\left({\mathbit{u}}^{k-\mathrm{1}}\right)-\mathit{\eta }\left({\mathbit{u}}^{k-\mathrm{1}}\right)\right]{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{kk}\left({\mathbit{u}}^{k}\right){\mathit{\delta }}_{ij}-P\left({\mathbit{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}}±\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., ) 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}×{\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 uk is the fully converged solution .

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 uk−1, and the strain rates at the numerator from the latest iterate uk. Let us consider that one calculates the stresses only based on the latest iterate uk, i.e., the viscous coefficients ζ and η, and the replacement pressure are functions of uk instead of uk−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; ), 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 ($\stackrel{\mathrm{˙}}{\mathit{ϵ}}=\frac{\partial u}{\partial x}$) is given by

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

where $\mathit{\zeta }=\frac{{P}_{p}}{\mathrm{2}|\stackrel{\mathrm{˙}}{\mathit{ϵ}}|}$ and P=Pp for a large plastic deformation.

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

$\begin{array}{}\text{(12)}& \mathit{\sigma }=\frac{{P}_{p}}{\mathrm{2}|{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}^{k-\mathrm{1}}|}{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}^{k}-\frac{{P}_{p}}{\mathrm{2}},\end{array}$

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

On the other hand, expressing both ζ and $\stackrel{\mathrm{˙}}{\mathit{ϵ}}$ as a function of uk leads to

$\begin{array}{}\text{(13)}& \mathit{\sigma }=\frac{{P}_{p}}{\mathrm{2}|\stackrel{\mathrm{˙}}{{\mathit{ϵ}}^{k}}|}{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}^{k}-\frac{{P}_{p}}{\mathrm{2}},\end{array}$

which is always equal to Pp whatever the velocity field uk used.

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

$\begin{array}{}\text{(14)}& \mathrm{1}={\left[\frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+P}{P}\right]}^{\mathrm{2}}+{e}^{\mathrm{2}}{\left[\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 1Principal stresses normalized by the ice strength Pp 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 2Principal stresses after two nonlinear iterations calculated only from uk and normalized by the ice strength Pp. The solution appears to be numerically converged because the σij are only a function of uk. ${\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 3Principal 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).

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 or with the hyperbolic tangent of 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 ), 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({\mathbit{u}}_{l}\right){\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{ij}\left({\mathbit{u}}^{k}\right)+\left[\mathit{\zeta }\left({\mathbit{u}}_{l}\right)-\mathit{\eta }\left({\mathbit{u}}_{l}\right)\right]{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{kk}\left({\mathbit{u}}^{k}\right){\mathit{\delta }}_{ij}-P\left({\mathbit{u}}_{l}\right){\mathit{\delta }}_{ij}/\mathrm{2}$ with ${\mathbit{u}}_{l}=\left({\mathbit{u}}^{k-\mathrm{1}}+{\mathbit{u}}^{k-\mathrm{2}}\right)/\mathrm{2}$.

If one does not use a replacement pressure, the stresses in step 2 should be calculated the same way with P=Pp. 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}$ .

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 . This different constitutive law is obtained by expressing the viscous coefficients and the replacement pressure as

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

where ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}={\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}+{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}$ is the divergence, $|{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}^{*}|=max\left(|{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}|,{d}_{\mathrm{min}}\right)$ with dmin a small deformation similar to Δmin, ϕ is the angle of friction, ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{s}}^{*}=max\left({\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{s}},{s}_{\mathrm{min}}\right)$ with ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{s}}={\left[{\left(\frac{{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{11}}-{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{22}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}$ the maximum shear strain rate and smin another small deformation, here set equal to dmin. In terms of the stress invariants, the Mohr–Coulomb failure criterion is simply written as ${\mathit{\sigma }}_{\mathrm{II}}=-{\mathit{\sigma }}_{\mathrm{I}}\mathrm{sin}\mathit{\varphi }$. This Mohr–Coulomb implementation assumes a pure shear flow rule. Divergence (larger than dmin) 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}×{\mathrm{10}}^{\mathrm{3}}$ Nm−2, ${d}_{\mathrm{min}}=\mathrm{2}×{\mathrm{10}}^{-\mathrm{9}}$ s−1 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 uk is used for calculating the σij. Figure 4Stress invariants after two nonlinear iterations calculated only from uk and normalized by the ice strength Pp. The solution appears to be numerically converged because the σij are only a function of uk. 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}×{\mathrm{10}}^{-\mathrm{8}}$), P*, dmin and sin ϕ were respectively set to 5×102 Nm−2, $\mathrm{1}×{\mathrm{10}}^{-\mathrm{8}}$ s−1 and 0.01. Consistent with the results obtained with the ellipse, the converged stresses normalized by Pp 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}=\frac{{\mathit{\sigma }}_{\mathrm{I}}}{P}=\frac{\mathit{\zeta }{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}}{P}-\frac{P}{\mathrm{2}P}=\frac{{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}}{\mathrm{2}|{\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}|}-\frac{\mathrm{1}}{\mathrm{2}},\end{array}$

which is equal to zero (tip of the triangle) when ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}>\mathrm{0}$ and equal to −1 (short side of the triangle) when ${\stackrel{\mathrm{˙}}{\mathit{ϵ}}}_{\mathrm{I}}<\mathrm{0}$. Figure 5Fully 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 uk−1 and the strain rates at the numerator from the latest iterate uk. 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 .

Code availability

Revision 333 of the McGill sea ice model, with some modifications for including the Mohr–Coulomb rheology, was used for the numerical experiments described in this article. The code is available on Zenodo at https://doi.org/10.5281/zenodo.3629542 . The Zenodo repository also includes the output files for the normalized stresses and the Matlab routine used for plotting.

Author contributions

JFL and FD derived the equations and designed the numerical experiments. JFL ran the experiments and wrote the article.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Philippe Blain, Amélie Bouchat, Bruno Tremblay and two anonymous reviewers for their helpful comments on the manuscript.

Review statement

This paper was edited by Heiko Goelzer and reviewed by two anonymous referees.

References

Geiger, C. A., Hibler, W. D., and Ackley, S. F.: Large-scale sea ice drift and deformation: Comparison between models and observations in the western Weddell Sea during 1992, J. Geophys. Res., 103, 21893–21913, https://doi.org/10.1029/98JC01258, 1998. a, b

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V.: A new modeling framework for sea-ice mechanics based on elasto-brittle rheology, Ann. Glaciol., 52, 123–132, https://doi.org/10.3189/172756411795931499, 2011. a

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. a, b, c, d

Hibler, W. D. and Ackley, S. F.: Numerical simulation of the Weddell Sea pack ice, J. Geophys. Res., 88, 2873–2887, https://doi.org/10.1029/JC088iC05p02873, 1983. a

Hunke, E. C.: Viscous-plastic sea ice dynamics with the EVP model: linearization issues, J. Comput. Phys., 170, 18–38, https://doi.org/10.1006/jcph.2001.6710, 2001. a, b, c

Hunke, E. C. and Lipscomb, W. H.: CICE: the Los Alamos sea ice model documentation and software user's manual version 4.1, Tech. Rep. LA-CC-06-012, Los Alamos National Laboratory, 2010. a, b

Ip, C. F., Hibler, W. D., and Flato, G. M.: On the effect of rheology on seasonal sea-ice simulations, Ann. Glaciol., 15, 17–25, https://doi.org/10.3189/1991AoG15-1-17-25, 1991. a

Kimmritz, M., Danilov, S., and Losch, M.: On the convergence of the modified Elastic-Viscous-Plastic method for solving the sea ice momentum equation, J. Comput. Phys., 296, 90–100, https://doi.org/10.1016/j.jcp.2015.04.051, 2015. a, b

Kreyscher, M., Harder, M., Lemke, P., and Flato, G. M.: Results of the Sea Ice Model Intercomparison Project: Evaluation of sea ice rheology schemes for use in climate simulations, J. Geophys. Res., 105, 11299–11320, https://doi.org/10.1029/1999JC000016, 2000. a, b, c

Lemieux, J.-F. and Dupont, F.: McGill sea ice model (Version rev. 333mc), Zenodo, https://doi.org/10.5281/zenodo.3629542, 2019. a

Lemieux, J.-F. and Tremblay, B.: Numerical convergence of viscous-plastic sea ice models, J. Geophys. Res., 114, C05009, https://doi.org/10.1029/2008JC005017, 2009. a, b, c, d, e, f, g, h, i, j

Lemieux, J.-F., Knoll, D. A., Tremblay, B., Holland, D. M., and Losch, M.: A comparison of the Jacobian-free Newton Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: a serial algorithm study, J. Comput. Phys., 231, 5926–5944, https://doi.org/10.1016/j.jcp.2012.05.024, 2012. a, b

Losch, M., Menemenlis, D., Campin, J.-M., Heimbach, P., and Hill, C.: On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations, Ocean Model., 33, 129–144, https://doi.org/10.1016/j.ocemod.2009.12.008, 2010. a, b

Losch, M., Fuchs, A., Lemieux, J.-F., and Vanselow, A.: A parallel Jacobian-free Newton-Krylov solver for a coupled sea ice-ocean model, J. Comput. Phys., 257, 901–911, https://doi.org/10.1016/j.jcp.2013.09.026, 2014. a

Mehlmann, C. and Richter, T.: A modified global Newton solver for viscous-plastic sea ice models, Ocean Model., 116, 96–107, https://doi.org/10.1016/j.ocemod.2017.06.001, 2017. a

Notz, D., Jahn, A., Holland, M., Hunke, E., Massonnet, F., Stroeve, J., Tremblay, B., and Vancoppenolle, M.: The CMIP6 Sea-Ice Model Intercomparison Project (SIMIP): understanding sea ice through climate-model simulations, Geosci. Model Dev., 9, 3427–3446, https://doi.org/10.5194/gmd-9-3427-2016, 2016. a

Wang, K. and Wang, C.: Modeling linear kinematic features in pack ice, J. Geophys. Res., 114, C12011, https://doi.org/10.1029/2008JC005217, 2009. a, b

Zhang, J. and Hibler, W. D.: On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res., 102, 8691–8702, https://doi.org/10.1029/96JC03744, 1997. a, b, c