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

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

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


Revision 1.2 - (show annotations) (download) (as text)
Fri Jul 25 14:57:43 2008 UTC (16 years, 11 months ago) by mlosch
Branch: MAIN
CVS Tags: end_of_ceaice_monolith, HEAD
Changes since 1.1: +12 -8 lines
File MIME type: application/x-tex
fix typo, add something about the effective conductivity in the
presence of snow.

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

  ViewVC Help
Powered by ViewVC 1.1.22