:tocdepth: 3 .. _dynam: Dynamics ======== There are now different rheologies available in the CICE code. The elastic-viscous-plastic (EVP) model represents a modification of the standard viscous-plastic (VP) model for sea ice dynamics :cite:`Hibler79`. The elastic-anisotropic-plastic (EAP) model, on the other hand, explicitly accounts for the observed sub-continuum anisotropy of the sea ice cover :cite:`Wilchinsky06,Weiss09`. If `kdyn` = 1 in the namelist then the EVP rheology is used (module **ice\_dyn\_evp.F90**), while `kdyn` = 2 is associated with the EAP rheology (**ice\_dyn\_eap.F90**). 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 :cite:`Wilchinsky06,Tsamados13`. 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 EVP sea ice dynamics model is thoroughly documented in :cite:`Hunke97`, :cite:`Hunke01`, :cite:`Hunke02` and :cite:`Hunke03` and the EAP dynamics in :cite:`Tsamados13`. 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 :cite:`Hunke99` and :cite:`Tsamados13`. Here we summarize the equations and direct the reader to the above references for details. The numerical implementation in this code release is that of :cite:`Hunke02` and :cite:`Hunke03`, with revisions to the numerical solver as in :cite:`Bouillon13`. The implementation of the EAP sea ice dynamics into CICE is described in detail in :cite:`Tsamados13`. .. _momentum: ******** Momentum ******** The force balance per unit area in the ice pack is given by a two-dimensional momentum equation :cite:`Hibler79`, obtained by integrating the 3D equation through the thickness of the ice in the vertical direction: .. math:: m{\partial {\bf u}\over\partial t} = \nabla\cdot{\bf \sigma} + \vec{\tau}_a+\vec{\tau}_w + \vec{\tau}_b - \hat{k}\times mf{\bf u} - mg\nabla H_\circ, :label: vpmom where :math:`m` is the combined mass of ice and snow per unit area and :math:`\vec{\tau}_a` and :math:`\vec{\tau}_w` are wind and ocean stresses, respectively. The term :math:`\vec{\tau}_b` is a seabed stress (also referred to as basal stress) that represents the grounding of pressure ridges in shallow water :cite:`Lemieux16`. The mechanical properties of the ice are represented by the internal stress tensor :math:`\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 :cite:`Hunke03` and :cite:`Connolley04`. The momentum equation is discretized in time as follows, for the classic EVP approach. First, for clarity, the two components of Equation :eq:`vpmom` are .. math:: \begin{aligned} m{\partial u\over\partial t} &=& {\partial\sigma_{1j}\over\partial x_j} + \tau_{ax} + a_i c_w \rho_w \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\cos\theta - \left(V_w-v\right)\sin\theta\right] -C_bu +mfv - mg{\partial H_\circ\over\partial x}, \\ m{\partial v\over\partial t} &=& {\partial\sigma_{2j}\over\partial x_j} + \tau_{ay} + a_i c_w \rho_w \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\sin\theta + \left(V_w-v\right)\cos\theta\right] -C_bv-mfu - mg{\partial H_\circ\over\partial y}. \end{aligned} In the code, :math:`{\tt vrel}=a_i c_w \rho_w\left|{\bf U}_w - {\bf u}^k\right|` and :math:`C_b=T_b \left( \sqrt{(u^k)^2+(v^k)^2}+u_0 \right)^{-1}`, where :math:`k` denotes the subcycling step. The following equations illustrate the time discretization and define some of the other variables used in the code. .. math:: \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta\ + C_b \right)}_{\tt cca} u^{k+1} - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t_e}u^k, :label: umom .. math:: \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} + \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta + C_b \right)}_{\tt cca}v^{k+1} = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t_e}v^k, :label: vmom and vrel\ :math:`\cdot`\ waterx(y) = taux(y). We solve this system of equations analytically for :math:`u^{k+1}` and :math:`v^{k+1}`. Define .. math:: \hat{u} = F_u + \tau_{ax} - mg{\partial H_\circ\over\partial x} + {\tt vrel} \left(U_w\cos\theta - V_w\sin\theta\right) + {m\over\Delta t_e}u^k :label: cevpuhat .. math:: \hat{v} = F_v + \tau_{ay} - mg{\partial H_\circ\over\partial y} + {\tt vrel} \left(U_w\sin\theta + V_w\cos\theta\right) + {m\over\Delta t_e}v^k, :label: cevpvhat where :math:`{\bf F} = \nabla\cdot\sigma^{k+1}`. Then .. math:: \begin{aligned} \left({m\over\Delta t_e} +{\tt vrel}\cos\theta\ + C_b \right)u^{k+1} - \left(mf + {\tt vrel}\sin\theta\right) v^{k+1} &=& \hat{u} \\ \left(mf + {\tt vrel}\sin\theta\right) u^{k+1} + \left({m\over\Delta t_e} +{\tt vrel}\cos\theta + C_b \right)v^{k+1} &=& \hat{v}.\end{aligned} Solving simultaneously for :math:`u^{k+1}` and :math:`v^{k+1}`, .. math:: \begin{aligned} u^{k+1} = {a \hat{u} + b \hat{v} \over a^2 + b^2} \\ v^{k+1} = {a \hat{v} - b \hat{u} \over a^2 + b^2}, \end{aligned} where .. math:: a = {m\over\Delta t_e} + {\tt vrel}\cos\theta + C_b \\ :label: cevpa .. math:: b = mf + {\tt vrel}\sin\theta. :label: cevpb When the subcycling is finished for each (thermodynamic) time step, the ice–ocean stress must be constructed from `taux(y)` and the terms containing `vrel` on the left hand side of the equations. The Hibler-Bryan form for the ice-ocean stress :cite:`Hibler87` is included in **ice\_dyn\_shared.F90** but is currently commented out, pending further testing. .. _seabed-stress: *************** Seabed stress *************** The parameterization for the seabed stress is described in :cite:`Lemieux16`. The components of the basal seabed stress are :math:`\tau_{bx}=C_bu` and :math:`\tau_{by}=C_bv`, where :math:`C_b` is a coefficient expressed as .. math:: C_b= k_2 \max [0,(h_u - h_{cu})] e^{-\alpha_b * (1 - a_u)} (\sqrt{u^2+v^2}+u_0)^{-1}, \\ :label: Cb where :math:`k_2` determines the maximum seabed stress that can be sustained by the grounded parameterized ridge(s), :math:`u_0` is a small residual velocity and :math:`\alpha_b=20` is a parameter to ensure that the seabed stress quickly drops when the ice concentration is smaller than 1. In the code, :math:`k_2 \max [0,(h_u - h_{cu})] e^{-\alpha_b * (1 - a_u)}` is defined as :math:`T_b`. The quantities :math:`h_u`, :math:`a_{u}` and :math:`h_{cu}` are calculated at the 'u' point based on local ice conditions (surrounding tracer points). They are respectively given by .. math:: h_u=\max[v_i(i,j),v_i(i+1,j),v_i(i,j+1),v_i(i+1,j+1)], \\ :label: hu .. math:: a_u=\max[a_i(i,j),a_i(i+1,j),a_i(i,j+1),a_i(i+1,j+1)]. \\ :label: au .. math:: h_{cu}=a_u h_{wu} / k_1, \\ :label: hcu where the :math:`a_i` and :math:`v_i` are the total ice concentrations and ice volumes around the :math:`u` point :math:`i,j` and :math:`k_1` is a parameter that defines the critical ice thickness :math:`h_{cu}` at which the parameterized ridge(s) reaches the seafloor for a water depth :math:`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 :math:`C_b` in equation :eq:`Cb`, the seabed stress components are non-zero only when :math:`h_u > h_{cu}`. The maximum seabed stress depends on the weigth of the ridge above hydrostatic balance and the value of :math:`k_2`. It is, however, the parameter :math:`k_1` that has the most notable impact on the simulated extent of landfast ice. The value of :math:`k_1` can be changed at runtime using the namelist variable ``k1``. The grounding scheme can be turned on or off using the namelist logical basalstress. Note that the user must provide a bathymetry field for using this grounding scheme. It is suggested to have a bathymetry field with water depths larger than 5 m that represents well shallow water regions such as the Laptev Sea and the East Siberian Sea. To prevent unrealistic grounding, :math:`T_b` is set to zero when :math:`h_{wu}` is larger than 30 m. This maximum value is chosen based on observations of large keels in the Arctic Ocean :cite:`Amundrud04`. .. _internal-stress: *************** Internal stress *************** For convenience we formulate the stress tensor :math:`\bf \sigma` in terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, :math:`\sigma_2=\sigma_{11}-\sigma_{22}`, and introduce the divergence, :math:`D_D`, and the horizontal tension and shearing strain rates, :math:`D_T` and :math:`D_S` respectively. CICE now outputs the internal ice pressure which is an important field to support navigation in ice-infested water. The internal ice pressure :math:`(sigP)` is the average of the normal stresses multiplied by :math:`-1` and is therefore simply equal to :math:`-\sigma_1/2`. *Elastic-Viscous-Plastic* In the EVP model the internal stress tensor is determined from a regularized version of the VP constitutive law. Following the approach of :cite:`Konig10` (see also :cite:`Lemieux16`), the elliptical yield curve can be modified such that the ice has isotropic tensile strength. The tensile strength :math:`T_p` is expressed as a fraction of the ice strength :math:`P`, that is :math:`T_p=k_t P` where :math:`k_t` should be set to a value between 0 and 1 (this can be changed at runtime with the namelist parameter ``Ktens``). The constitutive law is therefore .. math:: {1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} + {P_R(1-k_t)\over 2\zeta} = D_D, \\ :label: sig1 .. math:: {1\over E}{\partial\sigma_2\over\partial t} + {\sigma_2\over 2\eta} = D_T, :label: sig2 .. math:: {1\over E}{\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 2\eta} = {1\over 2}D_S, :label: sig12 where .. math:: D_D = \dot{\epsilon}_{11} + \dot{\epsilon}_{22}, .. math:: D_T = \dot{\epsilon}_{11} - \dot{\epsilon}_{22}, .. math:: D_S = 2\dot{\epsilon}_{12}, .. math:: \dot{\epsilon}_{ij} = {1\over 2}\left({{\partial u_i}\over{\partial x_j}} + {{\partial u_j}\over{\partial x_i}}\right), .. math:: \zeta = {P(1+k_t)\over 2\Delta}, .. math:: \eta = {P(1+k_t)\over {2\Delta e^2}}, .. math:: \Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2}, and :math:`P_R` is a “replacement pressure” (see :cite:`Geiger98`, for example), which serves to prevent residual ice motion due to spatial variations of :math:`P` when the rates of strain are exactly zero. The ice strength :math:`P` is a function of the ice thickness and concentration as it is described in the `Icepack Documentation `_. The parameteter :math:`e` is the ratio of the major and minor axes of the elliptical yield curve, also called the ellipse aspect ratio. It can be changed using the namelist parameter ``e_ratio``. Viscosities are updated during the subcycling, so that the entire dynamics component is subcycled within the time step, and the elastic parameter :math:`E` is defined in terms of a damping timescale :math:`T` for elastic waves, :math:`\Delta t_e < T < \Delta t`, as .. math:: E = {\zeta\over T}, where :math:`T=E_\circ\Delta t` and :math:`E_\circ` (eyc) is a tunable parameter less than one. Including the modification proposed by :cite:`Bouillon13` for equations :eq:`sig2` and :eq:`sig12` in order to improve numerical convergence, the stress equations become .. math:: \begin{aligned} {\partial\sigma_1\over\partial t} + {\sigma_1\over 2T} + {P_R(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta} D_D, \\ {\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {P(1+k_t)\over 2Te^2\Delta} D_T,\\ {\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 2T} &=& {P(1+k_t)\over 4Te^2\Delta}D_S.\end{aligned} Once discretized in time, these last three equations are written as .. math:: \begin{aligned} {(\sigma_1^{k+1}-\sigma_1^{k})\over\Delta t_e} + {\sigma_1^{k+1}\over 2T} + {P_R^k(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta^k} D_D^k, \\ {(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {P(1+k_t)\over 2Te^2\Delta^k} D_T^k,\\ {(\sigma_{12}^{k+1}-\sigma_{12}^{k})\over\Delta t_e} + {\sigma_{12}^{k+1}\over 2T} &=& {P(1+k_t)\over 4Te^2\Delta^k}D_S^k,\end{aligned} :label: sigdisc where :math:`k` denotes again the subcycling step. All coefficients on the left-hand side are constant except for :math:`P_R`. This modification compensates for the decreased efficiency of including the viscosity terms in the subcycling. (Note that the viscosities do not appear explicitly.) Choices of the parameters used to define :math:`E`, :math:`T` and :math:`\Delta t_e` are discussed in Sections :ref:`revp` and :ref:`parameters`. The bilinear discretization used for the stress terms :math:`\partial\sigma_{ij}/\partial x_j` in the momentum equation is now used, which enabled 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 :cite:`Hunke02`. *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 :ref:`fig-EAP`). 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 :math:`\mathbf r` in :ref:`fig-EAP`), 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 :math:`\mathbf{A}` defined by .. math:: {\mathbf A}=\int_{\mathbb{S}}\vartheta(\mathbf r)\mathbf r\mathbf r d\mathbf r\label{structuretensor}. where :math:`\mathbb{S}` is a unit-radius circle; **A** is a unit trace, 2\ :math:`\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 :math:`\vartheta(\mathbf r )` to be Gaussian, :math:`\vartheta(z)=\omega_{1}\exp(-\omega_{2}z^{2})`, where :math:`z` is the ice floe inclination with respect to the axis :math:`x_{1}` of preferential alignment of ice floes (see :ref:`fig-EAP`), :math:`\vartheta(z)` is periodic with period :math:`\pi`, and the positive coefficients :math:`\omega_{1}` and :math:`\omega_{2}` are calculated to ensure normalization of :math:`\vartheta(z)`, i.e. :math:`\int_{0}^{2\pi}\vartheta(z)dz=1`. The ratio of the principal components of :math:`\mathbf{A}`, :math:`A_{1}/A_{2}`, are derived from the phenomenological evolution equation for the structure tensor :math:`\mathbf A`, .. math:: \frac{D\mathbf{A}}{D t}=\mathbf{F}_{iso}(\mathbf{A})+\mathbf{F}_{frac}(\mathbf{A},\boldsymbol\sigma), :label: evolutionA where :math:`t` is the time, and :math:`D/Dt` is the co-rotational time derivative accounting for advection and rigid body rotation (:math:`D\mathbf A/Dt = d\mathbf A/dt -\mathbf W \cdot \mathbf A -\mathbf A \cdot \mathbf W^{T}`) with :math:`\mathbf W` being the vorticity tensor. :math:`\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. :math:`\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 :cite:`Wilchinsky06`, based on laboratory experiments by :cite:`Schulson01` 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 :math:`\sigma_{1}` and :math:`\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, (:math:`\sigma_{1}/\sigma_{2}