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
 ^{1}Recherche en Prévision Numérique Environnementale/Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
 ^{2}Service Météorologique Canadien, Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
 ^{1}Recherche en Prévision Numérique Environnementale/Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
 ^{2}Service Météorologique Canadien, Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
Correspondence: JeanFrançois Lemieux (jeanfrancois.lemieux@canada.ca)
Hide author detailsCorrespondence: JeanFrançois Lemieux (jeanfrancois.lemieux@canada.ca)
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).
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 elasticVP (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.
With the Hibler (1979) VP rheology, the components σ_{ij} of the stress tensor are given by
where δ_{ij} is the Kronecker delta, ${\dot{\mathit{\u03f5}}}_{ij}$ are the strain rates defined by ${\dot{\mathit{\u03f5}}}_{\mathrm{11}}=\frac{\partial u}{\partial x}$, ${\dot{\mathit{\u03f5}}}_{\mathrm{22}}=\frac{\partial v}{\partial y}$ and ${\dot{\mathit{\u03f5}}}_{\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{\u03f5}}}_{kk}={\dot{\mathit{\u03f5}}}_{\mathrm{11}}+{\dot{\mathit{\u03f5}}}_{\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):
where $\mathrm{\Delta}={\left[({\dot{\mathit{\u03f5}}}_{\mathrm{11}}+{\dot{\mathit{\u03f5}}}_{\mathrm{22}}{)}^{\mathrm{2}}+{e}^{\mathrm{2}}({\dot{\mathit{\u03f5}}}_{\mathrm{11}}{\dot{\mathit{\u03f5}}}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{e}^{\mathrm{2}}{\dot{\mathit{\u03f5}}}_{\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
where ${\mathrm{\Delta}}^{*}=max(\mathrm{\Delta},{\mathrm{\Delta}}_{\mathrm{min}})$ with ${\mathrm{\Delta}}_{\mathrm{min}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{9}}$ s^{−1}. 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
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).
Using Eqs. (1), (3), (4), (5) and the definition of Δ, one can obtain
Introducing the principal stresses σ_{1} and σ_{2} given by
Eq. (6) becomes
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
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
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).
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 pseudotime 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^{3} Nm^{−2} and e=2. The model was restarted on 1 January 2002 at 12:00 UTC from a longterm 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^{0}do
$k=\mathrm{1},\phantom{\rule{0.33em}{0ex}}{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\phantom{\rule{0.25em}{0ex}}\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}.
The steps for calculating and plotting the normalized stresses are given below.

Solve the nonlinear system of equations for
u^{k}∼u 
Calculate
${\mathit{\sigma}}_{ij}=\mathrm{2}\mathit{\eta}\left({\mathit{u}}^{k\mathrm{1}}\right){\dot{\mathit{\u03f5}}}_{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{\u03f5}}}_{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

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}}}$ 
Plot the
${\mathit{\sigma}}_{\mathrm{1}}^{n},{\mathit{\sigma}}_{\mathrm{2}}^{n}$using symbols such as circles

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 icecovered 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 Cgrid. Step 2 should be omitted for the standard and modified EVP; the timestepped 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 rateindependent 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{\u03f5}}=\frac{\partial u}{\partial x}$) is given by
where $\mathit{\zeta}=\frac{{P}_{p}}{\mathrm{2}\left\dot{\mathit{\u03f5}}\right}$ and P=P_{p} for a large plastic deformation.
Correctly expressing ζ as a function of u^{k−1} and $\dot{\mathit{\u03f5}}$ as a function of u^{k} (with ${\dot{\mathit{\u03f5}}}^{k}=\dot{\mathit{\u03f5}}\left({\mathit{u}}^{k}\right)$), we obtain
which is equal to −P_{p} only once the numerical solution has converged.
On the other hand, expressing both ζ and $\dot{\mathit{\u03f5}}$ as a function of u^{k} leads to
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^{2}, we get
Defining ${\mathit{\sigma}}_{\mathrm{1}}^{n}={\mathit{\sigma}}_{\mathrm{1}}/P$ and ${\mathit{\sigma}}_{\mathrm{2}}^{n}={\mathit{\sigma}}_{\mathrm{2}}/P$, we obtain
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.
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{\u03f5}}}_{ij}\left({\mathit{u}}^{k}\right)+\left[\mathit{\zeta}\right({\mathit{u}}_{l})\mathit{\eta}({\mathit{u}}_{l}\left)\right]{\dot{\mathit{\u03f5}}}_{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
where ${\dot{\mathit{\u03f5}}}_{\mathrm{I}}={\dot{\mathit{\u03f5}}}_{\mathrm{11}}+{\dot{\mathit{\u03f5}}}_{\mathrm{22}}$ is the divergence, $\left{\dot{\mathit{\u03f5}}}_{\mathrm{I}}^{*}\right=max\left(\right{\dot{\mathit{\u03f5}}}_{\mathrm{I}},{d}_{\mathrm{min}})$ with d_{min} a small deformation similar to Δ_{min}, ϕ is the angle of friction, ${\dot{\mathit{\u03f5}}}_{\mathrm{s}}^{*}=max({\dot{\mathit{\u03f5}}}_{\mathrm{s}},{s}_{\mathrm{min}})$ with ${\dot{\mathit{\u03f5}}}_{\mathrm{s}}={\left[{\left(\frac{{\dot{\mathit{\u03f5}}}_{\mathrm{11}}{\dot{\mathit{\u03f5}}}_{\mathrm{22}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\dot{\mathit{\u03f5}}}_{\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}}\mathrm{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^{−2}, ${d}_{\mathrm{min}}=\mathrm{2}\times {\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 u^{k} is used for calculating the σ_{ij}.
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^{2} Nm^{−2}, $\mathrm{1}\times {\mathrm{10}}^{\mathrm{8}}$ s^{−1} 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
which is equal to zero (tip of the triangle) when ${\dot{\mathit{\u03f5}}}_{\mathrm{I}}>\mathrm{0}$ and equal to −1 (short side of the triangle) when ${\dot{\mathit{\u03f5}}}_{\mathrm{I}}<\mathrm{0}$.
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 timestepping 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 seaice intercomparison project (Notz et al., 2016).
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 (Lemieux and Dupont, 2019). The Zenodo repository also includes the output files for the normalized stresses and the Matlab routine used for plotting.
JFL and FD derived the equations and designed the numerical experiments. JFL ran the experiments and wrote the article.
The authors declare that they have no conflict of interest.
We thank Philippe Blain, Amélie Bouchat, Bruno Tremblay and two anonymous reviewers for their helpful comments on the manuscript.
This paper was edited by Heiko Goelzer and reviewed by two anonymous referees.
Geiger, C. A., Hibler, W. D., and Ackley, S. F.: Largescale 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 seaice mechanics based on elastobrittle 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.: Viscousplastic 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. LACC06012, Los Alamos National Laboratory, 2010. a, b
Ip, C. F., Hibler, W. D., and Flato, G. M.: On the effect of rheology on seasonal seaice simulations, Ann. Glaciol., 15, 17–25, https://doi.org/10.3189/1991AoG1511725, 1991. a
Kimmritz, M., Danilov, S., and Losch, M.: On the convergence of the modified ElasticViscousPlastic 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 viscousplastic 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 Jacobianfree Newton Krylov method and the EVP model for solving the sea ice momentum equation with a viscousplastic 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 seaice 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 Jacobianfree NewtonKrylov solver for a coupled sea iceocean 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 viscousplastic 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 SeaIce Model Intercomparison Project (SIMIP): understanding sea ice through climatemodel simulations, Geosci. Model Dev., 9, 3427–3446, https://doi.org/10.5194/gmd934272016, 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