/[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.9 by mlosch, Wed Jun 4 13:33:45 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 simulations, 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}, EVP
16      \citep{hunke97};
17    \item ice-ocean stress can be formulated as in \citet{hibler87};
18    \item ice variables are advected by sophisticated advection schemes;
19    \item growth and melt parameterization have been refined and extended
20      in order to allow for automatic differentiation of the code.
21    \end{itemize}
22    The model equations and their numerical realization are summarized
23    below.
24    
25  \subsection{Dynamics}  \subsection{Dynamics}
26  \label{sec:dynamics}  \label{sec:dynamics}
27    
# Line 43  wind/current vectors. Line 64  wind/current vectors.
64    
65  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
66  be related to the ice strain rate and strength by a nonlinear  be related to the ice strain rate and strength by a nonlinear
67  viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}:  viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
68  \begin{equation}  \begin{equation}
69    \label{eq:vpequation}    \label{eq:vpequation}
70    \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 98  minor axis $e$ equal to $2$; they are gi
98      \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)      \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
99      (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 +      (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 +
100      2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})      2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
101    \right]^{-\frac{1}{2}}    \right]^{\frac{1}{2}}
102  \end{align*}  \end{align*}
103  The bulk viscosities are bounded above by imposing both a minimum  The bulk viscosities are bounded above by imposing both a minimum
104  $\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 147  timestep. According to the recommendatio
147  EVP-model is stepped forward in time 120 times within the physical  EVP-model is stepped forward in time 120 times within the physical
148  ocean model time step (although this parameter is under debate), to  ocean model time step (although this parameter is under debate), to
149  allow for elastic waves to disappear.  Because the scheme does not  allow for elastic waves to disappear.  Because the scheme does not
150  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
151    timestep and simple to implement on parallel computers
152  \citep{hunke97}. For completeness, we repeat the equations for the  \citep{hunke97}. For completeness, we repeat the equations for the
153  components of the stress tensor $\sigma_{1} =  components of the stress tensor $\sigma_{1} =
154  \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 156  $\sigma_{12}$. Introducing the divergenc
156  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
157  and shearing strain rates, $D_T =  and shearing strain rates, $D_T =
158  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
159  2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,  2\dot{\epsilon}_{12}$, respectively, and using the above
160  the equations can be written as:  abbreviations, the equations\refeq{evpequation} can be written as:
161  \begin{align}  \begin{align}
162    \label{eq:evpstresstensor1}    \label{eq:evpstresstensor1}
163    \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +    \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
# Line 205  advect tracers, is really an average ove Line 227  advect tracers, is really an average ove
227  velocity and the ice velocity leading to an inconsistency as the ice  velocity and the ice velocity leading to an inconsistency as the ice
228  temperature and salinity are different from the oceanic variables.  temperature and salinity are different from the oceanic variables.
229    
230  Sea ice distributions are characterized by sharp gradients and edges.  %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
231  For this reason, we employ positive, multidimensional 2nd and 3rd-order  %\begin{itemize}
232  advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the  %\item transition from existing B-Grid to C-Grid
233  thermodynamic variables discussed in the next section.  %\item boundary conditions: no-slip, free-slip, half-slip
234    %\item fancy (multi dimensional) advection schemes
235  \subparagraph{boundary conditions: no-slip, free-slip, half-slip}  %\item VP vs.\ EVP \citep{hunke97}
236    %\item ocean stress formulation \citep{hibler87}
237  \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}  
238    
239  \subsection{Thermodynamics}  \subsection{Thermodynamics}
240  \label{sec:thermodynamics}  \label{sec:thermodynamics}
# Line 227  In the original formulation the sea ice Line 243  In the original formulation the sea ice
243  uses simple thermodynamics following the appendix of  uses simple thermodynamics following the appendix of
244  \citet{semtner76}. This formulation does not allow storage of heat  \citet{semtner76}. This formulation does not allow storage of heat
245  (heat capacity of ice is zero, and this type of model is often refered  (heat capacity of ice is zero, and this type of model is often refered
246  to as a ``zero-layer'' model). Upward heat flux is parameterized  to as a ``zero-layer'' model). Upward conductive heat flux is parameterized
247  assuming a linear temperature profile and together with a constant ice  assuming a linear temperature profile and together with a constant ice
248  conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is  conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is
249  the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the  the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the
250  difference between water and ice surface temperatures. The surface  difference between water and ice surface temperatures. The surface
251  heat budget is computed in a similar way to that of  heat flux is computed in a similar way to that of \citet{parkinson79}
252  \citet{parkinson79} and \citet{manabe79}.  and \citet{manabe79}.
253    
254    The conductive heat flux depends strongly on the ice thickness $h$.
255    However, the ice thickness in the model represents a mean over a
256    potentially very heterogeneous thickness distribution.  In order to
257    parameterize this sub-grid scale distribution for heat flux
258    computations, the mean ice thickness $h$ is split into seven thickness
259    categories $H_{n}$ that are equally distributed between $2h$ and
260    a minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
261    \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
262    thickness category area averaged to give the total heat flux. \ml{[I
263      don't have citation for that, anyone?]}
264    
265    The atmospheric heat flux is balanced by an oceanic heat flux from
266    below.  The oceanic flux is proportional to
267    $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are
268    the density and heat capacity of sea water and $T_{fr}$ is the local
269    freezing point temperature that is a function of salinity. Contrary to
270    \citet{menemenlis05}, this flux is not assumed to instantaneously melt
271    or create ice, but a time scale of three days is used to relax $T_{w}$
272    to the freezing point.
273    
274    The parameterization of lateral and vertical growth of sea ice follows
275    that of \citet{hibler79, hibler80}.
276    
277    On top of the ice there is a layer of snow that modifies the heat flux
278    and the albedo \citep{zhang98}. If enough snow accumulates so that its
279    weight submerges the ice and the snow is flooded, a simple mass
280    conserving parameterization of snowice formation (a flood-freeze
281    algorithm following Archimedes' principle) turns snow into ice until
282    the ice surface is back at $z=0$ \citep{leppaeranta83}.
283    
284    Effective ice thickness (ice volume per unit area,
285    $c\cdot{h}$), concentration $c$ and effective snow thickness
286    ($c\cdot{h}_{s}$) are advected by ice velocities:
287    \begin{equation}
288      \label{eq:advection}
289      \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
290      \Gamma_{X} + D_{X}
291    \end{equation}
292    where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
293    diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
294    %
295    From the various advection scheme that are available in the MITgcm
296    \citep{mitgcm02}, we choose flux-limited schemes
297    \citep[multidimensional 2nd and 3rd-order advection scheme with flux
298    limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges
299    that are typical of sea ice distributions and to rule out unphysical
300    over- and undershoots (negative thickness or concentration). These
301    scheme conserve volume and horizontal area and are unconditionally
302    stable, so that we can set $D_{X}=0$.  \ml{[do we need to proove that?
303      can we proove that? citation?]}
304    
305  There is considerable doubt about the reliability of such a simple  There is considerable doubt about the reliability of such a simple
306  thermodynamic model---\citet{semtner84} found significant errors in  thermodynamic model---\citet{semtner84} found significant errors in
# Line 251  the two ice layers and the thickness of Line 318  the two ice layers and the thickness of
318    flux limiter \citep{roe85}.}    flux limiter \citep{roe85}.}
319    
320    
321  \subsection{C-grid}  %\subsection{C-grid}
322  \begin{itemize}  %\begin{itemize}
323  \item no-slip vs. free-slip for both lsr and evp;  %\item no-slip vs. free-slip for both lsr and evp;
324    "diagnostics" to look at and use for comparison  %  "diagnostics" to look at and use for comparison
325    \begin{itemize}  %  \begin{itemize}
326    \item ice transport through Fram Strait/Denmark Strait/Davis  %  \item ice transport through Fram Strait/Denmark Strait/Davis
327      Strait/Bering strait (these are general)  %    Strait/Bering strait (these are general)
328    \item ice transport through narrow passages, e.g.\ Nares-Strait  %  \item ice transport through narrow passages, e.g.\ Nares-Strait
329    \end{itemize}  %  \end{itemize}
330  \item compare different advection schemes (if lsr turns out to be more  %\item compare different advection schemes (if lsr turns out to be more
331    effective, then with lsr otherwise I prefer evp), eg.  %  effective, then with lsr otherwise I prefer evp), eg.
332    \begin{itemize}  %  \begin{itemize}
333    \item default 2nd-order with diff1=0.002  %  \item default 2nd-order with diff1=0.002
334    \item 1st-order upwind with diff1=0.  %  \item 1st-order upwind with diff1=0.
335    \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)  %  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
336    \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)  %  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
337    \end{itemize}  %  \end{itemize}
338    That should be enough. Here, total ice mass and location of ice edge  %  That should be enough. Here, total ice mass and location of ice edge
339    is interesting. However, this comparison can be done in an idealized  %  is interesting. However, this comparison can be done in an idealized
340    domain, may not require full Arctic Domain?  %  domain, may not require full Arctic Domain?
341  \item  %\item
342  Do a little study on the parameters of LSR and EVP  %Do a little study on the parameters of LSR and EVP
343  \begin{enumerate}  %\begin{enumerate}
344  \item convergence of LSR, how many iterations do you need to get a  %\item convergence of LSR, how many iterations do you need to get a
345    true elliptic yield curve  %  true elliptic yield curve
346  \item vary deltaTevp and the relaxation parameter for EVP and see when  %\item vary deltaTevp and the relaxation parameter for EVP and see when
347    the EVP solution breaks down (relative to the forcing time scale).  %  the EVP solution breaks down (relative to the forcing time scale).
348    For this, it is essential that the evp solver gives use "stripeless"  %  For this, it is essential that the evp solver gives use "stripeless"
349    solutions, that is your dtevp = 1sec solutions/or 10sec solutions  %  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
350    with SEAICE\_evpDampC = 615.  %  with SEAICE\_evpDampC = 615.
351  \end{enumerate}  %\end{enumerate}
352  \end{itemize}  
353    %\end{itemize}
354    
355    %%% Local Variables:
356    %%% mode: latex
357    %%% TeX-master: "ceaice"
358    %%% End:
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22