:tocdepth: 3 .. _dynam: Dynamics ======== There are now different dynamical solvers 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 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**). 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`. The EVP 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`. The VP solver implementation mostly follows :cite:`Lemieux07`, with FGMRES :cite:`Saad93` as the linear solver and GMRES as the preconditioner. Here we summarize the equations and direct the reader to the above references for details. .. _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`. 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} :label: momsys A bilinear discretization is used for the stress terms :math:`\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 :cite:`Hunke02`. .. _evp-momentum: Elastic-Viscous-Plastic ~~~~~~~~~~~~~~~~~~~~~~~ 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 :math:`{\tt vrel}\ \cdot\ {\tt waterx(y)}= {\tt 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 .. _vp-momentum: Viscous-Plastic ~~~~~~~~~~~~~~~ In the VP approach, equation :eq:`momsys` is discretized implicitly using a Backward Euler approach, and stresses are not computed explicitly: .. math:: \begin{align} m\frac{(u^{n}-u^{n-1})}{\Delta t} &= \frac{\partial \sigma_{1j}^n}{\partial x_j} - \tau_{w,x}^n + \tau_{b,x}^n + mfv^n + r_{x}^n, \\ m\frac{(v^{n}-v^{n-1})}{\Delta t} &= \frac{\partial \sigma^{n} _{2j}}{\partial x_j } - \tau_{w,y}^n + \tau_{b,y}^n -mfu^{n} + r_{y}^n \end{align} :label: u_sit where :math:`r = (r_x,r_y)` contains all terms that do not depend on the velocities :math:`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 :eq:`u_sit` are :math:`u^n` and :math:`v^n`. Once discretized in space, equation :eq:`u_sit` leads to a system of :math:`N` nonlinear equations with :math:`N` unknowns that can be concisely written as .. math:: \mathbf{A}(\mathbf{u})\mathbf{u} = \mathbf{b}(\mathbf{u}), :label: nonlin_sys where :math:`\mathbf{A}` is an :math:`N\times N` matrix and :math:`\mathbf{u}` and :math:`\mathbf{b}` are vectors of size :math:`N`. Note that we have dropped the time level index :math:`n`. The vector :math:`\mathbf{u}` is formed by stacking first the :math:`u` components, followed by the :math:`v` components of the discretized ice velocity. The vector :math:`\mathbf{b}` is a function of the velocity vector :math:`\mathbf{u}` because of the water and seabed stress terms as well as parts of the rheology term that depend non-linearly on :math:`\mathbf{u}`. The nonlinear system :eq:`nonlin_sys` is solved using a Picard iteration method. Starting from a previous iterate :math:`\mathbf{u}_{k-1}`, the nonlinear system is linearized by substituting :math:`\mathbf{u}_{k-1}` in the expression of the matrix :math:`\mathbf{A}` and the vector :math:`\mathbf{b}`: .. math:: \mathbf{A}(\mathbf{u}_{k-1})\mathbf{u}_{k} = \mathbf{b}(\mathbf{u}_{k-1}) :label: picard The resulting linear system is solved using the Flexible Generalized Minimum RESidual (FGMRES, :cite:`Saad93`) method and this process is repeated iteratively. Ice-Ocean stress ~~~~~~~~~~~~~~~~ At the end of each (thermodynamic) time step, the ice–ocean stress must be constructed from :math:`{\tt taux(y)}` and the terms containing :math:`{\tt 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 weight 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: .. 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}, where .. math:: \dot{\epsilon}_{ij} = {1\over 2}\left({{\partial u_i}\over{\partial x_j}} + {{\partial u_j}\over{\partial x_i}}\right) 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 multiplied by :math:`-1` and is therefore simply equal to :math:`-\sigma_1/2`. 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``). .. _stress-vp: Viscous-Plastic ~~~~~~~~~~~~~~~ The VP constitutive law is given by .. math:: \sigma_{ij} = 2 \eta \dot{\epsilon}_{ij} + (\zeta - \eta) D_D - P_R(1 - k_t)\frac{\delta_{ij}}{2} :label: vp-const where :math:`\eta` and :math:`\zeta` are the bulk and shear viscosities. An elliptical yield curve is used, with the viscosities given by .. math:: \zeta = {P(1+k_t)\over 2\Delta}, .. math:: \eta = {P(1+k_t)\over {2\Delta e^2}}, where .. 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 described in the `Icepack Documentation `_. The parameter :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``. .. _stress-evp: Elastic-Viscous-Plastic ~~~~~~~~~~~~~~~~~~~~~~~ In the EVP model the internal stress tensor is determined from a regularized version of the VP constitutive law :eq:`vp-const`. 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 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`. .. _stress-eap: 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}