/[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.1 by dimitri, Tue Feb 26 19:27:26 2008 UTC revision 1.12 by dimitri, Fri Aug 15 14:16:35 2008 UTC
# Line 1  Line 1 
1  \section{Model}  \section{Model Formulation}
2  \label{sec:model}  \label{sec:model}
3    
4    The MITgcm sea ice model (MITsim) is based on a variant of the
5    viscous-plastic (VP) dynamic-thermodynamic sea ice model
6    \citep{zhang97} first introduced by \citet{hibler79, hibler80}. In
7    order to adapt this model to the requirements of coupled
8    ice-ocean state estimation, many important aspects of the original code have
9    been modified and improved:
10    \begin{itemize}
11    \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
13      and free-slip lateral boundary conditions;
14    \item two different solution methods for solving the nonlinear
15      momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
16      \citep{hunke97};
17    \item ice-ocean stress can be formulated as in \citet{hibler87} or as in \citet{cam08};
18    \item ice variables \ml{can be} advected by sophisticated, \ml{conservative}
19      advection schemes \ml{with flux limiting};
20    \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}
23    
24    The sea ice model is tightly coupled to the ocean compontent of the MITgcm
25    \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    
28    
29  \subsection{Dynamics}  \subsection{Dynamics}
30  \label{sec:dynamics}  \label{app:dynamics}
31    
32  The momentum equation of the sea-ice model is  The momentum equation of the sea-ice model is
33  \begin{equation}  \begin{equation}
# Line 19  $\vtau_{air}$ and $\vtau_{ocean}$ are th Line 44  $\vtau_{air}$ and $\vtau_{ocean}$ are th
44  respectively;  respectively;
45  $g$ is the gravity accelation;  $g$ is the gravity accelation;
46  $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;  $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
47  $\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential  $\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}$ is the sea surface
48  in response to ocean dynamics ($g\eta$) and to atmospheric pressure  height potential in response to ocean dynamics ($g\eta$), to
49  loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);  atmospheric pressure loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a
50  and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress  reference density) and a term due to snow and ice loading \citep{campin08};
51  tensor $\sigma_{ij}$.  and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice
52  When using the rescaled vertical coordinate system, z$^\ast$, of  stress tensor $\sigma_{ij}$. %
 \citet{cam08}, $\phi(0)$ also includes a term due to snow and ice  
 loading, $mg/\rho_{0}$.  
53  Advection of sea-ice momentum is neglected. The wind and ice-ocean stress  Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
54  terms are given by  terms are given by
55  \begin{align*}  \begin{align*}
# Line 43  wind/current vectors. Line 66  wind/current vectors.
66    
67  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
68  be related to the ice strain rate and strength by a nonlinear  be related to the ice strain rate and strength by a nonlinear
69  viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}:  viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
70  \begin{equation}  \begin{equation}
71    \label{eq:vpequation}    \label{eq:vpequation}
72    \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}    \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
# Line 77  minor axis $e$ equal to $2$; they are gi Line 100  minor axis $e$ equal to $2$; they are gi
100      \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)      \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
101      (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 +      (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 +
102      2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})      2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
103    \right]^{-\frac{1}{2}}    \right]^{\frac{1}{2}}.
104  \end{align*}  \end{align*}
105  The bulk viscosities are bounded above by imposing both a minimum  The bulk viscosities are bounded above by imposing both a minimum
106  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
# Line 126  timestep. According to the recommendatio Line 149  timestep. According to the recommendatio
149  EVP-model is stepped forward in time 120 times within the physical  EVP-model is stepped forward in time 120 times within the physical
150  ocean model time step (although this parameter is under debate), to  ocean model time step (although this parameter is under debate), to
151  allow for elastic waves to disappear.  Because the scheme does not  allow for elastic waves to disappear.  Because the scheme does not
152  require a matrix inversion it is fast in spite of the small timestep  require a matrix inversion it is fast in spite of the small internal
153    timestep and simple to implement on parallel computers
154  \citep{hunke97}. For completeness, we repeat the equations for the  \citep{hunke97}. For completeness, we repeat the equations for the
155  components of the stress tensor $\sigma_{1} =  components of the stress tensor $\sigma_{1} =
156  \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and  \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
# Line 134  $\sigma_{12}$. Introducing the divergenc Line 158  $\sigma_{12}$. Introducing the divergenc
158  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
159  and shearing strain rates, $D_T =  and shearing strain rates, $D_T =
160  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
161  2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,  2\dot{\epsilon}_{12}$, respectively, and using the above
162  the equations can be written as:  abbreviations, the equations\refeq{evpequation} can be written as:
163  \begin{align}  \begin{align}
164    \label{eq:evpstresstensor1}    \label{eq:evpstresstensor1}
165    \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +    \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
# Line 153  $T=E_{0}\Delta{t}$ with the tunable para Line 177  $T=E_{0}\Delta{t}$ with the tunable para
177  the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend  the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend
178  $E_{0} = \frac{1}{3}$.  $E_{0} = \frac{1}{3}$.
179    
180  For details of the spatial discretization, the reader is referred to  Our discretization differs from \citet{zhang98, zhang03} in the
181  \citet{zhang98, zhang03}. Our discretization differs only (but  underlying grid, namely the Arakawa C-grid, but is otherwise
182  importantly) in the underlying grid, namely the Arakawa C-grid, but is  straightforward. The EVP model, in particular, is discretized
 otherwise straightforward. The EVP model, in particular, is discretized  
183  naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the  naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
184  center points and $\sigma_{12}$ on the corner (or vorticity) points of  center points and $\sigma_{12}$ on the corner (or vorticity) points of
185  the grid. With this choice all derivatives are discretized as central  the grid. With this choice all derivatives are discretized as central
# Line 205  advect tracers, is really an average ove Line 228  advect tracers, is really an average ove
228  velocity and the ice velocity leading to an inconsistency as the ice  velocity and the ice velocity leading to an inconsistency as the ice
229  temperature and salinity are different from the oceanic variables.  temperature and salinity are different from the oceanic variables.
230    
231  Sea ice distributions are characterized by sharp gradients and edges.  %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
232  For this reason, we employ positive, multidimensional 2nd and 3rd-order  %\begin{itemize}
233  advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the  %\item transition from existing B-Grid to C-Grid
234  thermodynamic variables discussed in the next section.  %\item boundary conditions: no-slip, free-slip, half-slip
235    %\item fancy (multi dimensional) advection schemes
236  \subparagraph{boundary conditions: no-slip, free-slip, half-slip}  %\item VP vs.\ EVP \citep{hunke97}
237    %\item ocean stress formulation \citep{hibler87}
238  \begin{itemize}  %\end{itemize}
 \item transition from existing B-Grid to C-Grid  
 \item boundary conditions: no-slip, free-slip, half-slip  
 \item fancy (multi dimensional) advection schemes  
 \item VP vs.\ EVP \citep{hunke97}  
 \item ocean stress formulation \citep{hibler87}  
 \end{itemize}  
239    
240  \subsection{Thermodynamics}  \subsection{Thermodynamics}
241  \label{sec:thermodynamics}  \label{app:thermodynamics}
242    
243  In the original formulation the sea ice model \citep{menemenlis05}  In its original formulation the sea ice model \citep{menemenlis05}
244  uses simple thermodynamics following the appendix of  uses simple thermodynamics following the appendix of
245  \citet{semtner76}. This formulation does not allow storage of heat  \citet{semtner76}. This formulation does not allow storage of heat,
246  (heat capacity of ice is zero, and this type of model is often refered  that is, the heat capacity of ice is zero. Upward conductive heat flux
247  to as a ``zero-layer'' model). Upward heat flux is parameterized  is parameterized assuming a linear temperature profile and together
248  assuming a linear temperature profile and together with a constant ice  with a constant ice conductivity. It is expressed as
249  conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is  $(K/h)(T_{w}-T_{0})$, where $K$ is the ice conductivity, $h$ the ice
250  the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the  thickness, and $T_{w}-T_{0}$ the difference between water and ice
251  difference between water and ice surface temperatures. The surface  surface temperatures. This type of model is often refered to as a
252  heat budget is computed in a similar way to that of  ``zero-layer'' model. The surface heat flux is computed in a similar
253  \citet{parkinson79} and \citet{manabe79}.  way to that of \citet{parkinson79} and \citet{manabe79}.
254    
255    The conductive heat flux depends strongly on the ice thickness $h$.
256    However, the ice thickness in the model represents a mean over a
257    potentially very heterogeneous thickness distribution.  In order to
258    parameterize a sub-grid scale distribution for heat flux
259    computations, the mean ice thickness $h$ is split into seven thickness
260    categories $H_{n}$ that are equally distributed between $2h$ and a
261    minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
262    \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
263    thickness category is area-averaged to give the total heat flux
264    \citep{hibler84}.
265    
266    \ml{[This is Ian Fenty's work and we may want to remove this paragraph
267      from the paper]: %
268    The atmospheric heat flux is balanced by an oceanic heat flux from
269    below.  The oceanic flux is proportional to
270    $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are
271    the density and heat capacity of sea water and $T_{fr}$ is the local
272    freezing point temperature that is a function of salinity. Contrary to
273    \citet{menemenlis05}, this flux is not assumed to instantaneously melt
274    or create ice, but a time scale of three days is used to relax $T_{w}$
275    to the freezing point.}
276    %
277    The parameterization of lateral and vertical growth of sea ice follows
278    that of \citet{hibler79, hibler80}.
279    
280    On top of the ice there is a layer of snow that modifies the heat flux
281    and the albedo \citep{zhang98}. Snow modifies the effective
282    conductivity according to
283    \[\frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},\]
284    where $K_s$ is the conductivity of snow and $h_s$ the snow thickness.
285    If enough snow accumulates so that its weight submerges the ice and
286    the snow is flooded, a simple mass conserving parameterization of
287    snowice formation (a flood-freeze algorithm following Archimedes'
288    principle) turns snow into ice until the ice surface is back at $z=0$
289    \citep{leppaeranta83}.
290    
291    Effective ice thickness (ice volume per unit area,
292    $c\cdot{h}$), concentration $c$ and effective snow thickness
293    ($c\cdot{h}_{s}$) are advected by ice velocities:
294    \begin{equation}
295      \label{eq:advection}
296      \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
297      \Gamma_{X} + D_{X}
298    \end{equation}
299    where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
300    diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
301    %
302    From the various advection scheme that are available in the MITgcm
303    \citep{mitgcm02}, we choose flux-limited schemes
304    \citep[multidimensional 2nd and 3rd-order advection scheme with flux
305    limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges
306    that are typical of sea ice distributions and to rule out unphysical
307    over- and undershoots (negative thickness or concentration). These
308    scheme conserve volume and horizontal area and are unconditionally
309    stable, so that we can set $D_{X}=0$.  \ml{[do we need to proove that?
310      can we proove that? citation?]}
311    
312  There is considerable doubt about the reliability of such a simple  There is considerable doubt about the reliability of such a simple
313  thermodynamic model---\citet{semtner84} found significant errors in  thermodynamic model---\citet{semtner84} found significant errors in
# Line 250  the two ice layers and the thickness of Line 324  the two ice layers and the thickness of
324    this can only be accomplished with a 2nd-order advection scheme with    this can only be accomplished with a 2nd-order advection scheme with
325    flux limiter \citep{roe85}.}    flux limiter \citep{roe85}.}
326    
327    Further documentation and model code can be found at
328    \url{http://mitgcm.org}.
329    
330  \subsection{C-grid}  
331  \begin{itemize}  %\subsection{C-grid}
332  \item no-slip vs. free-slip for both lsr and evp;  %\begin{itemize}
333    "diagnostics" to look at and use for comparison  %\item no-slip vs. free-slip for both lsr and evp;
334    \begin{itemize}  %  "diagnostics" to look at and use for comparison
335    \item ice transport through Fram Strait/Denmark Strait/Davis  %  \begin{itemize}
336      Strait/Bering strait (these are general)  %  \item ice transport through Fram Strait/Denmark Strait/Davis
337    \item ice transport through narrow passages, e.g.\ Nares-Strait  %    Strait/Bering strait (these are general)
338    \end{itemize}  %  \item ice transport through narrow passages, e.g.\ Nares-Strait
339  \item compare different advection schemes (if lsr turns out to be more  %  \end{itemize}
340    effective, then with lsr otherwise I prefer evp), eg.  %\item compare different advection schemes (if lsr turns out to be more
341    \begin{itemize}  %  effective, then with lsr otherwise I prefer evp), eg.
342    \item default 2nd-order with diff1=0.002  %  \begin{itemize}
343    \item 1st-order upwind with diff1=0.  %  \item default 2nd-order with diff1=0.002
344    \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)  %  \item 1st-order upwind with diff1=0.
345    \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)  %  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
346    \end{itemize}  %  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
347    That should be enough. Here, total ice mass and location of ice edge  %  \end{itemize}
348    is interesting. However, this comparison can be done in an idealized  %  That should be enough. Here, total ice mass and location of ice edge
349    domain, may not require full Arctic Domain?  %  is interesting. However, this comparison can be done in an idealized
350  \item  %  domain, may not require full Arctic Domain?
351  Do a little study on the parameters of LSR and EVP  %\item
352  \begin{enumerate}  %Do a little study on the parameters of LSR and EVP
353  \item convergence of LSR, how many iterations do you need to get a  %\begin{enumerate}
354    true elliptic yield curve  %\item convergence of LSR, how many iterations do you need to get a
355  \item vary deltaTevp and the relaxation parameter for EVP and see when  %  true elliptic yield curve
356    the EVP solution breaks down (relative to the forcing time scale).  %\item vary deltaTevp and the relaxation parameter for EVP and see when
357    For this, it is essential that the evp solver gives use "stripeless"  %  the EVP solution breaks down (relative to the forcing time scale).
358    solutions, that is your dtevp = 1sec solutions/or 10sec solutions  %  For this, it is essential that the evp solver gives use "stripeless"
359    with SEAICE\_evpDampC = 615.  %  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
360  \end{enumerate}  %  with SEAICE\_evpDampC = 615.
361  \end{itemize}  %\end{enumerate}
362    
363    %\end{itemize}
364    
365    %%% Local Variables:
366    %%% mode: latex
367    %%% TeX-master: "ceaice"
368    %%% End:
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22