/[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.9 - (hide annotations) (download) (as text)
Wed Jun 4 13:33:45 2008 UTC (17 years, 1 month ago) by mlosch
Branch: MAIN
Changes since 1.8: +38 -37 lines
File MIME type: application/x-tex
minor changes, typos, comment out large parts of sketch

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 mlosch 1.9 \item growth and melt parameterization have been refined and extended
20 mlosch 1.4 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 mlosch 1.7 \right]^{\frac{1}{2}}
102 dimitri 1.1 \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 mlosch 1.9 require a matrix inversion it is fast in spite of the small internal
151     timestep and simple to implement on parallel computers
152 dimitri 1.1 \citep{hunke97}. For completeness, we repeat the equations for the
153     components of the stress tensor $\sigma_{1} =
154     \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
155     $\sigma_{12}$. Introducing the divergence $D_D =
156     \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
157     and shearing strain rates, $D_T =
158     \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
159 mlosch 1.9 2\dot{\epsilon}_{12}$, respectively, and using the above
160     abbreviations, the equations\refeq{evpequation} can be written as:
161 dimitri 1.1 \begin{align}
162     \label{eq:evpstresstensor1}
163     \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
164     \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
165     \label{eq:evpstresstensor2}
166     \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
167     &= \frac{P}{2T\Delta} D_T \\
168     \label{eq:evpstresstensor12}
169     \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
170     &= \frac{P}{4T\Delta} D_S
171     \end{align}
172     Here, the elastic parameter $E$ is redefined in terms of a damping timescale
173     $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
174     $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and
175     the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend
176     $E_{0} = \frac{1}{3}$.
177    
178     For details of the spatial discretization, the reader is referred to
179     \citet{zhang98, zhang03}. Our discretization differs only (but
180     importantly) in the underlying grid, namely the Arakawa C-grid, but is
181     otherwise straightforward. The EVP model, in particular, is discretized
182     naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
183     center points and $\sigma_{12}$ on the corner (or vorticity) points of
184     the grid. With this choice all derivatives are discretized as central
185     differences and averaging is only involved in computing $\Delta$ and
186     $P$ at vorticity points.
187    
188     For a general curvilinear grid, one needs in principle to take metric
189     terms into account that arise in the transformation of a curvilinear
190     grid on the sphere. For now, however, we can neglect these metric
191     terms because they are very small on the \ml{[modify following
192     section3:] cubed sphere grids used in this paper; in particular,
193     only near the edges of the cubed sphere grid, we expect them to be
194     non-zero, but these edges are at approximately 35\degS\ or 35\degN\
195     which are never covered by sea-ice in our simulations. Everywhere
196     else the coordinate system is locally nearly cartesian.} However, for
197     last-glacial-maximum or snowball-earth-like simulations the question
198     of metric terms needs to be reconsidered. Either, one includes these
199     terms as in \citet{zhang03}, or one finds a vector-invariant
200     formulation for the sea-ice internal stress term that does not require
201     any metric terms, as it is done in the ocean dynamics of the MITgcm
202     \citep{adcroft04:_cubed_sphere}.
203    
204     Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
205     the tangential velocities points lie on the boundary. For C-grids, the
206     lateral boundary condition for tangential velocities is realized via
207     ``ghost points'', allowing alternatively no-slip or free-slip
208     conditions. In ocean models free-slip boundary conditions in
209     conjunction with piecewise-constant (``castellated'') coastlines have
210     been shown to reduce in effect to no-slip boundary conditions
211     \citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
212     lateral boundary conditions have not yet been studied.
213    
214     Moving sea ice exerts a stress on the ocean which is the opposite of
215     the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
216     applied directly to the surface layer of the ocean model. An
217     alternative ocean stress formulation is given by \citet{hibler87}.
218     Rather than applying $\vtau_{ocean}$ directly, the stress is derived
219     from integrating over the ice thickness to the bottom of the oceanic
220     surface layer. In the resulting equation for the \emph{combined}
221     ocean-ice momentum, the interfacial stress cancels and the total
222     stress appears as the sum of windstress and divergence of internal ice
223     stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
224     Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
225     now the velocity in the surface layer of the ocean that is used to
226     advect tracers, is really an average over the ocean surface
227     velocity and the ice velocity leading to an inconsistency as the ice
228     temperature and salinity are different from the oceanic variables.
229    
230 mlosch 1.4 %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
231     %\begin{itemize}
232     %\item transition from existing B-Grid to C-Grid
233     %\item boundary conditions: no-slip, free-slip, half-slip
234     %\item fancy (multi dimensional) advection schemes
235     %\item VP vs.\ EVP \citep{hunke97}
236     %\item ocean stress formulation \citep{hibler87}
237     %\end{itemize}
238 dimitri 1.1
239     \subsection{Thermodynamics}
240     \label{sec:thermodynamics}
241    
242     In the original formulation the sea ice model \citep{menemenlis05}
243     uses simple thermodynamics following the appendix of
244     \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
246 mlosch 1.3 to as a ``zero-layer'' model). Upward conductive heat flux is parameterized
247 dimitri 1.1 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
249     the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the
250     difference between water and ice surface temperatures. The surface
251 mlosch 1.3 heat flux is computed in a similar way to that of \citet{parkinson79}
252     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 mlosch 1.9 a minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
261 mlosch 1.4 \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 mlosch 1.3
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 mlosch 1.8 Effective ice thickness (ice volume per unit area,
285 mlosch 1.3 $c\cdot{h}$), concentration $c$ and effective snow thickness
286 mlosch 1.4 ($c\cdot{h}_{s}$) are advected by ice velocities:
287 mlosch 1.5 \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 mlosch 1.4 where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
293 mlosch 1.5 diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
294 mlosch 1.4 %
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 dimitri 1.1
305     There is considerable doubt about the reliability of such a simple
306     thermodynamic model---\citet{semtner84} found significant errors in
307     phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
308     such models---, so that today many sea ice models employ more complex
309     thermodynamics. A popular thermodynamics model of \citet{winton00} is
310     based on the 3-layer model of \citet{semtner76} and treats brine
311     content by means of enthalphy conservation. This model requires in
312     addition to ice-thickness and compactness (fractional area) additional
313     state variables to be advected by ice velocities, namely enthalphy of
314     the two ice layers and the thickness of the overlying snow layer.
315     \ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
316     quantities in order to ensure conservation of enthalphy. Currently
317     this can only be accomplished with a 2nd-order advection scheme with
318     flux limiter \citep{roe85}.}
319    
320    
321 mlosch 1.9 %\subsection{C-grid}
322     %\begin{itemize}
323     %\item no-slip vs. free-slip for both lsr and evp;
324     % "diagnostics" to look at and use for comparison
325     % \begin{itemize}
326     % \item ice transport through Fram Strait/Denmark Strait/Davis
327     % Strait/Bering strait (these are general)
328     % \item ice transport through narrow passages, e.g.\ Nares-Strait
329     % \end{itemize}
330     %\item compare different advection schemes (if lsr turns out to be more
331     % effective, then with lsr otherwise I prefer evp), eg.
332     % \begin{itemize}
333     % \item default 2nd-order with diff1=0.002
334     % \item 1st-order upwind with diff1=0.
335     % \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
336     % \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
337     % \end{itemize}
338     % 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
340     % domain, may not require full Arctic Domain?
341     %\item
342     %Do a little study on the parameters of LSR and EVP
343     %\begin{enumerate}
344     %\item convergence of LSR, how many iterations do you need to get a
345     % true elliptic yield curve
346     %\item vary deltaTevp and the relaxation parameter for EVP and see when
347     % the EVP solution breaks down (relative to the forcing time scale).
348     % For this, it is essential that the evp solver gives use "stripeless"
349     % solutions, that is your dtevp = 1sec solutions/or 10sec solutions
350     % with SEAICE\_evpDampC = 615.
351     %\end{enumerate}
352 dimitri 1.2
353 mlosch 1.9 %\end{itemize}
354 mlosch 1.3
355     %%% Local Variables:
356     %%% mode: latex
357     %%% TeX-master: "ceaice"
358     %%% End:

  ViewVC Help
Powered by ViewVC 1.1.22