On the calculation of normalized viscous–plastic sea ice stresses
- 1Recherche en Prévision Numérique Environnementale/Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
- 2Service Météorologique Canadien, Environnement et Changement Climatique Canada, 2121 route Transcanadienne, Dorval, Quebec, Canada
Correspondence: Jean-François Lemieux (firstname.lastname@example.org)
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 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.
where δij is the Kronecker delta, are the strain rates defined by , and with u and v the components of the horizontal sea ice velocity vector, , ζ is the bulk viscosity, η is the shear viscosity and Pp is the ice strength (we follow the notation of Kreyscher et al., 2000).
where , 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
A drawback of the standard VP rheology is that the term in Eq. (1) can cause the ice to deform even in the absence of forcing. To remedy this problem, is replaced by , 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 Pp 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).
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 Pp, which leads to
Defining and , we obtain
which describes a family of ellipses that depend on the ratio for their size and on the ratio P∕Pp for their center. Equation (10) with 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 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×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
with a linear
3. Stop if
where is the residual at iteration k, the symbol 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 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 (Hibler and Ackley, 1983). Lemieux and Tremblay (2009) also linearized the rheology term with . For faster convergence of the Picard solver, we recommend to use 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.
The steps for calculating and plotting the normalized stresses are given below.
Solve the nonlinear system of equations foruk∼u
using symbols such as circles
Plot the normalized yield curve
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 (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 and .
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 . 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 (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 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; 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 () is given by
where and P=Pp for a large plastic deformation.
Correctly expressing ζ as a function of uk−1 and as a function of uk (with ), we obtain
which is equal to −Pp only once the numerical solution has converged.
On the other hand, expressing both ζ and as a function of uk leads to
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
Defining and , 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 ). 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 with .
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 (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 is the divergence, with dmin a small deformation similar to Δmin, ϕ is the angle of friction, with 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 . 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 .
It is observed that with this new rheology, the Picard solver really struggles to obtain a numerically converged solution. With Nm−2, s−1 and sin ϕ=0.5 (i.e., ), 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.
To obtain a fully converged solution (with ), P*, dmin and sin ϕ were respectively set to 5×102 Nm−2, 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
which is equal to zero (tip of the triangle) when and equal to −1 (short side of the triangle) when .
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 (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.: 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
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., 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
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