/[MITgcm]/MITgcm_contrib/articles/ceaice/ceaice_model.tex
ViewVC logotype

Diff of /MITgcm_contrib/articles/ceaice/ceaice_model.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.8 by mlosch, Tue Apr 29 14:03:31 2008 UTC revision 1.11 by dimitri, Thu Aug 14 16:12:41 2008 UTC
# Line 5  The MITgcm sea ice model (MITsim) is bas Line 5  The MITgcm sea ice model (MITsim) is bas
5  viscous-plastic (VP) dynamic-thermodynamic sea ice model  viscous-plastic (VP) dynamic-thermodynamic sea ice model
6  \citep{zhang97} first introduced by \citet{hibler79, hibler80}. In  \citep{zhang97} first introduced by \citet{hibler79, hibler80}. In
7  order to adapt this model to the requirements of coupled  order to adapt this model to the requirements of coupled
8  ice-ocean simulations, many important aspects of the original code have  ice-ocean state estimation, many important aspects of the original code have
9  been modified and improved:  been modified and improved:
10  \begin{itemize}  \begin{itemize}
11  \item the code has been rewritten for an Arakawa C-grid, both B- and  \item the code has been rewritten for an Arakawa C-grid, both B- and
12    C-grid variants are available; the C-grid code allows for no-slip    C-grid variants are available; the C-grid code allows for no-slip
13    and free-slip lateral boundary conditions;    and free-slip lateral boundary conditions;
14  \item two different solution methods for solving the nonlinear  \item two different solution methods for solving the nonlinear
15    momentum equations have been adopted: LSOR \citep{zhang97}, EVP    momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
16    \citep{hunke97};    \citep{hunke97};
17  \item ice-ocean stress can be formulated as in \citet{hibler87};  \item ice-ocean stress can be formulated as in \citet{hibler87} or as in \citet{cam08};
18  \item ice variables are advected by sophisticated advection schemes;  \item ice variables \ml{can be} advected by sophisticated, \ml{conservative}
19  \item growth and melt parameterizaion have been refined and extended    advection schemes \ml{with flux limiting};
20    in order to allow for automatic differentiation of the code.  \item growth and melt parameterizations have been refined and extended
21      in order to allow for more stable automatic differentiation of the code.
22  \end{itemize}  \end{itemize}
 The model equations and their numerical realization are summarized  
 below.  
23    
24  \subsection{Dynamics}  The sea ice model is tightly coupled to the ocean compontent of the MITgcm
25  \label{sec:dynamics}  \citep{mar97a}.  Heat, fresh water fluxes and surface stresses are computed
26    from the atmospheric state and modified by the ice model at every time step.
27    The model equations and details of their numerical realization are summarized
28    in the appendix. Further documentation and model code can be found at
29    \url{http://mitgcm.org}.
30    
31  The momentum equation of the sea-ice model is  %\subsection{C-grid}
 \begin{equation}  
   \label{eq:momseaice}  
   m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +  
   \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},  
 \end{equation}  
 where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;  
 $\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;  
 $\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$  
 directions, respectively;  
 $f$ is the Coriolis parameter;  
 $\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,  
 respectively;  
 $g$ is the gravity accelation;  
 $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;  
 $\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential  
 in response to ocean dynamics ($g\eta$) and to atmospheric pressure  
 loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);  
 and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress  
 tensor $\sigma_{ij}$.  
 When using the rescaled vertical coordinate system, z$^\ast$, of  
 \citet{cam08}, $\phi(0)$ also includes a term due to snow and ice  
 loading, $mg/\rho_{0}$.  
 Advection of sea-ice momentum is neglected. The wind and ice-ocean stress  
 terms are given by  
 \begin{align*}  
   \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|  
                    R_{air}  (\vek{U}_{air}  -\vek{u}), \\  
   \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|  
                    R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\  
 \end{align*}  
 where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere  
 and surface currents of the ocean, respectively; $C_{air/ocean}$ are  
 air and ocean drag coefficients; $\rho_{air/ocean}$ are reference  
 densities; and $R_{air/ocean}$ are rotation matrices that act on the  
 wind/current vectors.  
   
 For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can  
 be related to the ice strain rate and strength by a nonlinear  
 viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:  
 \begin{equation}  
   \label{eq:vpequation}  
   \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}  
   + \left[\zeta(\dot{\epsilon}_{ij},P) -  
     \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}    
   - \frac{P}{2}\delta_{ij}.  
 \end{equation}  
 The ice strain rate is given by  
 \begin{equation*}  
   \dot{\epsilon}_{ij} = \frac{1}{2}\left(  
     \frac{\partial{u_{i}}}{\partial{x_{j}}} +  
     \frac{\partial{u_{j}}}{\partial{x_{i}}}\right).  
 \end{equation*}  
 The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on  
 both thickness $h$ and compactness (concentration) $c$:  
 \begin{equation}  
   P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},  
 \label{eq:icestrength}  
 \end{equation}  
 with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear  
 viscosities $\eta$ and $\zeta$ are functions of ice strain rate  
 invariants and ice strength such that the principal components of the  
 stress lie on an elliptical yield curve with the ratio of major to  
 minor axis $e$ equal to $2$; they are given by:  
 \begin{align*}  
   \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},  
    \zeta_{\max}\right) \\  
   \eta =& \frac{\zeta}{e^2} \\  
   \intertext{with the abbreviation}  
   \Delta = & \left[  
     \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)  
     (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 +  
     2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})  
   \right]^{\frac{1}{2}}  
 \end{align*}  
 The bulk viscosities are bounded above by imposing both a minimum  
 $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a  
 maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where  
 $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress  
 tensor computation the replacement pressure $P = 2\,\Delta\zeta$  
 \citep{hibler95} is used so that the stress state always lies on the  
 elliptic yield curve by definition.  
   
 In the so-called truncated ellipse method the shear viscosity $\eta$  
 is capped to suppress any tensile stress \citep{hibler97, geiger98}:  
 \begin{equation}  
   \label{eq:etatem}  
   \eta = \min\left(\frac{\zeta}{e^2},  
   \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}  
   {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2  
       +4\dot{\epsilon}_{12}^2}}\right).  
 \end{equation}  
   
 In the current implementation, the VP-model is integrated with the  
 semi-implicit line successive over relaxation (LSOR)-solver of  
 \citet{zhang98}, which allows for long time steps that, in our case,  
 are limited by the explicit treatment of the Coriolis term. The  
 explicit treatment of the Coriolis term does not represent a severe  
 limitation because it restricts the time step to approximately the  
 same length as in the ocean model where the Coriolis term is also  
 treated explicitly.  
   
 \citet{hunke97}'s introduced an elastic contribution to the strain  
 rate in order to regularize Eq.\refeq{vpequation} in such a way that  
 the resulting elastic-viscous-plastic (EVP) and VP models are  
 identical at steady state,  
 \begin{equation}  
   \label{eq:evpequation}  
   \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +  
   \frac{1}{2\eta}\sigma_{ij}  
   + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}    
   + \frac{P}{4\zeta}\delta_{ij}  
   = \dot{\epsilon}_{ij}.  
 \end{equation}  
 %In the EVP model, equations for the components of the stress tensor  
 %$\sigma_{ij}$ are solved explicitly. Both model formulations will be  
 %used and compared the present sea-ice model study.  
 The EVP-model uses an explicit time stepping scheme with a short  
 timestep. According to the recommendation of \citet{hunke97}, the  
 EVP-model is stepped forward in time 120 times within the physical  
 ocean model time step (although this parameter is under debate), to  
 allow for elastic waves to disappear.  Because the scheme does not  
 require a matrix inversion it is fast in spite of the small timestep  
 \citep{hunke97}. For completeness, we repeat the equations for the  
 components of the stress tensor $\sigma_{1} =  
 \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and  
 $\sigma_{12}$. Introducing the divergence $D_D =  
 \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension  
 and shearing strain rates, $D_T =  
 \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =  
 2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,  
 the equations can be written as:  
 \begin{align}  
   \label{eq:evpstresstensor1}  
   \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +  
   \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\  
   \label{eq:evpstresstensor2}  
   \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}  
   &= \frac{P}{2T\Delta} D_T \\  
   \label{eq:evpstresstensor12}  
   \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}  
   &= \frac{P}{4T\Delta} D_S  
 \end{align}  
 Here, the elastic parameter $E$ is redefined in terms of a damping timescale  
 $T$ for elastic waves \[E=\frac{\zeta}{T}.\]  
 $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and  
 the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend  
 $E_{0} = \frac{1}{3}$.  
   
 For details of the spatial discretization, the reader is referred to  
 \citet{zhang98, zhang03}. Our discretization differs only (but  
 importantly) in the underlying grid, namely the Arakawa C-grid, but is  
 otherwise straightforward. The EVP model, in particular, is discretized  
 naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the  
 center points and $\sigma_{12}$ on the corner (or vorticity) points of  
 the grid. With this choice all derivatives are discretized as central  
 differences and averaging is only involved in computing $\Delta$ and  
 $P$ at vorticity points.  
   
 For a general curvilinear grid, one needs in principle to take metric  
 terms into account that arise in the transformation of a curvilinear  
 grid on the sphere. For now, however, we can neglect these metric  
 terms because they are very small on the \ml{[modify following  
   section3:] cubed sphere grids used in this paper; in particular,  
 only near the edges of the cubed sphere grid, we expect them to be  
 non-zero, but these edges are at approximately 35\degS\ or 35\degN\  
 which are never covered by sea-ice in our simulations.  Everywhere  
 else the coordinate system is locally nearly cartesian.}  However, for  
 last-glacial-maximum or snowball-earth-like simulations the question  
 of metric terms needs to be reconsidered.  Either, one includes these  
 terms as in \citet{zhang03}, or one finds a vector-invariant  
 formulation for the sea-ice internal stress term that does not require  
 any metric terms, as it is done in the ocean dynamics of the MITgcm  
 \citep{adcroft04:_cubed_sphere}.  
   
 Lateral boundary conditions are naturally ``no-slip'' for B-grids, as  
 the tangential velocities points lie on the boundary. For C-grids, the  
 lateral boundary condition for tangential velocities is realized via  
 ``ghost points'', allowing alternatively no-slip or free-slip  
 conditions. In ocean models free-slip boundary conditions in  
 conjunction with piecewise-constant (``castellated'') coastlines have  
 been shown to reduce in effect to no-slip boundary conditions  
 \citep{adcroft98:_slippery_coast}; for sea-ice models the effects of  
 lateral boundary conditions have not yet been studied.  
   
 Moving sea ice exerts a stress on the ocean which is the opposite of  
 the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is  
 applied directly to the surface layer of the ocean model. An  
 alternative ocean stress formulation is given by \citet{hibler87}.  
 Rather than applying $\vtau_{ocean}$ directly, the stress is derived  
 from integrating over the ice thickness to the bottom of the oceanic  
 surface layer. In the resulting equation for the \emph{combined}  
 ocean-ice momentum, the interfacial stress cancels and the total  
 stress appears as the sum of windstress and divergence of internal ice  
 stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also  
 Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that  
 now the velocity in the surface layer of the ocean that is used to  
 advect tracers, is really an average over the ocean surface  
 velocity and the ice velocity leading to an inconsistency as the ice  
 temperature and salinity are different from the oceanic variables.  
   
 %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}  
32  %\begin{itemize}  %\begin{itemize}
33  %\item transition from existing B-Grid to C-Grid  %\item no-slip vs. free-slip for both lsr and evp;
34  %\item boundary conditions: no-slip, free-slip, half-slip  %  "diagnostics" to look at and use for comparison
35  %\item fancy (multi dimensional) advection schemes  %  \begin{itemize}
36  %\item VP vs.\ EVP \citep{hunke97}  %  \item ice transport through Fram Strait/Denmark Strait/Davis
37  %\item ocean stress formulation \citep{hibler87}  %    Strait/Bering strait (these are general)
38  %\end{itemize}  %  \item ice transport through narrow passages, e.g.\ Nares-Strait
39    %  \end{itemize}
40  \subsection{Thermodynamics}  %\item compare different advection schemes (if lsr turns out to be more
41  \label{sec:thermodynamics}  %  effective, then with lsr otherwise I prefer evp), eg.
42    %  \begin{itemize}
43  In the original formulation the sea ice model \citep{menemenlis05}  %  \item default 2nd-order with diff1=0.002
44  uses simple thermodynamics following the appendix of  %  \item 1st-order upwind with diff1=0.
45  \citet{semtner76}. This formulation does not allow storage of heat  %  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
46  (heat capacity of ice is zero, and this type of model is often refered  %  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
47  to as a ``zero-layer'' model). Upward conductive heat flux is parameterized  %  \end{itemize}
48  assuming a linear temperature profile and together with a constant ice  %  That should be enough. Here, total ice mass and location of ice edge
49  conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is  %  is interesting. However, this comparison can be done in an idealized
50  the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the  %  domain, may not require full Arctic Domain?
51  difference between water and ice surface temperatures. The surface  %\item
52  heat flux is computed in a similar way to that of \citet{parkinson79}  %Do a little study on the parameters of LSR and EVP
53  and \citet{manabe79}.  %\begin{enumerate}
54    %\item convergence of LSR, how many iterations do you need to get a
55    %  true elliptic yield curve
56    %\item vary deltaTevp and the relaxation parameter for EVP and see when
57    %  the EVP solution breaks down (relative to the forcing time scale).
58    %  For this, it is essential that the evp solver gives use "stripeless"
59    %  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
60    %  with SEAICE\_evpDampC = 615.
61    %\end{enumerate}
62    
63  The conductive heat flux depends strongly on the ice thickness $h$.  %\end{itemize}
 However, the ice thickness in the model represents a mean over a  
 potentially very heterogeneous thickness distribution.  In order to  
 parameterize this sub-grid scale distribution for heat flux  
 computations, the mean ice thickness $h$ is split into seven thickness  
 categories $H_{n}$ that are equally distributed between $2h$ and  
 minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=  
 \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each  
 thickness category area averaged to give the total heat flux. \ml{[I  
   don't have citation for that, anyone?]}  
   
 The atmospheric heat flux is balanced by an oceanic heat flux from  
 below.  The oceanic flux is proportional to  
 $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are  
 the density and heat capacity of sea water and $T_{fr}$ is the local  
 freezing point temperature that is a function of salinity. Contrary to  
 \citet{menemenlis05}, this flux is not assumed to instantaneously melt  
 or create ice, but a time scale of three days is used to relax $T_{w}$  
 to the freezing point.  
   
 The parameterization of lateral and vertical growth of sea ice follows  
 that of \citet{hibler79, hibler80}.  
   
 On top of the ice there is a layer of snow that modifies the heat flux  
 and the albedo \citep{zhang98}. If enough snow accumulates so that its  
 weight submerges the ice and the snow is flooded, a simple mass  
 conserving parameterization of snowice formation (a flood-freeze  
 algorithm following Archimedes' principle) turns snow into ice until  
 the ice surface is back at $z=0$ \citep{leppaeranta83}.  
   
 Effective ice thickness (ice volume per unit area,  
 $c\cdot{h}$), concentration $c$ and effective snow thickness  
 ($c\cdot{h}_{s}$) are advected by ice velocities:  
 \begin{equation}  
   \label{eq:advection}  
   \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +  
   \Gamma_{X} + D_{X}  
 \end{equation}  
 where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the  
 diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.  
 %  
 From the various advection scheme that are available in the MITgcm  
 \citep{mitgcm02}, we choose flux-limited schemes  
 \citep[multidimensional 2nd and 3rd-order advection scheme with flux  
 limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges  
 that are typical of sea ice distributions and to rule out unphysical  
 over- and undershoots (negative thickness or concentration). These  
 scheme conserve volume and horizontal area and are unconditionally  
 stable, so that we can set $D_{X}=0$.  \ml{[do we need to proove that?  
   can we proove that? citation?]}  
   
 There is considerable doubt about the reliability of such a simple  
 thermodynamic model---\citet{semtner84} found significant errors in  
 phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in  
 such models---, so that today many sea ice models employ more complex  
 thermodynamics. A popular thermodynamics model of \citet{winton00} is  
 based on the 3-layer model of \citet{semtner76} and treats brine  
 content by means of enthalphy conservation. This model requires in  
 addition to ice-thickness and compactness (fractional area) additional  
 state variables to be advected by ice velocities, namely enthalphy of  
 the two ice layers and the thickness of the overlying snow layer.  
 \ml{[Jean-Michel, your turf: ]Care must be taken in advecting these  
   quantities in order to ensure conservation of enthalphy. Currently  
   this can only be accomplished with a 2nd-order advection scheme with  
   flux limiter \citep{roe85}.}  
   
   
 \subsection{C-grid}  
 \begin{itemize}  
 \item no-slip vs. free-slip for both lsr and evp;  
   "diagnostics" to look at and use for comparison  
   \begin{itemize}  
   \item ice transport through Fram Strait/Denmark Strait/Davis  
     Strait/Bering strait (these are general)  
   \item ice transport through narrow passages, e.g.\ Nares-Strait  
   \end{itemize}  
 \item compare different advection schemes (if lsr turns out to be more  
   effective, then with lsr otherwise I prefer evp), eg.  
   \begin{itemize}  
   \item default 2nd-order with diff1=0.002  
   \item 1st-order upwind with diff1=0.  
   \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)  
   \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)  
   \end{itemize}  
   That should be enough. Here, total ice mass and location of ice edge  
   is interesting. However, this comparison can be done in an idealized  
   domain, may not require full Arctic Domain?  
 \item  
 Do a little study on the parameters of LSR and EVP  
 \begin{enumerate}  
 \item convergence of LSR, how many iterations do you need to get a  
   true elliptic yield curve  
 \item vary deltaTevp and the relaxation parameter for EVP and see when  
   the EVP solution breaks down (relative to the forcing time scale).  
   For this, it is essential that the evp solver gives use "stripeless"  
   solutions, that is your dtevp = 1sec solutions/or 10sec solutions  
   with SEAICE\_evpDampC = 615.  
 \end{enumerate}  
   
 \end{itemize}  
64    
65  %%% Local Variables:  %%% Local Variables:
66  %%% mode: latex  %%% mode: latex

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.11

  ViewVC Help
Powered by ViewVC 1.1.22