2.5. Dynamics¶
The force balance per unit area in the ice pack is given by a two-dimensional momentum equation [15], obtained by integrating the 3D equation through the thickness of the ice in the vertical direction:
where \(m\) is the combined mass of ice and snow per unit area and \(\vec{\tau}_a\) and \(\vec{\tau}_w\) are wind and ocean stresses, respectively. The term \(\vec{\tau}_b\) is a seabed stress (also referred to as basal stress) that represents the grounding of pressure ridges in shallow water [34]. The mechanical properties of the ice are represented by the internal stress tensor \(\sigma_{ij}\). The other two terms on the right hand side are stresses due to Coriolis effects and the sea surface slope. The parameterization for the wind and ice–ocean stress terms must contain the ice concentration as a multiplicative factor to be consistent with the formal theory of free drift in low ice concentration regions. A careful explanation of the issue and its continuum solution is provided in [23] and [8].
For clarity, the two components of Equation (1) are
On the B grid, the equations above are solved at the U point for the collocated u and v components (see figure Schematic of CICE B-grid.). On the C grid, however, the two components are not collocated: the u component is at the E point while the v component is at the N point.
The B grid spatial discretization is based on a variational method described in [21] and [22]. A bilinear discretization is used for the stress terms \(\partial\sigma_{ij}/\partial x_j\), which enables the discrete equations to be derived from the continuous equations written in curvilinear coordinates. In this manner, metric terms associated with the curvature of the grid are incorporated into the discretization explicitly. Details pertaining to the spatial discretization are found in [22]
On the C grid, however, a finite difference approach is used for the spatial discretization. The C grid discretization is based on [7], [6] and [29].
There are different approaches in the CICE code for representing sea ice
rheology and for solving the sea ice momentum equation: the viscous-plastic (VP) rheology [15] with an implicit method,
the elastic-viscous-plastic (EVP) [21] model which represents a modification of the
VP model, the revised EVP (rEVP) approach [35][6] and the elastic-anisotropic-plastic (EAP) model which explicitly accounts for the sub-continuum
anisotropy of the sea ice cover [65][63]. If
kdyn = 1 in the namelist then the EVP model is used (module
ice_dyn_evp.F90), while kdyn = 2 is associated with the EAP
model (ice_dyn_eap.F90), and kdyn = 3 is associated with the
VP model (ice_dyn_vp.F90). The rEVP approach can be used by setting kdyn = 1 and revised_evp = true in the namelist.
At times scales associated with the wind forcing, the EVP model reduces to the VP model while the EAP model reduces to the anisotropic rheology described in detail in [65][60]. At shorter time scales the adjustment process takes place in both models by a numerically more efficient elastic wave mechanism. While retaining the essential physics, this elastic wave modification leads to a fully explicit numerical scheme which greatly improves the model’s computational efficiency. The rEVP is also a fully explicit scheme which by construction should lead to the VP solution.
The EVP sea ice dynamics model is thoroughly documented in [21], [20], [22] and [23] and the EAP dynamics in [60]. Simulation results and performance of the EVP and EAP models have been compared with the VP model and with each other in realistic simulations of the Arctic respectively in [25] and [60].
The EVP numerical implementation in this code release is that of [22] and [23], with revisions to the numerical solver as in [6]. Details about the rEVP solver can be found in [35], [6], [28] and [30]. The implementation of the EAP sea ice dynamics into CICE is described in detail in [60].
The VP solver implementation mostly follows [33], with
FGMRES [53] as the linear solver and GMRES as the preconditioner.
Note that the VP solver has not yet been tested on the tx1 grid.
The EVP, rEVP, EAP and VP approaches are all available with the B grid. However, at the moment, only the EVP and rEVP schemes are possible with the C grid.
Here we summarize the equations and direct the reader to the above references for details.
2.5.1. Momentum time stepping¶
2.5.1.1. EVP time discretization and solution¶
The momentum equation is discretized in time as follows, for the classic EVP approach. In the code, \({\tt vrel}=a_i c_w \rho_w\left|{\bf U}_w - {\bf u}^k\right|\) and \(C_b=T_b \left( \sqrt{(u^k)^2+(v^k)^2}+u_0 \right)^{-1}\), where \(k\) denotes the subcycling step. The following equations illustrate the time discretization and define some of the other variables used in the code.
where \({\tt vrel}\ \cdot\ {\tt waterx(y)}= {\tt taux(y)}\) and the definitions of \(u^{l}\) and \(v^{l}\) vary depending on the grid.
As \(u\) and \(v\) are collocated on the B grid, \(u^{l}\) and \(v^{l}\) are respectively \(u^{k+1}\) and \(v^{k+1}\) such that this system of equations can be solved as follows. Define
where \({\bf F} = \nabla\cdot\sigma^{k+1}\). Then
Solving simultaneously for \(u^{k+1}\) and \(v^{k+1}\),
where
Note that the time discretization and solution method for the EAP is exactly the same as for the B grid EVP. More details on the EAP model are given in Section Elastic-Anisotropic-Plastic.
However, on the C grid, \(u\) and \(v\) are not collocated. When solving the \(u\) momentum equation for \(u^{k+1}\) (at the E point), \(v^{l}=v^{k}_{int}\) where \(v^{k}_{int}\) is \(v^{k}\) from the surrounding N points interpolated to the E point. The same approach is used for the \(v\) momentum equation. With this explicit treatment of the off-diagonal terms [29], \(u^{k+1}\) and \(v^{k+1}\) are obtained by solving
2.5.1.2. Revised EVP time discretization and solution¶
The revised EVP approach is based on a pseudo-time iterative scheme [35], [6], [28]. By construction, the revised EVP approach should lead to the VP solution (given the right numerical parameters and a sufficiently large number of iterations). To do so, the inertial term is formulated such that it matches the backward Euler approach of implicit solvers and there is an additional term for the pseudo-time iteration. Hence, with the revised approach, the discretized momentum equations (3) and (4) become
where \(\beta^*\) is a numerical parameter and \(u^n, v^n\) are the components of the previous time level solution. With \(\beta=\beta^* \Delta t \left( m \Delta t_e \right)^{-1}\) [6], these equations can be written as
At this point, the solutions \(u^{k+1}\) and \(v^{k+1}\) for the B or the C grids are obtained in the same manner as for the standard EVP approach (see Section EVP time discretization and solution for details).
2.5.1.3. Implicit (VP) time discretization and solution¶
In the VP approach, equation (2) is discretized implicitly using a Backward Euler approach, and stresses are not computed explicitly:
where \(r = (r_x,r_y)\) contains all terms that do not depend on the velocities \(u^n, v^n\) (namely the sea surface tilt and the wind stress). As the water drag, seabed stress and rheology term depend on the velocity field, the only unknowns in equation (13) are \(u^n\) and \(v^n\).
Once discretized in space, equation (13) leads to a system of \(N\) nonlinear equations with \(N\) unknowns that can be concisely written as
where \(\mathbf{A}\) is an \(N\times N\) matrix and \(\mathbf{u}\) and \(\mathbf{b}\) are vectors of size \(N\). Note that we have dropped the time level index \(n\). The vector \(\mathbf{u}\) is formed by stacking first the \(u\) components, followed by the \(v\) components of the discretized ice velocity. The vector \(\mathbf{b}\) is a function of the velocity vector \(\mathbf{u}\) because of the water and seabed stress terms as well as parts of the rheology term that depend non-linearly on \(\mathbf{u}\).
The nonlinear system (14) is solved using a Picard iteration method. Starting from a previous iterate \(\mathbf{u}_{k-1}\), the nonlinear system is linearized by substituting \(\mathbf{u}_{k-1}\) in the expression of the matrix \(\mathbf{A}\) and the vector \(\mathbf{b}\):
The resulting linear system is solved using the Flexible Generalized Minimum RESidual (FGMRES, [53]) method and this process is repeated iteratively.
The maximum number of Picard iterations can be set using the namelist flag maxits_nonlin.
The relative tolerance for the Picard solver can be set using the namelist flag reltol_nonlin.
The Picard iterative process stops when \(\left\lVert \mathbf{u}_{k} \right\rVert_2 < {\tt reltol\_nonlin} \cdot \left\lVert\mathbf{u}_{0}\right\rVert_2\) or when maxits_nonlin is reached.
Parameters for the FGMRES linear solver and the preconditioner can be controlled using additional namelist flags (see dynamics_nml).
2.5.2. Surface stress terms¶
The formulation for the wind stress is described in Icepack Documentation. Below, some details about the ice-ocean stress and the seabed stress are given.
2.5.2.1. Ice-Ocean stress¶
At the end of each (thermodynamic) time step, the ice–ocean stress must be constructed from \({\tt taux(y)}\) and the terms containing \({\tt vrel}\) on the left hand side of the equations.
The Hibler-Bryan form for the ice-ocean stress [17] is included in ice_dyn_shared.F90 but is currently commented out, pending further testing.
2.5.2.2. Seabed stress¶
CICE includes two options for calculating the seabed stress,
i.e. the term in the momentum equation that represents the interaction
between grounded ice keels and the seabed. The seabed stress can be
activated by setting seabed_stress to true in the namelist. The seabed stress (or basal
stress) parameterization of [34] is chosen if seabed_stress_method = LKD while the approach based on the probability of contact between the ice and the seabed is used if seabed_stress_method = probabilistic.
For both parameterizations, the components of the seabed stress are expressed as \(\tau_{bx}=C_bu\) and \(\tau_{by}=C_bv\), where \(C_b\) is a seabed stress coefficient.
The two parameterizations differ in their calculation of the \(C_b\) coefficients.
Note that the user must provide a bathymetry field for using these grounding schemes. It is suggested to have a bathymetry field with water depths larger than 5 m that represents well shallow water (less than 30 m) regions such as the Laptev Sea and the East Siberian Sea.
Seabed stress based on linear keel draft (LKD)
This parameterization for the seabed stress is described in [34]. It assumes that the largest keel draft varies linearly with the mean thickness in a grid cell (i.e. sea ice volume). The \(C_b\) coefficients are expressed as
where \(k_2\) determines the maximum seabed stress that can be sustained by the grounded parameterized ridge(s), \(u_0\) is a small residual velocity and \(\alpha_b\) is a parameter to ensure that the seabed stress quickly drops when the ice concentration is smaller than 1. In the code, \(k_2 \max [0,(h - h_{c})] e^{-\alpha_b * (1 - a)}\) is defined as \(T_b\).
On the B grid, the quantities \(h\), \(a\) and \(h_{c}\) are calculated at the U point and are referred to as \(h_u\), \(a_{u}\) and \(h_{cu}\). They are respectively given by
where the \(a_i\) and \(v_i\) are the total ice concentrations and ice volumes around the U point \(i,j\) and \(k_1\) is a parameter that defines the critical ice thickness \(h_{cu}\) at which the parameterized ridge(s) reaches the seafloor for a water depth \(h_{wu}=\min[h_w(i,j),h_w(i+1,j),h_w(i,j+1),h_w(i+1,j+1)]\). Given the formulation of \(C_b\) in equation (16), the seabed stress components are non-zero only when \(h_u > h_{cu}\).
As \(u\) and \(v\) are not collocated on the C grid, \(T_b\) is calculated at E and N points. For example, at the E point, \(h_e\), \(a_{e}\) and \(h_{ce}\) are respectively
where \(h_{we}=\min[h_w(i,j),h_w(i+1,j)]\). Similar calculations are done at the N points.
To prevent unrealistic grounding, \(T_b\) is set to zero when \(h_{wu}\) is larger than 30 m (same idea on the C grid depending on \(h_{we}\) and \(h_{wn}\)). This maximum value is chosen based on observations of large keels in the Arctic Ocean [1].
The maximum seabed stress depends on the weight of the ridge
above hydrostatic balance and the value of \(k_2\). It is, however, the parameter \(k_1\) that has the most notable impact on the simulated extent of landfast ice.
The value of \(k_1\) can be changed at runtime using the namelist variable k1.
Seabed stress based on probabilistic approach
This more sophisticated grounding parameterization computes the seabed stress based on the probability of contact between the ice thickness distribution (ITD) and the seabed [11]. Multi-thickness category models such as CICE typically use a few thickness categories (5-10). This crude representation of the ITD does not resolve the tail of the ITD, which is crucial for grounding events.
To represent the tail of the distribution, the simulated ITD is converted to a positively skewed probability function \(f(x)\) with \(x\) the sea ice thickness. The mean and variance are set equal to the ones of the original ITD. A log-normal distribution is used for \(f(x)\).
It is assumed that the bathymetry \(y\) (at the ‘t’ point) follows a normal distribution \(b(y)\). The mean of \(b(y)\) comes from the user’s bathymetry field and the standard deviation \(\sigma_b\) is currently fixed to 2.5 m. Two possible improvements would be to specify a distribution based on high resolution bathymetry data and to take into account variations of the water depth due to changes in the sea surface height.
Assuming hydrostatic balance and neglecting the impact of snow, the draft of floating ice of thickness \(x\) is \(D(x)=\rho_i x / \rho_w\) where \(\rho_i\) is the sea ice density. Hence, the probability of contact (\(P_c\)) between the ITD and the seabed is given by
\(T_b\) is first calculated at the T point (referred to as \(T_{bt}\)). \(T_{bt}\) depends on the weight of the ridge in excess of hydrostatic balance. The parameterization first calculates
and then obtains \(T_{bt}\) by multiplying \(T_{bt}^*\) by \(e^{-\alpha_b * (1 - a_i)}\) (similar to what is done for seabed_stress_method = LKD).
To calculate \(T_{bt}^*\) in equation (23), \(f(x)\) and \(b(y)\) are discretized using many small categories (100). \(f(x)\) is discretized between 0 and 50 m while \(b(y)\) is truncated at plus and minus three \(\sigma_b\). \(f(x)\) is also modified by setting it to zero after a certain percentile of the log-normal distribution. This percentile, which is currently set to 99.7%, notably affects the simulation of landfast ice and is used as a tuning parameter. Its impact is similar to the one of the parameter \(k_1\) for the LKD method.
On the B grid, \(T_b\) at the U point is calculated from the T point values around it according to
Following again the LKD method, the seabed stress coefficients are finally expressed as
On the C grid, \(T_b\) is needs to be calculated at the E and N points. \(T_{be}\) and \(T_{bn}\) are respectively given by
The \(C_{b}\) are different at the E and N points and are respectively \(T_{be} (\sqrt{u^2+v^2_{int}}+u_0)^{-1}\) and \(T_{bn} (\sqrt{u^2_{int} + v^2}+u_0)^{-1}\) where \(v_{int}\) (\(u_{int}\)) is \(v\) ( \(u\)) interpolated to the E (N) point.
2.5.3. Rheology¶
For convenience we formulate the stress tensor \(\bf \sigma\) in
terms of \(\sigma_1=\sigma_{11}+\sigma_{22}\) (stressp),
\(\sigma_2=\sigma_{11}-\sigma_{22}\) (stressm), and introduce the
divergence, \(D_D\), and the horizontal tension and shearing
strain rates, \(D_T\) and \(D_S\) respectively:
where
Note that \(\sigma_1\) and \(\sigma_2\) are not to be confused with the normalized principal stresses,
\(\sigma_{n,1}\) and \(\sigma_{n,2}\) (sig1 and sig2), which are defined as:
where \(P\) is the ice strength.
In addition to the normalized principal stresses, CICE can output the internal ice pressure which is an important field to support navigation in ice-infested water.
The internal ice pressure (sigP) is the average of the normal stresses (\(\sigma_{11}\), \(\sigma_{22}\)) multiplied by \(-1\) and
is therefore simply equal to \(-\sigma_1/2\).
2.5.3.1. Viscous-Plastic¶
The VP constitutive law is given by
where \(\eta\) and \(\zeta\) are the bulk and shear viscosities and \(P_R\) is a “replacement pressure” (see [13], for example), which serves to prevent residual ice motion due to spatial variations of the ice strength \(P\) when the strain rates are exactly zero.
An elliptical yield curve is used, with the viscosities given by
where
When the deformation \(\Delta\) tends toward zero, the viscosities tend toward infinity. To avoid this issue, \(\Delta\) needs to be limited and is replaced by \(\Delta^*\) in equation (29). Two methods for limiting \(\Delta\) (or for capping the viscosities) are available in the code. If the namelist parameter capping_method is set to max, \(\Delta^*=max(\Delta, \Delta_{min})\) [15] while with capping_method set to sum, the smoother formulation \(\Delta^*=(\Delta + \Delta_{min})\) of [31] is used.
The ice strength \(P\) is a function of the ice thickness distribution as described in the Icepack Documentation.
Two other modifications to the standard VP rheology of [15] are available.
First, following the approach of [3] (see also [34]), the
elliptical yield curve can be modified such that the ice has isotropic tensile strength.
The tensile strength is expressed as a fraction of \(P\), that is \(k_t P\)
where \(k_t\) should be set to a value between 0 and 1 (this can
be changed at runtime with the namelist parameter Ktens).
Second, while \(e_f\) is the ratio of the major and minor axes of the elliptical yield curve, the parameter
\(e_g\) characterizes the plastic potential, i.e. another ellipse that decouples the flow rule from the
yield curve ([46]). \(e_f\) and \(e_g\) are respectively called e_yieldcurve and e_plasticpot in the code and
can be set in the namelist. The plastic potential can lead to more realistic fracture angles between linear kinematic features. [46] suggest to set \(e_f\) to a value larger than 1 and to have \(e_g < e_f\).
By default, the namelist parameters are set to \(e_f=e_g=2\) and \(k_t=0\) which correspond to the standard VP rheology.
There are four options in the code for solving the sea ice momentum equation with a VP formulation: the standard EVP approach, a 1d EVP solver, the revised EVP approach and an implicit Picard solver. The choice of the capping method for the viscosities and the modifications to the yield curve and to the flow rule described above are available for these four different solution methods. Note that only the EVP and revised EVP methods are currently available if one chooses the C grid.
2.5.3.2. Elastic-Viscous-Plastic¶
In the EVP model the internal stress tensor is determined from a regularized version of the VP constitutive law (28). The constitutive law is therefore
Viscosities are updated during the subcycling, so that the entire dynamics component is subcycled within the time step, and the elastic parameter \(E\) is defined in terms of a damping timescale \(T\) for elastic waves, \(\Delta t_e < T < \Delta t\), as
where \(T=E_\circ\Delta t\) and \(E_\circ\) (elasticDamp) is a tunable parameter less than one. Including the modification proposed by [6] for equations (33) and (34) in order to improve numerical convergence, the stress equations become
Once discretized in time, these last three equations are written as
where \(k\) denotes again the subcycling step. All coefficients on the left-hand side are constant except for \(P_R\). This modification compensates for the decreased efficiency of including the viscosity terms in the subcycling. Choices of the parameters used to define \(E\), \(T\) and \(\Delta t_e\) are discussed in Sections Revised EVP approach and Choosing an appropriate time step.
On the B grid, the stresses \(\sigma_{1}\), \(\sigma_{2}\) and \(\sigma_{12}\) are collocated at the U point. To calculate these stresses, the viscosities \(\zeta\) and \(\eta\) and the replacement pressure \(P_R\) are also defined at the U point.
However, on the C grid, \(\sigma_{1}\) and \(\sigma_{2}\) are collocated at the T point while \(\sigma_{12}\) is defined at the U point. During a subcycling step, \(\zeta\), \(\eta\) and \(P_R\) are first calculated at the T point. To do so, \(\Delta\) given by equation (31) is calculated following the approach of [6] (see also [29] for details). With this approach, \(D_S^2\) at the T point is obtained by calculating \(D_S^2\) at the U points and interpolating these values to the T point. As \(\sigma_{12}\) is calculated at the U point, \(\eta\) also needs to be computed as these locations. If visc_method in the namelist is set to avg_zeta (the default value), \(\eta\) at the U point is obtained by interpolating T point values to this location. This corresponds to the approach used by [6] and the one associated with the C1 configuration of [29]. On the other hand, if visc_method = avg_strength, the strength \(P\) calculated at T points is interpolated to the U point and \(\Delta\) is calculated at the U point in order to obtain \(\eta\) following equations (29) and (30). This latter approach is the one used in the C2 configuration of [29].
2.5.3.3. 1d EVP solver¶
The standard EVP solver iterates hundreds of times, where each iteration includes a communication through MPI and a limited number of calculations. This limits how much the solver can be optimized as the speed is primarily determined by the communication. The 1d EVP solver avoids the communication by utilizing shared memory, which removes the requirement for calls to the MPI communicator. As a consequence of this the potential scalability of the code is improved. The performance is best on shared memory but the solver is also functional on MPI and hybrid MPI/OpenMP setups as it will run on the master processor alone.
The scalability of geophysical models is in general terms limited by the memory usage. In order to optimize this the 1d EVP solver solves the same equations that are outlined in the section Elastic-Viscous-Plastic but it transforms all matrices to vectors (1d matrices) as this compiles better with the computer hardware. The vectorization and the contiguous placement of arrays in the memory makes it easier for the compiler to optimize the code and pass pointers instead of copying the vectors. The 1d solver is not supported for tripole grids and the code will abort if this combination is attempted.
2.5.3.4. Revised EVP approach¶
Introducing the numerical parameter \(\alpha=2T \Delta t_e ^{-1}\) [6], the stress equations in (35) become
where as opposed to the classic EVP, the second term in each equation is at iteration \(k\) [6]. Also, contrary to the classic EVP,
\(\Delta t_e\) times the number of subcycles (or iterations) does not need to be equal to the advective time step \(\Delta t\).
Finally, as with the classic EVP approach, the stresses are initialized using the previous time level values.
The revised EVP is activated by setting the namelist parameter revised_evp = true.
In the code \(\alpha\) is arlx and \(\beta\) is brlx (introduced in Section Revised EVP time discretization and solution). The values of arlx and brlx can be set in the namelist.
It is recommended to use large values of these parameters and to set \(\alpha=\beta\) [28].
2.5.3.5. Elastic-Anisotropic-Plastic¶
In the EAP model the internal stress tensor is related to the geometrical properties and orientation of underlying virtual diamond shaped floes (see Diamond-shaped floes). In contrast to the isotropic EVP rheology, the anisotropic plastic yield curve within the EAP rheology depends on the relative orientation of the diamond shaped floes (unit vector \(\mathbf r\) in Diamond-shaped floes), with respect to the principal direction of the deformation rate (not shown). Local anisotropy of the sea ice cover is accounted for by an additional prognostic variable, the structure tensor \(\mathbf{A}\) defined by
where \(\mathbb{S}\) is a unit-radius circle; A is a unit trace, 2\(\times\)2 matrix. From now on we shall describe the orientational distribution of floes using the structure tensor. For simplicity we take the probability density function \(\vartheta(\mathbf r )\) to be Gaussian, \(\vartheta(z)=\omega_{1}\exp(-\omega_{2}z^{2})\), where \(z\) is the ice floe inclination with respect to the axis \(x_{1}\) of preferential alignment of ice floes (see Diamond-shaped floes), \(\vartheta(z)\) is periodic with period \(\pi\), and the positive coefficients \(\omega_{1}\) and \(\omega_{2}\) are calculated to ensure normalization of \(\vartheta(z)\), i.e. \(\int_{0}^{2\pi}\vartheta(z)dz=1\). The ratio of the principal components of \(\mathbf{A}\), \(A_{1}/A_{2}\), are derived from the phenomenological evolution equation for the structure tensor \(\mathbf A\),
where \(t\) is the time, and \(D/Dt\) is the co-rotational time derivative accounting for advection and rigid body rotation (\(D\mathbf A/Dt = d\mathbf A/dt -\mathbf W \cdot \mathbf A -\mathbf A \cdot \mathbf W^{T}\)) with \(\mathbf W\) being the vorticity tensor. \(\mathbf F_{iso}\) is a function that accounts for a variety of processes (thermal cracking, melting, freezing together of floes) that contribute to a more isotropic nature to the ice cover. \(\mathbf F_{frac}\) is a function determining the ice floe re-orientation due to fracture, and explicitly depends upon sea ice stress (but not its magnitude). Following [65], based on laboratory experiments by [54] we consider four failure mechanisms for the Arctic sea ice cover. These are determined by the ratio of the principal values of the sea ice stress \(\sigma_{1}\) and \(\sigma_{2}\): (i) under biaxial tension, fractures form across the perpendicular principal axes and therefore counteract any apparent redistribution of the floe orientation; (ii) if only one of the principal stresses is compressive, failure occurs through axial splitting along the compression direction; (iii) under biaxial compression with a low confinement ratio, (\(\sigma_{1}/\sigma_{2}<R\)), sea ice fails Coulombically through formation of slip lines delineating new ice floes oriented along the largest compressive stress; and finally (iv) under biaxial compression with a large confinement ratio, (\(\sigma_{1}/\sigma_{2}\ge R\)), the ice is expected to fail along both principal directions so that the cumulative directional effect balances to zero.
Diamond-shaped floes¶
Figure Diamond-shaped floes shows geometry of interlocking diamond-shaped floes (taken from [65]). \(\phi\) is half of the acute angle of the diamonds. \(L\) is the edge length. \(\boldsymbol n_{1}\), \(\boldsymbol n_{2}\) and \(\boldsymbol\tau_{1}\), \(\boldsymbol\tau_{2}\) are respectively the normal and tangential unit vectors along the diamond edges. \(\mathbf v=L\boldsymbol\tau_{2}\cdot\dot{\boldsymbol\epsilon}\) is the relative velocity between the two floes connected by the vector \(L \boldsymbol \tau_{2}\). \(\mathbf r\) is the unit vector along the main diagonal of the diamond. Note that the diamonds illustrated here represent one possible realisation of all possible orientations. The angle \(z\) represents the rotation of the diamonds’ main axis relative to their preferential orientation along the axis \(x_1\).
The new anisotropic rheology requires solving the evolution Equation (36) for the structure tensor in addition to the momentum and stress equations. The evolution equation for \(\mathbf{A}\) is solved within the EVP subcycling loop, and consistently with the momentum and stress evolution equations, we neglect the advection term for the structure tensor. Equation (36) then reduces to the system of two equations:
where the first terms on the right hand side correspond to the isotropic contribution, \(F_{iso}\), and \(M_{11}\) and \(M_{12}\) are the components of the term \(F_{frac}\) in Equation (36) that are given in [65] and [60]. These evolution equations are discretized semi-implicitly in time. The degree of anisotropy is measured by the largest eigenvalue (\(A_{1}\)) of this tensor (\(A_{2}=1-A_{1}\)). \(A_{1}=1\) corresponds to perfectly aligned floes and \(A_{1}=0.5\) to a uniform distribution of floe orientation. Note that while we have specified the aspect ratio of the diamond floes, through prescribing \(\phi\), we make no assumption about the size of the diamonds so that formally the theory is scale invariant.
As described in greater detail in [65], the internal ice stress for a single orientation of the ice floes can be calculated explicitly and decomposed, for an average ice thickness \(h\), into its ridging (r) and sliding (s) contributions
where \(P_{r}\) and \(P_{s}\) are the ridging and sliding strengths and the ridging and sliding stresses are functions of the angle \(\theta= \arctan(\dot\epsilon_{II}/\dot\epsilon_{I})\), the angle \(y\) between the major principal axis of the strain rate tensor (not shown) and the structure tensor (\(x_1\) axis in Diamond-shaped floes, and the angle \(z\) defined in Diamond-shaped floes. In the stress expressions above the underlying floes are assumed parallel, but in a continuum-scale sea ice region the floes can possess different orientations in different places and we take the mean sea ice stress over a collection of floes to be given by the average
where we have introduced the friction parameter \(k=P_{s}/P_{r}\) and where we identify the ridging ice strength \(P_{r}(h)\) with the strength \(P\) described in section 1 and used within the EVP framework.
As is the case for the EVP rheology, elasticity is included in the EAP description not to describe any physical effect, but to make use of the efficient, explicit numerical algorithm used to solve the full sea ice momentum balance. We use the analogous EAP stress equations,
where the anisotropic stress \(\boldsymbol\sigma^{EAP}\) is defined in a look-up table for the current values of strain rate and structure tensor. The look-up table is constructed by computing the stress (normalized by the strength) from Equations (39)–(41) for discrete values of the largest eigenvalue of the structure tensor, \(\frac{1}{2}\le A_{1}\le 1\), the angle \(0\le\theta\le2\pi\), and the angle \(-\pi/2\le y\le\pi/2\) between the major principal axis of the strain rate tensor and the structure tensor [60]. The updated stress, after the elastic relaxation, is then passed to the momentum equation and the sea ice velocities are updated in the usual manner within the subcycling loop of the EVP rheology. The structure tensor evolution equations are solved implicitly at the same frequency, \(\Delta t_{e}\), as the ice velocities and internal stresses. Finally, to be coherent with our new rheology we compute the area loss rate due to ridging as \(\vert\dot{\boldsymbol\epsilon}\vert\alpha_{r}(\theta)\), with \(\alpha_r(\theta)\) and \(\alpha_s(\theta)\) given by [64],
Both ridging rate and sea ice strength are computed in the outer loop of the dynamics.