ceaice_split_version/ceaice_part1/ceaice_model.tex

\section{Sea ice model formulation}
\label{sec:model}

%Traditionally, probably for historical reasons and the ease of
%treating the Coriolis term, most standard sea-ice models are
%discretized on Arakawa~B grids \citep[e.g.,][]{hibler79, harder99,
%  kreyscher00, zhang98, hunke97}, although there are sea ice models
%discretized on a C~grid \citep[e.g.,][]{ip91, tremblay97,
%  lemieux08}. %
%\ml{[there is also MI-IM, but I only found this as a reference:
%  \url{http://retro.met.no/english/r_and_d_activities/method/num_mod/MI-IM-Documentation.pdf}]}
%From the perspective of coupling a sea ice-model to a C-grid ocean
%model, the exchange of fluxes of heat and fresh-water pose no
%difficulty for a B-grid sea-ice model \citep[e.g.,][]{timmermann02a}.
%However, surface stress is defined at velocity points and thus needs
%to be interpolated between a B-grid sea-ice model and a C-grid ocean
%model. Smoothing implicitly associated with this interpolation may
%mask grid scale noise and may contribute to stabilizing the solution.
%On the other hand, by smoothing the stress signals are damped which
%could lead to reduced variability of the system. By choosing a C-grid
%for the sea-ice model, we circumvent this difficulty altogether and
%render the stress coupling as consistent as the buoyancy coupling.

%A further advantage of the C-grid formulation is apparent in narrow
%straits. In the limit of only one grid cell between coasts there is no
%flux allowed for a B-grid (with no-slip lateral boundary conditions),
%and models have used topographies with artificially widened straits to
%avoid this problem \citep{holloway07}. The C-grid formulation on the
%other hand allows a flux of sea-ice through narrow passages if
%free-slip along the boundaries is allowed. We demonstrate this effect
%in the Canadian Arctic Archipelago (CAA).

The MITgcm sea ice model is based on a variant of the
viscous-plastic (VP) dynamic-thermodynamic sea-ice model of
\citet{zhang97} first introduced by \citet{hibler79, hibler80}.
Many aspects of the original codes have been
adapted 
\begin{itemize}
\item the model has been rewritten for an Arakawa~C grid, both B- and
  C-grid variants are available; the C-grid code allows for no-slip
  and free-slip lateral boundary conditions,
\item two different solution methods for solving the nonlinear
  momentum equations, LSOR \citep{zhang97} and EVP
  \citep{hunke97, hunke01}, have been adopted,
\item ice-ocean stress can be formulated as in \citet{hibler87},
\item ice concentration and thickness, snow, and ice salinity or enthalpy can
  be advected by sophisticated, conservative
  advection schemes with flux limiters, and
\item growth and melt parameterizations have been refined and extended
  in order to allow for automatic differentiation of the code.
\end{itemize}
The sea ice model is tightly coupled to the ocean component of the
MITgcm \citep{marshall97:_hydros_quasi_hydros_nonhy,marshall97:_finit_volum_incom_navier_stokes}.
Heat, fresh water fluxes and surface stresses are computed from the
atmospheric state and modified by the ice model at every time step.
The remainder of this section describes the model equations and
details of their numerical realization. Further documentation and
model code can be found at \url{http://mitgcm.org}.

\subsection{Dynamics}
\label{sec:dynamics}

\ml{[Sergey says: a kind of repetition; need to revise previous section]}
Sea-ice motion is driven by ice-atmosphere, ice-ocean and internal
stresses. The internal stresses are evaluated following a
viscous-plastic (VP) constitutive law with an elliptic yield curve as
in \citet{hibler79}. The full momentum equations for the sea-ice model
and the solution by line successive over-relaxation (LSOR) are
described in \citet{zhang97}.  Alternatively, the momentum equation
can be solved with an elastic-viscous-plastic (EVP) solver following
\citet{hunke01}. In this technique, the evolution equations for the
internal stress tensor components are solved by sub-cycling the sea ice
momentum solver within one ocean model time step.

In both cases, the bulk viscosities can be bounded from above (if
required for numerical reasons). For stress tensor computations the
replacement pressure \citep{hibler95} is used so that the stress state
always lies within the elliptic yield curve by definition. Alternatively,
in the so-called truncated ellipse method (TEM) the shear viscosity is
capped to suppress any tensile stress \citep{hibler97, geiger98}.

The horizontal gradient of the ocean's surface is estimated directly
from ocean sea surface height and pressure loading from atmosphere,
ice and snow \citep{campin08}.  Ice does not float on top of the
ocean, instead it depresses the ocean surface according to its thickness and
buoyancy.

Lateral boundary conditions are naturally ``no-slip'' for B~grids, as
the tangential velocities points lie on the boundary. For C~grids, the
lateral boundary condition for tangential velocities
allow alternatively no-slip or free-slip
conditions. In ocean models free-slip boundary conditions in
conjunction with piecewise-constant (``castellated'') coastlines have
been shown to reduce to no-slip boundary conditions
\citep{adcroft98:_slippery_coast}; for coupled ocean sea-ice models the effects of
lateral boundary conditions have not yet been studied (as far as we
know).

Moving sea ice exerts a surface stress on the ocean. In coupling the
sea-ice model to the ocean model, this stress is applied directly to
the surface layer of the ocean model. An alternative ocean stress
formulation is given by \citet{hibler87}. Rather than applying the
interfacial stress directly, the stress is derived from integrating
over the ice thickness to the bottom of the oceanic surface layer. In
the resulting equation for the combined ocean-ice momentum, the
interfacial stress cancels and the total stress appears as the sum of
wind stress and divergence of internal ice stresses \citep[see also
Eq.\,2 of][]{hibler87}. While this formulation tightly embeds the
sea ice into the surface layer of the ocean, its disadvantage is that
the velocity in the surface layer of the ocean that is used to
advect ocean tracers is really an average over the ocean surface
velocity and the ice velocity, leading to an inconsistency as the ice
temperature and salinity are different from the oceanic variables.
Both stress coupling options are available for a direct comparison of
their effects on the sea-ice solution.

The discretization of the momentum equation is straightforward. It is
similar to that of \citet{zhang98, zhang03}, but differs fundamentally
in the underlying grid, namely the Arakawa~C grid. The EVP model, in
particular, is discretized naturally on the C~grid with the diagonal
components of the stress tensor on the center points and the
off-diagonal term on the corner (or vorticity) points of the grid.
With this choice all derivatives are discretized as central
differences and very little averaging is involved. Apart from the
standard C-grid implementation, the original B-grid implementation of
\citet{zhang97} is also available as an option in the code.

\subsection{Thermodynamics}
\label{sec:thermodynamics}

Upward conductive heat flux
is parameterized assuming a linear temperature profile and a constant
ice conductivity. This type of model is often referred to as a
``zero-layer'' model \citet{semtner76}. The surface heat flux is computed in a similar
way to that of \citet{parkinson79} and \citet{manabe79}.

The conductive heat flux depends strongly on the ice thickness $h$.
However, the ice thickness in the model represents a mean over a
potentially very heterogeneous thickness distribution.  In order to
parameterize a sub-grid scale distribution for heat flux computations,
the mean ice thickness $h$ is split into seven thickness categories
$H_{n}$ that are equally distributed between $2h$ and a minimum
imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
for $n\in[1,7]$. The heat fluxes computed for each thickness category
is area-averaged to give the total heat flux \citep{hibler84}.

The atmospheric heat flux is balanced by an oceanic heat flux.
The oceanic flux is proportional to the difference between
ocean surface temperature and the freezing point temperature of sea
water, which is a function of salinity.
This flux is not assumed to instantaneously melt
or create ice, but a time scale of three days is used to relax the
ocean temperature to the freezing point. While this
parameterization is not new \citep[it follows the ideas of,
e.g.,][]{mellor86, mcphee92, lohmann98, notz03}, it avoids a
discontinuity in the functional relationship between model variables,
which is crucial for making the code differentiable for adjoint code
generation (see companion, part 2, paper).
%\ml{[ONCE IT IS SUBMITTED, otherwise pers. communcations:]}\citep{fen09}
The parameterization of lateral and vertical growth of sea ice follows
that of \citet{hibler79, hibler80}.

On top of the ice there is a layer of snow that modifies the heat flux
and the albedo as in \citet{zhang98}.  If enough snow accumulates so
that its weight submerges the ice and the snow is flooded, a simple
mass conserving parameterization of snow ice formation (a flood-freeze
algorithm following Archimedes' principle) turns snow into ice until
the ice surface is back %at $z=0$
to sea-level \citep{leppaeranta83}.

The concentration $c$, effective ice thickness (ice volume per unit
area, $c\cdot{h}$), effective snow thickness ($c\cdot{h}_{s}$), and effective
ice salinity (in g\,m$^{-2}$) are advected by ice velocities.
%
From the various advection scheme that are available in the MITgcm
\citep{mitgcm02}, we choose flux-limited schemes, i.e., multidimensional 2nd
and 3rd-order advection schemes with flux limiters
\citep{roe85, hundsdorfer94}, to preserve sharp gradients and
edges that are typical of sea ice distributions and to rule out
unphysical over- and undershoots (negative thickness or
concentration). These schemes, conserve volume and horizontal area and
are unconditionally stable, so that no extra diffusion is required. 
%\ml{[Sergey says: any FV scheme will do this]}

There is considerable doubt about the reliability of a
``zero-layer'' thermodynamic model --- \citet{semtner84} found
significant errors in phase (one month lead) and amplitude
($\approx$50\%\,overestimate) in such models --- so that today many
sea ice models employ more complex thermodynamics. The MITgcm
sea ice model provides the option to use the thermodynamics model of
\citet{winton00}, which in turn
is based on the 3-layer model of \citet{semtner76} and which treats brine
content by means of enthalpy conservation. This scheme requires
additional state variables, namely the enthalpy of the two ice
layers (instead of effective ice salinity), to be advected by ice velocities.
%  \ml{[Jean-Michel, your
%  turf: ]Care must be taken in advecting these quantities in order to
%  ensure conservation of enthalpy. Currently this can only be
%  accomplished with a 2nd-order advection scheme with flux limiter
%  \citep{roe85}.}
The internal sea ice temperature is inferred from ice enthalpy.
To avoid unphysical (negative) values for ice thickness and
concentration, a positive 2nd-order advection scheme with a SuperBee
flux limiter \citep{roe85}
is used in this study to advect all sea-ice-related
quantities of the \citet{winton00} thermodynamic model.
Because of the non-linearity of the advection scheme,
care must be taken in advecting these quantities: when simply using 
ice velocity to advect enthalpy, the total energy (i.e., the volume
integral of enthalpy) is not conserved. Alternatively, one can advect 
the energy content (i.e., product of ice-volume and enthalpy)
but then false enthalpy extrema can occur,
which then leads to unrealistic ice temperature.
In the currently implemented solution, the sea-ice mass flux is used
to advect the enthalpy in order to ensure conservation of enthalpy
and to prevent false enthalpy extrema.

In \refsec{globalmodel} and \ref{sec:arcticmodel}
we exercise and compare several
of the options, which have been discussed above, intercompare
the impact of the different formulations (all of which are widely
used in sea ice modeling today) on Arctic sea ice simulation 
\citep{proshutinsky07:_aomip}.
%% Got to here..... more later
%% Add reference to JGR special issue here.....
%\citep{prosh07:_aomipspecial}.


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice_part1"
%%% End: 
1	cnh	1.16	\section{Sea ice model formulation}
2	heimbach	1.1	\label{sec:model}
3
4	mlosch	1.2	%Traditionally, probably for historical reasons and the ease of
5			%treating the Coriolis term, most standard sea-ice models are
6	mlosch	1.15	%discretized on Arakawa~B grids \citep[e.g.,][]{hibler79, harder99,
7	mlosch	1.2	% kreyscher00, zhang98, hunke97}, although there are sea ice models
8	mlosch	1.15	%discretized on a C~grid \citep[e.g.,][]{ip91, tremblay97,
9	mlosch	1.5	% lemieux08}. %
10	mlosch	1.2	%\ml{[there is also MI-IM, but I only found this as a reference:
11			% \url{http://retro.met.no/english/r_and_d_activities/method/num_mod/MI-IM-Documentation.pdf}]}
12			%From the perspective of coupling a sea ice-model to a C-grid ocean
13			%model, the exchange of fluxes of heat and fresh-water pose no
14			%difficulty for a B-grid sea-ice model \citep[e.g.,][]{timmermann02a}.
15			%However, surface stress is defined at velocity points and thus needs
16			%to be interpolated between a B-grid sea-ice model and a C-grid ocean
17			%model. Smoothing implicitly associated with this interpolation may
18			%mask grid scale noise and may contribute to stabilizing the solution.
19			%On the other hand, by smoothing the stress signals are damped which
20			%could lead to reduced variability of the system. By choosing a C-grid
21			%for the sea-ice model, we circumvent this difficulty altogether and
22			%render the stress coupling as consistent as the buoyancy coupling.
23
24			%A further advantage of the C-grid formulation is apparent in narrow
25			%straits. In the limit of only one grid cell between coasts there is no
26	jmc	1.9	%flux allowed for a B-grid (with no-slip lateral boundary conditions),
27	mlosch	1.2	%and models have used topographies with artificially widened straits to
28			%avoid this problem \citep{holloway07}. The C-grid formulation on the
29			%other hand allows a flux of sea-ice through narrow passages if
30			%free-slip along the boundaries is allowed. We demonstrate this effect
31	jmc	1.9	%in the Canadian Arctic Archipelago (CAA).
32	heimbach	1.1
33	dimitri	1.10	The MITgcm sea ice model is based on a variant of the
34	mlosch	1.2	viscous-plastic (VP) dynamic-thermodynamic sea-ice model of
35	cnh	1.16	\citet{zhang97} first introduced by \citet{hibler79, hibler80}.
36			Many aspects of the original codes have been
37			adapted
38	heimbach	1.1	\begin{itemize}
39	mlosch	1.15	\item the model has been rewritten for an Arakawa~C grid, both B- and
40	heimbach	1.1	C-grid variants are available; the C-grid code allows for no-slip
41	dimitri	1.11	and free-slip lateral boundary conditions,
42	heimbach	1.1	\item two different solution methods for solving the nonlinear
43	dimitri	1.11	momentum equations, LSOR \citep{zhang97} and EVP
44			\citep{hunke97, hunke01}, have been adopted,
45			\item ice-ocean stress can be formulated as in \citet{hibler87},
46			\item ice concentration and thickness, snow, and ice salinity or enthalpy can
47			be advected by sophisticated, conservative
48			advection schemes with flux limiters, and
49	heimbach	1.1	\item growth and melt parameterizations have been refined and extended
50			in order to allow for automatic differentiation of the code.
51			\end{itemize}
52	jmc	1.9	The sea ice model is tightly coupled to the ocean component of the
53	cnh	1.16	MITgcm \citep{marshall97:_hydros_quasi_hydros_nonhy,marshall97:_finit_volum_incom_navier_stokes}.
54	heimbach	1.1	Heat, fresh water fluxes and surface stresses are computed from the
55			atmospheric state and modified by the ice model at every time step.
56	mlosch	1.2	The remainder of this section describes the model equations and
57			details of their numerical realization. Further documentation and
58			model code can be found at \url{http://mitgcm.org}.
59
60			\subsection{Dynamics}
61			\label{sec:dynamics}
62
63	mlosch	1.18	\ml{[Sergey says: a kind of repetition; need to revise previous section]}
64	mlosch	1.2	Sea-ice motion is driven by ice-atmosphere, ice-ocean and internal
65			stresses. The internal stresses are evaluated following a
66			viscous-plastic (VP) constitutive law with an elliptic yield curve as
67			in \citet{hibler79}. The full momentum equations for the sea-ice model
68			and the solution by line successive over-relaxation (LSOR) are
69			described in \citet{zhang97}. Alternatively, the momentum equation
70			can be solved with an elastic-viscous-plastic (EVP) solver following
71	mlosch	1.12	\citet{hunke01}. In this technique, the evolution equations for the
72	dimitri	1.6	internal stress tensor components are solved by sub-cycling the sea ice
73			momentum solver within one ocean model time step.
74	mlosch	1.2
75			In both cases, the bulk viscosities can be bounded from above (if
76			required for numerical reasons). For stress tensor computations the
77			replacement pressure \citep{hibler95} is used so that the stress state
78	cnh	1.16	always lies within the elliptic yield curve by definition. Alternatively,
79	mlosch	1.2	in the so-called truncated ellipse method (TEM) the shear viscosity is
80			capped to suppress any tensile stress \citep{hibler97, geiger98}.
81
82			The horizontal gradient of the ocean's surface is estimated directly
83			from ocean sea surface height and pressure loading from atmosphere,
84	dimitri	1.11	ice and snow \citep{campin08}. Ice does not float on top of the
85			ocean, instead it depresses the ocean surface according to its thickness and
86			buoyancy.
87	mlosch	1.2
88	mlosch	1.15	Lateral boundary conditions are naturally ``no-slip'' for B~grids, as
89			the tangential velocities points lie on the boundary. For C~grids, the
90	cnh	1.16	lateral boundary condition for tangential velocities
91			allow alternatively no-slip or free-slip
92	mlosch	1.2	conditions. In ocean models free-slip boundary conditions in
93			conjunction with piecewise-constant (``castellated'') coastlines have
94			been shown to reduce to no-slip boundary conditions
95	cnh	1.16	\citep{adcroft98:_slippery_coast}; for coupled ocean sea-ice models the effects of
96	mlosch	1.2	lateral boundary conditions have not yet been studied (as far as we
97			know).
98
99			Moving sea ice exerts a surface stress on the ocean. In coupling the
100	jmc	1.9	sea-ice model to the ocean model, this stress is applied directly to
101	mlosch	1.2	the surface layer of the ocean model. An alternative ocean stress
102			formulation is given by \citet{hibler87}. Rather than applying the
103			interfacial stress directly, the stress is derived from integrating
104			over the ice thickness to the bottom of the oceanic surface layer. In
105	mlosch	1.15	the resulting equation for the combined ocean-ice momentum, the
106	mlosch	1.2	interfacial stress cancels and the total stress appears as the sum of
107	dimitri	1.6	wind stress and divergence of internal ice stresses \citep[see also
108	mlosch	1.2	Eq.\,2 of][]{hibler87}. While this formulation tightly embeds the
109	mlosch	1.15	sea ice into the surface layer of the ocean, its disadvantage is that
110	dimitri	1.11	the velocity in the surface layer of the ocean that is used to
111			advect ocean tracers is really an average over the ocean surface
112			velocity and the ice velocity, leading to an inconsistency as the ice
113	mlosch	1.2	temperature and salinity are different from the oceanic variables.
114			Both stress coupling options are available for a direct comparison of
115	mlosch	1.18	their effects on the sea-ice solution.
116	mlosch	1.2
117	dimitri	1.6	The discretization of the momentum equation is straightforward. It is
118	mlosch	1.2	similar to that of \citet{zhang98, zhang03}, but differs fundamentally
119	mlosch	1.15	in the underlying grid, namely the Arakawa~C grid. The EVP model, in
120			particular, is discretized naturally on the C~grid with the diagonal
121	mlosch	1.2	components of the stress tensor on the center points and the
122			off-diagonal term on the corner (or vorticity) points of the grid.
123			With this choice all derivatives are discretized as central
124			differences and very little averaging is involved. Apart from the
125			standard C-grid implementation, the original B-grid implementation of
126			\citet{zhang97} is also available as an option in the code.
127
128			\subsection{Thermodynamics}
129			\label{sec:thermodynamics}
130
131	cnh	1.16	Upward conductive heat flux
132	mlosch	1.2	is parameterized assuming a linear temperature profile and a constant
133	jmc	1.9	ice conductivity. This type of model is often referred to as a
134	cnh	1.16	``zero-layer'' model \citet{semtner76}. The surface heat flux is computed in a similar
135	mlosch	1.2	way to that of \citet{parkinson79} and \citet{manabe79}.
136
137			The conductive heat flux depends strongly on the ice thickness $h$.
138			However, the ice thickness in the model represents a mean over a
139			potentially very heterogeneous thickness distribution. In order to
140			parameterize a sub-grid scale distribution for heat flux computations,
141			the mean ice thickness $h$ is split into seven thickness categories
142			$H_{n}$ that are equally distributed between $2h$ and a minimum
143			imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
144			for $n\in[1,7]$. The heat fluxes computed for each thickness category
145			is area-averaged to give the total heat flux \citep{hibler84}.
146
147	mlosch	1.12	The atmospheric heat flux is balanced by an oceanic heat flux.
148			The oceanic flux is proportional to the difference between
149	mlosch	1.3	ocean surface temperature and the freezing point temperature of sea
150	cnh	1.16	water, which is a function of salinity.
151			This flux is not assumed to instantaneously melt
152	mlosch	1.3	or create ice, but a time scale of three days is used to relax the
153	mlosch	1.8	ocean temperature to the freezing point. While this
154	dimitri	1.11	parameterization is not new \citep[it follows the ideas of,
155	mlosch	1.3	e.g.,][]{mellor86, mcphee92, lohmann98, notz03}, it avoids a
156			discontinuity in the functional relationship between model variables,
157	dimitri	1.6	which is crucial for making the code differentiable for adjoint code
158	cnh	1.16	generation (see companion, part 2, paper).
159	dimitri	1.13	%\ml{[ONCE IT IS SUBMITTED, otherwise pers. communcations:]}\citep{fen09}
160	mlosch	1.2	The parameterization of lateral and vertical growth of sea ice follows
161			that of \citet{hibler79, hibler80}.
162
163			On top of the ice there is a layer of snow that modifies the heat flux
164			and the albedo as in \citet{zhang98}. If enough snow accumulates so
165			that its weight submerges the ice and the snow is flooded, a simple
166	dimitri	1.11	mass conserving parameterization of snow ice formation (a flood-freeze
167	mlosch	1.2	algorithm following Archimedes' principle) turns snow into ice until
168	jmc	1.19	the ice surface is back %at $z=0$
169			to sea-level \citep{leppaeranta83}.
170	mlosch	1.2
171			The concentration $c$, effective ice thickness (ice volume per unit
172	dimitri	1.11	area, $c\cdot{h}$), effective snow thickness ($c\cdot{h}_{s}$), and effective
173			ice salinity (in g\,m$^{-2}$) are advected by ice velocities.
174	mlosch	1.2	%
175			From the various advection scheme that are available in the MITgcm
176	dimitri	1.11	\citep{mitgcm02}, we choose flux-limited schemes, i.e., multidimensional 2nd
177			and 3rd-order advection schemes with flux limiters
178			\citep{roe85, hundsdorfer94}, to preserve sharp gradients and
179	mlosch	1.2	edges that are typical of sea ice distributions and to rule out
180			unphysical over- and undershoots (negative thickness or
181	mlosch	1.18	concentration). These schemes, conserve volume and horizontal area and
182			are unconditionally stable, so that no extra diffusion is required.
183			%\ml{[Sergey says: any FV scheme will do this]}
184	mlosch	1.2
185	cnh	1.16	There is considerable doubt about the reliability of a
186	dimitri	1.11	``zero-layer'' thermodynamic model --- \citet{semtner84} found
187	mlosch	1.2	significant errors in phase (one month lead) and amplitude
188	dimitri	1.11	($\approx$50\%\,overestimate) in such models --- so that today many
189			sea ice models employ more complex thermodynamics. The MITgcm
190			sea ice model provides the option to use the thermodynamics model of
191			\citet{winton00}, which in turn
192			is based on the 3-layer model of \citet{semtner76} and which treats brine
193	cnh	1.16	content by means of enthalpy conservation. This scheme requires
194	jmc	1.9	additional state variables, namely the enthalpy of the two ice
195	dimitri	1.11	layers (instead of effective ice salinity), to be advected by ice velocities.
196	jmc	1.9	% \ml{[Jean-Michel, your
197			% turf: ]Care must be taken in advecting these quantities in order to
198			% ensure conservation of enthalpy. Currently this can only be
199			% accomplished with a 2nd-order advection scheme with flux limiter
200			% \citep{roe85}.}
201	dimitri	1.11	The internal sea ice temperature is inferred from ice enthalpy.
202	cnh	1.16	To avoid unphysical (negative) values for ice thickness and
203	mlosch	1.12	concentration, a positive 2nd-order advection scheme with a SuperBee
204			flux limiter \citep{roe85}
205			is used in this study to advect all sea-ice-related
206	dimitri	1.11	quantities of the \citet{winton00} thermodynamic model.
207	jmc	1.9	Because of the non-linearity of the advection scheme,
208			care must be taken in advecting these quantities: when simply using
209			ice velocity to advect enthalpy, the total energy (i.e., the volume
210			integral of enthalpy) is not conserved. Alternatively, one can advect
211			the energy content (i.e., product of ice-volume and enthalpy)
212			but then false enthalpy extrema can occur,
213			which then leads to unrealistic ice temperature.
214	mlosch	1.15	In the currently implemented solution, the sea-ice mass flux is used
215	cnh	1.17	to advect the enthalpy in order to ensure conservation of enthalpy
216	mlosch	1.15	and to prevent false enthalpy extrema.
217	heimbach	1.1
218	cnh	1.16	In \refsec{globalmodel} and \ref{sec:arcticmodel}
219			we exercise and compare several
220			of the options, which have been discussed above, intercompare
221			the impact of the different formulations (all of which are widely
222	cnh	1.17	used in sea ice modeling today) on Arctic sea ice simulation
223	mlosch	1.18	\citep{proshutinsky07:_aomip}.
224	cnh	1.16	%% Got to here..... more later
225			%% Add reference to JGR special issue here.....
226	mlosch	1.18	%\citep{prosh07:_aomipspecial}.
227
228
229	dimitri	1.11
230	heimbach	1.1	%%% Local Variables:
231			%%% mode: latex
232	mlosch	1.2	%%% TeX-master: "ceaice_part1"
233	heimbach	1.1	%%% End: