Teardrop and Parabolic Lens Yield Curves for Viscous-Plastic Sea Ice Models: New Constitutive Equations and Failure Angles

,

The VP rheological model requires the definition of a yield curve and a flow rule; in the following, we refer to a specific yield curve and flow rule as a rheology. The yield curve sets the stress limit at which sea ice deforms plastically. Several yield curve shapes have been proposed in the context of a viscous-plastic framework: elliptical (Hibler, 1979), the sine-wave lens (Bratchie, 1984), triangular or Mohr-Coulomb (Hibler & Schulson, 2000;Ip et al., 1991;Tremblay & Mysak, 1997), and TD or PL (Zhang & Rothrock, 2005). The flow rule then sets the relative amount of shear and divergence (positive or negative) along LKFs for a given stress state. The flow rule can be normal (Hibler, 1979;Zhang & Rothrock, 2005) or non-normal (Ip et al., 1991;Ringeisen et al., 2021) to the yield curve.
As a contribution to the Arctic Ice Dynamic Joint Experiment AIDJEX, Rothrock (1975) proposed two different yield curves: the TD and PL. Both yield curves have a more realistic dependence of shear strength dependence on normal stresses-in line with granular material behavior-and satisfy Drucker's convexity postulate for stability (Drucker, 1950;Palmer et al., 1967), contrary to the Mohr-Coulomb yield curve, typically used to simulate granular media. Akin to the elliptical yield curve, the TD and PL yield curves also represent both divergence and convergence, in contrast with the Mohr-Coulomb yield curve with a normal flow rule when the stress states are located on its limbs. These two yield curves have been implemented by Zhang and Rothrock (2005), and while being used in the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) model (Zhang, 2020), are not used elsewhere in the community.
A common problem to sea ice models using any rheological framework (e.g., VP, EAP, EB-family) is that they simulate wide intersection angles between conjugate pairs of LKFs when compared with observations -although some improvement can be made by changing the aspect ratio of the ellipse in the VP model (Bouchat & Tremblay, 2017;Hutter et al., 2022). The intersection angle between LKFs is important and has a direct impact on anisotropy in the ice field (an emergent property) and therefore future deformation. Idealized numerical experiments with the commonly used VP model with elliptical yield curves show that the intersection angles can be linked theoretically to the shape of the yield curve for normal or non-normal flow rules (Ringeisen et al., 2019(Ringeisen et al., , 2021. Measurements from laboratory experiments show yield stresses and small failure angles for saline ice samples that are in agreement with a Mohr-Coulomb or TD yield curve (Iliescu & Schulson, 2004;Richter-Menge & Jones, 1993;Schulson, 2004) and observations of a single angle of fracture (Erlingsson, 1988). Yet, the simulated failure angles for these non-elliptical yield curves have not been investigated.
Sea ice models are computationally expensive because of the highly nonlinear rheological equations (Koldunov et al., 2019;Lemieux & Tremblay, 2009;Losch & Danilov, 2012). For long-term climate simulations, sea ice models must be stable and efficient while producing a precise (converged) solution of the sea ice momentum equation The stability of the sea ice model can be discussed in terms of energy considerations (Dukowicz, 1997;Pritchard, 2005;Schulkes, 1996). For instance, negative viscosities in the constitutive equations are a spurious energy source that leads to numerical instabilities. Instabilities, in turn, prevent full convergence of the numerical solver, and therefore a proper representation of LKFs.
In this paper, we propose new formulations of the sea ice constitutive equations for the TD and Parabolic Lens (PL) yield curves with normal flow rules first derived by Zhang and Rothrock (2005). We then discuss how the new formulation solves three issues related to numerical convergence in the original formulation of the TD and PL yield curves. First and second, stress states on some parts of the yield curve had viscosity coefficients that are negative or zero, resulting in spurious energy sources of energy as opposed to sinks-The rheology is a dissipative term in the momentum equation. Third, the capping of the non-linear viscosities in the viscous regime Journal of Advances in Modeling Earth Systems RINGEISEN ET AL. 10.1029/2023MS003613 3 of 13 leads to reduced numerical convergence and stress states that lie outside of the yield curve. Interestingly, those stress states are both viscous and plastic in convergence and shear, or vice versa. The proposed new formulations significantly improve numerical convergence for sea ice models. We test the new formulations in an idealized uni-axial loading experiment and show that they lead to smaller angles of failure between conjugate pairs of LKFs and suggest that they are possible candidates to replace the elliptical yield curve.
The paper is structured as follows: Section 2 presents the sea ice model and the new derivation of the TD and PL viscous-plastic rheologies. Section 3 discusses the differences with the formulation from Zhang and Rothrock (2005) and how it created issues impeding numerical convergence. Section 4 presents the experiments used to compare the new and original rheology formulations. Section 5 compares the numerical convergence between the original and new formulations and investigates the failure angles created with the TD and PL yield curves. The discussion is presented in Section 6, followed by the Summary and conclusions in Section 7.

The Sea Ice Viscous-Plastic Rheological Model
The sea-ice momentum equation for a (vertically integrated) 2D viscous-plastic material can be written as: where ρ is the ice density, h is the mean sea ice thickness, u is the ice drift velocity field, f is the Coriolis parameter, k is the vertical unit, τ a is the surface air stress, τ o is the ocean drag, ∇Φ s is the acceleration due to the gradient of geopotential height (i.e., sea surface variations), and σ is the vertically integrated internal ice stress tensor defined by the sea ice VP constitutive equations.
The general form of the constitutive equations for non-elastic viscous material (or Reiner-Rivlin equations, Astarita & Marrucci, 1974) can be given by: where ζ, η, and p are the bulk viscosity, the shear viscosity, and an isotropic pressure term. The specific form of these terms is defined by details of the yield curve and the flow rule. The terms are the strain rates defined as . The nonlinear relationships between the viscosities and the strain rates take different forms to define the sea ice VP rheologies. Equation 2 can also be written in terms of invariants as where the stress invariants are I = and the strain rate invariants arėI The factor 1 2 for I and II follows the definition of the strain rate invariants given in Zhang and Rothrock (2005). The local isotropic ice compressive strength P, commonly used in the definition of ζ, η, and p, is defined as a function of the mean ice thickness h and concentration A, as where P ⋆ is the compressive strength of 1 m thick ice and C ⋆ is a model parameter defining the ice strength dependence on ice concentration.

Journal of Advances in Modeling Earth Systems
RINGEISEN ET AL.

Derivation of the TD and PL Viscosities
The TD and PL yield curves can be written as (Zhang & Rothrock, 2005): for the TD and PL yield curve, respectively, P is the local isotropic ice compressive strength given by Equation 7, k t = T/P is the tensile factor, T is the local isotropic ice tensile strength (König Beatty & Holland, 2010), and a in the original formulation is replaced by k t for consistency with more recent literature.
Equation 8, together with the associated normal flow rule conditions, constitute the TD and PL systems of three equations and three unknowns (σ I , σ II , λ). The solution of these systems can be written as and where =̇İ II . Note that a second root for I,TD is discarded because it leads to flow rules pointing inward of the yield curve (not shown).
From this point, our derivation differs from the one of Zhang and Rothrock (2005), see Section 3.1. Rewriting Equations 10 and 13 as a linear function of the divergence I and a strain-rate independent part (the pressure term p in Equation 3), we obtain, By comparing this formulation with that of Equation 3, the definition of ζ TD and ζ PL emerges from the system of equations, as For the shear viscosity η, we proceed as in Zhang and Rothrock (2005): we insert Equations 11 and 14 in Equation 4 to get for the TD, and for the PL. In the above, we write η in terms of σ I -as opposed to stain rates only-to simplify the notation.
The following regularizing conditions on σ I ensure that the TD and PL yield curves are C 0 -continuous and avoid a gap in the yield curve (see Section 3.2 and Figure where α should be near but less than 1 (e.g., α = 0.95) to avoid zero shear viscosity η, see Equations 20 and 21. A minor consequence of these new conditions is that no stress states are present near the tips of the yield curve. A major benefit is a significant improvement of the numerical convergence (Section 3.2). Similar conditions to avoid zero shear viscosity η at the tips were used with the TD and PL yield curve in Zhang and Rothrock (2005) but were not described in the paper (J. Zhang, personal communication).
When the deformation tends to zero, the viscosities tend to infinity. Following Hibler (1977), we cap the viscosities to some maximum values to transition from the plastic regime (stresses lie on the yield curve, and the viscosities are inversely proportional to strain rate) to the viscous regime (the stresses increase linearly with the strain rate and constant viscosities). In contrast with the elliptical yield curve, the shear viscosity is not a fixed fraction of the bulk viscosity and independent variable upper bounds for ζ and η are required: where ζ max and η max are free model parameter. In the following, we use ζ max = η max for simplicity. In the absence of these restrictions, the stress states can be plastic in shear and viscous in convergence (or vice versa). While this appears intuitively correct (e.g., viscous in convergence and plastic in shear along an LKF), this leads to a significant degradation of the numerical convergence. With the above relation, we ensure that all stress states inside the yield curve represent viscous deformations and that all stress states on the yield curve represent plastic deformations, see Section 3.3.

Journal of Advances in Modeling Earth Systems
RINGEISEN ET AL.
10.1029/2023MS003613 6 of 13 Zhang and Rothrock (2005) assumed the constitutive relations of Hibler (1979), Equation 26, is applicable to the (non-symmetrical, tensile-capable) TD and PL yield curves. Consequently, the first stress invariant and bulk viscosities for the TD yield curve are written as (using Equations 10 and 13):

Negative and Zero Bulk Viscosities
The shear viscosities η TD and η PL are computed the same way as Equations 20 and 21.
In the above, with the definition of I from Equations 10 and 13, the bulk viscosity ζ is negative when the numerator and denominator have different signs in the ranges that are marked in red regions on the TD and PL yield curves in Figure 1.  Inside region (1), the three colors indicate three behaviors: Stress states marked with a red bracket exceed the yield curve because of the viscous capping in ζ, and these stresses are displaced horizontally toward the 2 line. The green brackets indicate cases where the bulk viscosity is capped at ζ min (= 0) but the shear viscosity is not (no resistance to convergence and plastic in shear). The blue bracket (along the 2 line and inside of the yield curve) are cases where the bulk viscosity is capped to ζ min (=0) and the shear viscosity is capped at η max (no resistance to convergence and viscous in shear). Box (2) shows the gap in the yield curve created by the conditions described in Equations 31 and 32. (b) Trajectories of the stresses from before viscous capping to viscous capping in the original (orange, lower half plane) and new formulation (blue, upper half plane) (Section 3.3). In the original formulation (independent capping of η and ζ), the stresses migrate horizontally or vertically as the viscous bounds are applied independently and then in both directions when both coefficients are capped. In the new formulation, (synchronous capping of η and ζ, Equations 24 and 25), the stresses migrate toward the center of the TD yield curve A negative ζ is problematic as it makes the rheology term an energy source and the sea ice model instantly blows up. In Zhang and Rothrock (2005), the negative ζ are capped to ζ min = 0 (J. Zhang, personal communication). This leads to a gap in the yield curve before viscous bounds are prescribed on ζ (Gap in the green line, Region (1), red brackets, on Figure 2a), and stress states that lie along the vertical σ I = P/2 line both outside and inside of the yield curve near its top (Region (1), green and blue brackets, respectively, on Figure 2a). Recall that stress states must lie within (viscous regime) or on (plastic regime) the yield curve.
The root cause of the negative bulk viscosities is the assumption that the constitutive relation of the elliptical yield curve (with axial symmetry around I = − 2 ) and no tensile stress is also applicable to the non-symmetrical TD or symmetrical PL yield curves with normal flow rules but non-zero tensile strength. Note that the range of σ I between ( 2 ) (original formulation) and

Zero Shear Viscosity
In Zhang and Rothrock (2005), the following conditions are used for the tip(s) of the TD and PL yield curves: The reason why stress states for which the normal is such that with |l| ≥ 1 were not allowed is unclear. Note that such regions only exist when k t > 0. In addition, conditions 31 and 32 are not required for the mathematical derivation of Equations 10 and 13. In fact, they lead to another gap in principle stress space where no stress condition is defined (region (2) on Figure 2a).
Also, when the conditions for capping on σ I are met (Equations 31 and 32), the shear viscosity η and shear stress σ II are identically zero (Equations 20 and 21). While a zero shear viscosity is not an energy source, it still has a significant impact on the numerical stability given that most energy dissipation by the rheology is taking place in shear (Bouchat & Tremblay, 2014). Note that even small shear viscosities are sufficient to stabilize the model during large deformation events.

Capping of Bulk and Shear Viscosities
In the limit where deformations tend to zero, the viscous coefficients η and ζ become infinite. To avoid this situation, η and ζ must be bounded from above to some maximum allowable value (ζ max and η max ). In the original formulation of Zhang and Rothrock (2005), and akin to the viscous plastic sea ice model of Hibler (1979), ζ and η are capped independently to η max and ζ max (J. Zhang, personal communication). Consequently, the simulated state of stress can be viscous in convergence or divergence and plastic in shear, or vice versa. While this is not unphysical-for instance, sea ice deformation can be viscous in convergence and plastic in shear along an LKF-mixed mode of deformation in shear and convergence has a negative impact on numerical convergence. Conversely, when the shear and bulk viscosities are both capped at the same time, as done in the new formulation, an improvement in the residual of around half an order of magnitude can be achieved (results not shown).

Comparison of the Original and New Formulation
Figure 2 shows plastic and viscous stress states with the new formulation described in Section 2.2 and the original formulation of Zhang and Rothrock (2005) of the TD yield curve. To create these stress states, we create a random field of deformation rates , which contains all combinations of shear, compression, and tension. We then apply the constitutive equation to compute the σ ij stresses. Note that the stress states of both formulations are therefore calculated with the same strain rates. The magnitude of the random strain rates is set to ensure both viscous and plastic states.
With the original simulation, the stress states for which ζ = 0 are outside of the yield curve and gather along the I = 2 line-region (1) are independently capped-region (1), red and blue bracket respectivelyand there are no stress states on the yield curve close to the tip-region (2). These two features do not appear in our modified formulation of the yield curve. The comparison with the PL yield curve shows the same features (not shown).

Experimental Setup
Following Ringeisen et al. (2019), we test the new formulation for the TD and PL yield curve (and compare with the original formulation) in a simple uniaxial loading experiment using the MIT general circulation model (MITgcm, Campin et al., 2021). To this end, an ice floe of 60 km by 250 km is embedded in a 100 km by 260 km domain with a constant grid spacing of 1 km (see Figure 3). The initial condition for sea ice thickness, concentration and velocity are h = 1 m, A = 100% and u ice = 0. The ice floe is in direct contact with the southern boundary (no slip boundary condition) and centered in the domain laterally with two open water bands of 20 km on each side, and one of 10 km at the top of the domain to eliminate the potential effect of boundary conditions on the angle of fracture.
We use two types of forcing (surface and boundary) for two different numerical experiments (Circled numbers on Figure 3): 1. For the numerical convergence study (Section 5.1), we use a uniform southward surface stress of 0.15 N m −2 and a Picard solver for the solution of the full non-linear momentum equation with 15,000 outer-loop iterations (or pseudo-timesteps) (Zhang & Hibler, 1997). For each of the solver outer-loop iterations, we use a LSR solver for the linearized momentum equation until the solution reaches an accuracy LSR = max(| − −1|) = 10 −9 m s −1 or 15,000 iterations, whichever comes first. The timestep is 10 s and the total integration time is 200 s or 20 timesteps. In all simulations, we use a tensile factor k t = 0.05, unless specified otherwise. 2. For the study of the failure angles between conjugate pairs of linear kinematic features, we prescribe a linearly varying ice velocity at the northern border (from v i = 0 to v i = 0.1 m s −1 ) in order to have a spatially uniform stress state in the ice field (see also Ringeisen et al. (2019)). We also use a Picard solver with 1,500 outer-loop iterations (or pseudo-timesteps) (Zhang & Hibler, 1997), unless otherwise specified. And we use the LSR solver with a ɛ LSR = 10 −11 m s −1 accuracy or 1,500 iterations, whichever comes first. The timestep is 0.1 s and the experiment total length is 5 s.

Numerical Convergence
Hereafter, we refer to the L2-norm of the relative residuals of the non-linear Picard solver as the "non-linear residual." The non-linear residual for the new formulation decreases by nearly three orders of magnitude compared with that of the original formulation where the L2-norm increases despite the large maximum allowable number of outer-loop iterations (Figure 4). The new formulation, nevertheless, still has poorer convergence properties when compared with the more diffusive elliptical yield curve, that is, large internal ice stresses even for low mean normal stresses σ I . The stress states for the new TD formulation are mostly on or within the yield curve, a necessary but not sufficient condition indicating the convergence of the numerical solver (Lemieux & Tremblay, 2009), in contrast with the original formulation that has a large number of points far outside particularly for low mean normal stresses. Moreover, the zero shear stress states (σ II /P = 0) in the new formulation disappear ( Figure 5). We see similar improvements in numerical convergence and position of stress states within the yield curve with the PL yield curve (not shown).
These results are robust to the exact choice of the model domain, surface forcing, and solver parameters. For instance, when using less strict solver parameters (10 outer-loop iterations, 500 LSR iterations, and ɛ LSR = 10 −6 , typical for high-resolution (1-2 km) pan-Arctic sea ice simulations (e.g., Hutter et al., 2018), the reduction in non-linear residuals with the number of outer-loop iteration in both the original and new formulation is similar, and still larger than that of the elliptical yield curve (results not shown). It appears that the price to pay for a more realistic dependency of the shear stress to the normal stress is a higher number of outer-loop iterations.
Finally, we observe that the failure angle no longer depends on the exact rheological model when the numerical convergence criterion is relaxed (results not shown). In addition, with only two outer-loops (Zhang, 2020;Zhang & Hibler, 1997), 500 LSR iterations or ɛ = 10 −6 , the numerical convergence of both the original and new formulation is nearly the same, and the absolute non-linear residual remains high (Res ≃ 10 2 ); the number of linear iterations of the new formulation, however, is reduced (results not shown).

Theoretical Failure Angle
The formulation of the TD and PL yield curve use a normal flow rule, so both the Roscoe (1970) and Coulomb (1776) theories predict the same LKFs orientation (Vermeer, 1990). Because the loading is uni-axial and in compression, the first principle stress is equal to zero, and the two stress invariants are of opposite sign; therefore ′ II = − ′ I , where ′ I = I and ′ II = II (Ringeisen et al., 2019). For the TD yield curve, it means that This equation can be solved using Vieta's substitution for ′ I,s , the only real solution (Bronshtein et al., 2015). The reader is referred to Appendix A for a full derivation.
The theoretical failure angle is given by Ringeisen et al. (2019): where the slope of the TD yield curve is given by   Similarly, for the PL yield curve, the first stress invariant at failure, the slope of the yield curve at failure, and the angle of failure can be written as:

Sensitivity of the Failure Angles to Tensile Strength
We investigate the dependence of the failure angle θ on the tensile factor k t , the only free parameter in the TD and PL rheologies besides the isotropic compressive strength P. Diamond shape failure lines, or LKFs, are simulated using both rheologies, akin to other yield curves and flow rules (Ringeisen et al., 2019(Ringeisen et al., , 2021. For k t < 0.05, the uni-axial failure angles are below 30°, the minimum value for the elliptical yield curve with a normal flow rule for any ellipse aspect ratio. The angle of failure increases with k t , in agreement with the theoretical predictions ( Figure 6). The simulated angles of failure agree with the theoretical predictions derived in Section 5.2.1: For the TD rheology, the RMS error and the R 2 number between the simulated angles and theory are 0.46° and 0.992. For the PL rheology, the RMS error is 0.24°, and the R 2 is 0.998. Such an agreement is expected as the TD and PL yield curves both have a normal flow rule, as for the elliptical yield curve (Ringeisen et al., 2019). Note that for the same value of k t , the failure angles of the PL are always larger than those with the TD.

Discussion
We identified three issues in the derivation of the original TD and PL yields curve formulations that affect the convergence properties of the numerical solver if full convergence is attempted with a large number of outer loop iterations. To the authors' knowledge, however, the TD yield curve is only used in the PIOMAS model with (2 pseudo-timesteps, Zhang & Hibler, 1997), or outer loop iteration. Investigating the formulations of sea ice models in idealized simulations where numerical convergence can be attempted makes it easier to identify, investigate, and solve numerical issues with a given formulation.
The reasons for the slower convergence with the TD compared to the elliptical yield curve remains unclear. Possible causes include smaller (when compared with the elliptical yield curve) non-linear viscosities for small normal stress, a larger area of the yield curve that has low curvature leading to a smaller dependency of the internal stresses on the strain-rates (via the normal flow rule), and non-constant viscosities in the viscous regime inside the yield curve. For the elliptical yield curve with P = const, the stress states within the yield curve always have the same viscosities (ζ = ζ max and = max 2 ). With the TD and PL yield curves, the viscosities in the viscous regime are defined independently of one another and the transition from plastic to viscous deformations, as a consequence, also has to be treated independently. For small stresses (i.e., σ I and σ II ≃ 0), this leads to unphysical gradients of shear viscosity. These gradients lead to spatial stress variations, thus small-scale viscous deformations. These variations in deformations make the numerical convergence more difficult for the TD compared to the elliptical yield curve.
Both rheologies show good agreement with the theoretical prediction of failure angles with a normal flow rule. The TD and PL rheologies have two advantages: (a) They create angles smaller than the elliptical yield curve (Ringeisen et al., 2019) and in the range of observations (Erlingsson, 1988;Hutter et al., 2022) when used with a small tensile strength factor (k t < 5%) and (b) they feature a normal flow rule, which is simpler to solve numerically compared to the elliptical yield curve with a non-normal flow rule that was introduced to reduce angles between LKFs (Ringeisen et al., 2021). For the TD and PL yield curves, we see two drawbacks: (a) both yield curves only have two free parameters, the compressive strength P and the tensile factor k t , preventing adjustment of the shear strength of the material for a given normal stress, contrary to the elliptical yield curve where shear, compressive and tensile strength can be adjusted independently and (b) More surface of the TD yield curve is associated with sea ice divergence, while the elliptical and the PL yield curve have equal surface with sea ice divergence and convergence, in agreement with observations (Stern et al., 1995).

Summary and Conclusions
New formulations of the TD and PL yield curves with normal flow rule addresses three issues in their original formulations (Zhang & Rothrock, 2005). Two of these issues have a common cause: zero bulk or shear viscosity. The zero bulk viscosity ζ is a consequence of the assumption that the constitutive equation Equation 2 presented in Hibler (1979) is valid for all yield curves, including asymmetrical yield curves with respect to an average internal pressure 2 and yield curves with non-zero tensile strength. A similar but consistent constitutive equation for both TD and PL emerges from the yield curve equation and the normal flow rule conditions: the form of the viscous coefficient changes and the pressure term has a new scaling factor function of the tensile strength parameter k t . In this new formulation, all bulk and shear viscosities are positive removing the need for ad-hoc capping of the viscous coefficients.
The third identified issue is less important but may lead to unphysical states and poorer convergence (not shown). By consistently capping the bulk and shear viscosities for the transition from plastic to viscous states, we avoid stress states with mixed type (plastic and viscous) that while being physical lead to poorer numerical convergence properties. With this change, the stress states are guaranteed to lie on or within the yield curve and the numerical convergence is also improved.
We show that for small values of the tensile strength factors k t < 0.05 the TD and PL yield curves can create smaller angles than the elliptical yield curve in uniaxial compression. The shape of the TD yield curve resembles the shape of the Mohr-Coulomb yield curve and eliminates non-differentiable points when compared with other yield curves like the curved diamond yield curve (Wang, 2007) or the Coulombic yield curves (Hibler & Schulson, 2000). This makes the TD yield curve an interesting alternative for general use in the community.

Appendix A: Solving for the Angle for the TD Yield Curve
We solve the cubic equation that determines the coordinates of the intersection point between the teardrop yield curve and the σ II = −σ I axis. This equation of the cubic form x 3 + a 2 x 2 + a 1 x + a 0 is written ′3 I − 2 ′2 I + ( 2 − 2 ) ′ I − 2 = 0.
We solve this equation using Vieta's substitution (Bronshtein et al., 2015): we define Q and R If D < 0, then three distinct real solutions exist. Note that, because we compute the intersection of two curves (the axis σ II = −σ I , and the yield curve), the solutions are real and not complex. The three real solutions are defined by Among these three solutions, x 1 and x 3 are positive, and x 2 is negative (not shown). We search for a negative solution because we only consider a compressive state (i.e., σ I < 0), so ′ I = 2 is the correct intersection solution. For details, the first solution x 1 is the intersection of the yield curve and the axis defined by σ II = −σ I for σ I > k t , while the third solution x 3 is the intersection point of the yield curve symmetrical compared to the axis σ I and comes from the square power used between Equations 33 and 34, x 3 ∈ [0, k t ].

Data Availability Statement
The rheological models described in this paper are implemented in the sea ice package of the MIT general circulation model (MITgcm) version 67z (Campin et al., 2021). The model's code and simulations data are available in the Zenodo archive https://doi.org/10.5281/zenodo.6091690 (Ringeisen et al., 2022).