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

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

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


Revision 1.6 - (hide annotations) (download) (as text)
Fri Feb 29 16:47:45 2008 UTC (17 years, 4 months ago) by mlosch
Branch: MAIN
Changes since 1.5: +1 -1 lines
File MIME type: application/x-tex
fix a reference

1 mlosch 1.4 \section{Model Formulation}
2 dimitri 1.1 \label{sec:model}
3    
4 mlosch 1.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 mlosch 1.6 momentum equations have been adopted: LSOR \citep{zhang97}, EVP
16 mlosch 1.4 \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 parameterizaion 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 dimitri 1.1 \subsection{Dynamics}
26     \label{sec:dynamics}
27    
28     The momentum equation of the sea-ice model is
29     \begin{equation}
30     \label{eq:momseaice}
31     m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
32     \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
33     \end{equation}
34     where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
35     $\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
36     $\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
37     directions, respectively;
38     $f$ is the Coriolis parameter;
39     $\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
40     respectively;
41     $g$ is the gravity accelation;
42     $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
43     $\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential
44     in response to ocean dynamics ($g\eta$) and to atmospheric pressure
45     loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
46     and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
47     tensor $\sigma_{ij}$.
48     When using the rescaled vertical coordinate system, z$^\ast$, of
49     \citet{cam08}, $\phi(0)$ also includes a term due to snow and ice
50     loading, $mg/\rho_{0}$.
51     Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
52     terms are given by
53     \begin{align*}
54     \vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}|
55     R_{air} (\vek{U}_{air} -\vek{u}), \\
56     \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|
57     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\
58     \end{align*}
59     where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
60     and surface currents of the ocean, respectively; $C_{air/ocean}$ are
61     air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
62     densities; and $R_{air/ocean}$ are rotation matrices that act on the
63     wind/current vectors.
64    
65     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
67 mlosch 1.4 viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
68 dimitri 1.1 \begin{equation}
69     \label{eq:vpequation}
70     \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
71     + \left[\zeta(\dot{\epsilon}_{ij},P) -
72     \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}
73     - \frac{P}{2}\delta_{ij}.
74     \end{equation}
75     The ice strain rate is given by
76     \begin{equation*}
77     \dot{\epsilon}_{ij} = \frac{1}{2}\left(
78     \frac{\partial{u_{i}}}{\partial{x_{j}}} +
79     \frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
80     \end{equation*}
81     The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
82     both thickness $h$ and compactness (concentration) $c$:
83     \begin{equation}
84     P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},
85     \label{eq:icestrength}
86     \end{equation}
87     with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear
88     viscosities $\eta$ and $\zeta$ are functions of ice strain rate
89     invariants and ice strength such that the principal components of the
90     stress lie on an elliptical yield curve with the ratio of major to
91     minor axis $e$ equal to $2$; they are given by:
92     \begin{align*}
93     \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},
94     \zeta_{\max}\right) \\
95     \eta =& \frac{\zeta}{e^2} \\
96     \intertext{with the abbreviation}
97     \Delta = & \left[
98     \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
99     (1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 +
100     2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
101     \right]^{-\frac{1}{2}}
102     \end{align*}
103     The bulk viscosities are bounded above by imposing both a minimum
104     $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
105     maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
106     $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
107     tensor computation the replacement pressure $P = 2\,\Delta\zeta$
108     \citep{hibler95} is used so that the stress state always lies on the
109     elliptic yield curve by definition.
110    
111     In the so-called truncated ellipse method the shear viscosity $\eta$
112     is capped to suppress any tensile stress \citep{hibler97, geiger98}:
113     \begin{equation}
114     \label{eq:etatem}
115     \eta = \min\left(\frac{\zeta}{e^2},
116     \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
117     {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
118     +4\dot{\epsilon}_{12}^2}}\right).
119     \end{equation}
120    
121     In the current implementation, the VP-model is integrated with the
122     semi-implicit line successive over relaxation (LSOR)-solver of
123     \citet{zhang98}, which allows for long time steps that, in our case,
124     are limited by the explicit treatment of the Coriolis term. The
125     explicit treatment of the Coriolis term does not represent a severe
126     limitation because it restricts the time step to approximately the
127     same length as in the ocean model where the Coriolis term is also
128     treated explicitly.
129    
130     \citet{hunke97}'s introduced an elastic contribution to the strain
131     rate in order to regularize Eq.\refeq{vpequation} in such a way that
132     the resulting elastic-viscous-plastic (EVP) and VP models are
133     identical at steady state,
134     \begin{equation}
135     \label{eq:evpequation}
136     \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
137     \frac{1}{2\eta}\sigma_{ij}
138     + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}
139     + \frac{P}{4\zeta}\delta_{ij}
140     = \dot{\epsilon}_{ij}.
141     \end{equation}
142     %In the EVP model, equations for the components of the stress tensor
143     %$\sigma_{ij}$ are solved explicitly. Both model formulations will be
144     %used and compared the present sea-ice model study.
145     The EVP-model uses an explicit time stepping scheme with a short
146     timestep. According to the recommendation of \citet{hunke97}, the
147     EVP-model is stepped forward in time 120 times within the physical
148     ocean model time step (although this parameter is under debate), to
149     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
151     \citep{hunke97}. For completeness, we repeat the equations for the
152     components of the stress tensor $\sigma_{1} =
153     \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
154     $\sigma_{12}$. Introducing the divergence $D_D =
155     \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
156     and shearing strain rates, $D_T =
157     \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
158     2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,
159     the equations can be written as:
160     \begin{align}
161     \label{eq:evpstresstensor1}
162     \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
163     \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
164     \label{eq:evpstresstensor2}
165     \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
166     &= \frac{P}{2T\Delta} D_T \\
167     \label{eq:evpstresstensor12}
168     \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
169     &= \frac{P}{4T\Delta} D_S
170     \end{align}
171     Here, the elastic parameter $E$ is redefined in terms of a damping timescale
172     $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
173     $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and
174     the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend
175     $E_{0} = \frac{1}{3}$.
176    
177     For details of the spatial discretization, the reader is referred to
178     \citet{zhang98, zhang03}. Our discretization differs only (but
179     importantly) in the underlying grid, namely the Arakawa C-grid, but is
180     otherwise straightforward. The EVP model, in particular, is discretized
181     naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
182     center points and $\sigma_{12}$ on the corner (or vorticity) points of
183     the grid. With this choice all derivatives are discretized as central
184     differences and averaging is only involved in computing $\Delta$ and
185     $P$ at vorticity points.
186    
187     For a general curvilinear grid, one needs in principle to take metric
188     terms into account that arise in the transformation of a curvilinear
189     grid on the sphere. For now, however, we can neglect these metric
190     terms because they are very small on the \ml{[modify following
191     section3:] cubed sphere grids used in this paper; in particular,
192     only near the edges of the cubed sphere grid, we expect them to be
193     non-zero, but these edges are at approximately 35\degS\ or 35\degN\
194     which are never covered by sea-ice in our simulations. Everywhere
195     else the coordinate system is locally nearly cartesian.} However, for
196     last-glacial-maximum or snowball-earth-like simulations the question
197     of metric terms needs to be reconsidered. Either, one includes these
198     terms as in \citet{zhang03}, or one finds a vector-invariant
199     formulation for the sea-ice internal stress term that does not require
200     any metric terms, as it is done in the ocean dynamics of the MITgcm
201     \citep{adcroft04:_cubed_sphere}.
202    
203     Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
204     the tangential velocities points lie on the boundary. For C-grids, the
205     lateral boundary condition for tangential velocities is realized via
206     ``ghost points'', allowing alternatively no-slip or free-slip
207     conditions. In ocean models free-slip boundary conditions in
208     conjunction with piecewise-constant (``castellated'') coastlines have
209     been shown to reduce in effect to no-slip boundary conditions
210     \citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
211     lateral boundary conditions have not yet been studied.
212    
213     Moving sea ice exerts a stress on the ocean which is the opposite of
214     the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
215     applied directly to the surface layer of the ocean model. An
216     alternative ocean stress formulation is given by \citet{hibler87}.
217     Rather than applying $\vtau_{ocean}$ directly, the stress is derived
218     from integrating over the ice thickness to the bottom of the oceanic
219     surface layer. In the resulting equation for the \emph{combined}
220     ocean-ice momentum, the interfacial stress cancels and the total
221     stress appears as the sum of windstress and divergence of internal ice
222     stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
223     Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
224     now the velocity in the surface layer of the ocean that is used to
225     advect tracers, is really an average over the ocean surface
226     velocity and the ice velocity leading to an inconsistency as the ice
227     temperature and salinity are different from the oceanic variables.
228    
229 mlosch 1.4 %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
230     %\begin{itemize}
231     %\item transition from existing B-Grid to C-Grid
232     %\item boundary conditions: no-slip, free-slip, half-slip
233     %\item fancy (multi dimensional) advection schemes
234     %\item VP vs.\ EVP \citep{hunke97}
235     %\item ocean stress formulation \citep{hibler87}
236     %\end{itemize}
237 dimitri 1.1
238     \subsection{Thermodynamics}
239     \label{sec:thermodynamics}
240    
241     In the original formulation the sea ice model \citep{menemenlis05}
242     uses simple thermodynamics following the appendix of
243     \citet{semtner76}. This formulation does not allow storage of heat
244     (heat capacity of ice is zero, and this type of model is often refered
245 mlosch 1.3 to as a ``zero-layer'' model). Upward conductive heat flux is parameterized
246 dimitri 1.1 assuming a linear temperature profile and together with a constant ice
247     conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is
248     the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the
249     difference between water and ice surface temperatures. The surface
250 mlosch 1.3 heat flux is computed in a similar way to that of \citet{parkinson79}
251     and \citet{manabe79}.
252    
253     The conductive heat flux depends strongly on the ice thickness $h$.
254     However, the ice thickness in the model represents a mean over a
255     potentially very heterogeneous thickness distribution. In order to
256     parameterize this sub-grid scale distribution for heat flux
257     computations, the mean ice thickness $h$ is split into seven thickness
258     categories $H_{n}$ that are equally distributed between $2h$ and
259     minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
260 mlosch 1.4 \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
261     thickness category area averaged to give the total heat flux. \ml{[I
262     don't have citation for that, anyone?]}
263 mlosch 1.3
264     The atmospheric heat flux is balanced by an oceanic heat flux from
265     below. The oceanic flux is proportional to
266     $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are
267     the density and heat capacity of sea water and $T_{fr}$ is the local
268     freezing point temperature that is a function of salinity. Contrary to
269     \citet{menemenlis05}, this flux is not assumed to instantaneously melt
270     or create ice, but a time scale of three days is used to relax $T_{w}$
271     to the freezing point.
272    
273     The parameterization of lateral and vertical growth of sea ice follows
274     that of \citet{hibler79, hibler80}.
275    
276     On top of the ice there is a layer of snow that modifies the heat flux
277     and the albedo \citep{zhang98}. If enough snow accumulates so that its
278     weight submerges the ice and the snow is flooded, a simple mass
279     conserving parameterization of snowice formation (a flood-freeze
280     algorithm following Archimedes' principle) turns snow into ice until
281     the ice surface is back at $z=0$ \citep{leppaeranta83}.
282    
283     Effective ich thickness (ice volume per unit area,
284     $c\cdot{h}$), concentration $c$ and effective snow thickness
285 mlosch 1.4 ($c\cdot{h}_{s}$) are advected by ice velocities:
286 mlosch 1.5 \begin{equation}
287     \label{eq:advection}
288     \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
289     \Gamma_{X} + D_{X}
290     \end{equation}
291 mlosch 1.4 where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
292 mlosch 1.5 diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
293 mlosch 1.4 %
294     From the various advection scheme that are available in the MITgcm
295     \citep{mitgcm02}, we choose flux-limited schemes
296     \citep[multidimensional 2nd and 3rd-order advection scheme with flux
297     limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges
298     that are typical of sea ice distributions and to rule out unphysical
299     over- and undershoots (negative thickness or concentration). These
300     scheme conserve volume and horizontal area and are unconditionally
301     stable, so that we can set $D_{X}=0$. \ml{[do we need to proove that?
302     can we proove that? citation?]}
303 dimitri 1.1
304     There is considerable doubt about the reliability of such a simple
305     thermodynamic model---\citet{semtner84} found significant errors in
306     phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
307     such models---, so that today many sea ice models employ more complex
308     thermodynamics. A popular thermodynamics model of \citet{winton00} is
309     based on the 3-layer model of \citet{semtner76} and treats brine
310     content by means of enthalphy conservation. This model requires in
311     addition to ice-thickness and compactness (fractional area) additional
312     state variables to be advected by ice velocities, namely enthalphy of
313     the two ice layers and the thickness of the overlying snow layer.
314     \ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
315     quantities in order to ensure conservation of enthalphy. Currently
316     this can only be accomplished with a 2nd-order advection scheme with
317     flux limiter \citep{roe85}.}
318    
319    
320     \subsection{C-grid}
321     \begin{itemize}
322     \item no-slip vs. free-slip for both lsr and evp;
323     "diagnostics" to look at and use for comparison
324     \begin{itemize}
325     \item ice transport through Fram Strait/Denmark Strait/Davis
326     Strait/Bering strait (these are general)
327     \item ice transport through narrow passages, e.g.\ Nares-Strait
328     \end{itemize}
329     \item compare different advection schemes (if lsr turns out to be more
330     effective, then with lsr otherwise I prefer evp), eg.
331     \begin{itemize}
332     \item default 2nd-order with diff1=0.002
333     \item 1st-order upwind with diff1=0.
334     \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
335     \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
336     \end{itemize}
337     That should be enough. Here, total ice mass and location of ice edge
338     is interesting. However, this comparison can be done in an idealized
339     domain, may not require full Arctic Domain?
340     \item
341     Do a little study on the parameters of LSR and EVP
342     \begin{enumerate}
343     \item convergence of LSR, how many iterations do you need to get a
344     true elliptic yield curve
345     \item vary deltaTevp and the relaxation parameter for EVP and see when
346     the EVP solution breaks down (relative to the forcing time scale).
347     For this, it is essential that the evp solver gives use "stripeless"
348     solutions, that is your dtevp = 1sec solutions/or 10sec solutions
349     with SEAICE\_evpDampC = 615.
350     \end{enumerate}
351 dimitri 1.2
352 dimitri 1.1 \end{itemize}
353 mlosch 1.3
354     %%% Local Variables:
355     %%% mode: latex
356     %%% TeX-master: "ceaice"
357     %%% End:

  ViewVC Help
Powered by ViewVC 1.1.22