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

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

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


Revision 1.12 - (show annotations) (download) (as text)
Fri Aug 15 14:16:35 2008 UTC (16 years, 11 months ago) by dimitri
Branch: MAIN
CVS Tags: end_of_ceaice_monolith, HEAD
Changes since 1.11: +302 -2 lines
File MIME type: application/x-tex
Removing ceaice_adjoint.tex from ceaice.tex.  ceaice_adjoint.tex will become
basis for a second paper.  Also I suggest that there should be no appendix.
Instead a minimal description of model dynamics should be included in main
body of text but with references for everything that is already published,
including equations, etc., and only pointing out changes.  As a first step I
have brought back all the appendix to ceaice_model, then we can start cutting.

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

  ViewVC Help
Powered by ViewVC 1.1.22