ceaice_split_version/ceaice_part1/ceaice_model.tex

\section{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
%diretized 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 counditions),
%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 Candian Arctic Archipelago (CAA).

The MITgcm sea ice model (MITsim) is based on a variant of the
viscous-plastic (VP) dynamic-thermodynamic sea-ice model of
\citet{zhang97} first introduced by \citet{hibler79, hibler80}. In
order to adapt this model to the requirements of coupled ice-ocean
simulations, many important aspects of the original code have been
modified and improved:
\begin{itemize}
\item the code 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 have been adopted: LSOR \citep{zhang97}, EVP
  \citep{hunke97, hunke01};
\item ice-ocean stress can be formulated as in \citet{hibler87};
\item ice variables can be advected by sophisticated, conservative
  advection schemes with flux limiters;
\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 compontent of the
MITgcm \citep{marshall97:_finit_volum_incom_navier_stokes, mitgcm02}.
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}

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, 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 on 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, but within the ocean according to its 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 is realized via
``ghost points'', allowing 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 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 stess 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 \emph{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
now 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
the 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}

In its original formulation the sea ice model \citep{menemenlis05}
uses simple thermodynamics following the appendix of
\citet{semtner76}. This formulation does not allow storage of heat,
that is, the heat capacity of ice is zero. Upward conductive heat flux
is parameterized assuming a linear temperature profile and a constant
ice conductivity. This type of model is often refered to as a
``zero-layer'' model. 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 from
below.  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. Contrary to
\citet{menemenlis05}, 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 \citep{fen09}. 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 \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 snowice formation (a flood-freeze
algorithm following Archimedes' principle) turns snow into ice until
the ice surface is back at $z=0$ \citep{leppaeranta83}.

The concentration $c$, effective ice thickness (ice volume per unit
area, $c\cdot{h}$), and effective snow thickness ($c\cdot{h}_{s}$) are
advected by ice velocities.
%
From the various advection scheme that are available in the MITgcm
\citep{mitgcm02}, we choose flux-limited schemes
\citep[multidimensional 2nd and 3rd-order advection scheme with flux
limiter,][]{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.
The original 2nd order central differences scheme requires additional
diffusion and is used here for comparison.

Finally, 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 thermodynamics model of
\citet{winton00} is implemented in the MITsim.
It is based on the 3-layer model of \citet{semtner76} and treats brine
content by means of enthalphy conservation. This model requires
additional state variables, namely the enthalphy of the two ice
layers, 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 enthalphy. Currently this can only be
  accomplished with a 2nd-order advection scheme with flux limiter
  \citep{roe85}.}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice_part1"
%%% End: 
1	heimbach	1.1	\section{Model Formulation}
2			\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			%discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99,
7			% kreyscher00, zhang98, hunke97}, although there are sea ice models
8			%diretized 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			%flux allowed for a B-grid (with no-slip lateral boundary counditions),
27			%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			%in the Candian Arctic Archipelago (CAA).
32	heimbach	1.1
33			The MITgcm sea ice model (MITsim) is based on a variant of the
34	mlosch	1.2	viscous-plastic (VP) dynamic-thermodynamic sea-ice model of
35			\citet{zhang97} first introduced by \citet{hibler79, hibler80}. In
36			order to adapt this model to the requirements of coupled ice-ocean
37			simulations, many important aspects of the original code have been
38			modified and improved:
39	heimbach	1.1	\begin{itemize}
40			\item the code has been rewritten for an Arakawa C-grid, both B- and
41			C-grid variants are available; the C-grid code allows for no-slip
42			and free-slip lateral boundary conditions;
43			\item two different solution methods for solving the nonlinear
44			momentum equations have been adopted: LSOR \citep{zhang97}, EVP
45	mlosch	1.2	\citep{hunke97, hunke01};
46	heimbach	1.1	\item ice-ocean stress can be formulated as in \citet{hibler87};
47	mlosch	1.2	\item ice variables can be advected by sophisticated, conservative
48			advection schemes with flux limiters;
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
53			The sea ice model is tightly coupled to the ocean compontent of the
54			MITgcm \citep{marshall97:_finit_volum_incom_navier_stokes, mitgcm02}.
55			Heat, fresh water fluxes and surface stresses are computed from the
56			atmospheric state and modified by the ice model at every time step.
57	mlosch	1.2	The remainder of this section describes the model equations and
58			details of their numerical realization. Further documentation and
59			model code can be found at \url{http://mitgcm.org}.
60
61			\subsection{Dynamics}
62			\label{sec:dynamics}
63
64			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			\citet{hunke01}. In this technique, 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			always lies on the elliptic yield curve by definition. Alternatively,
79			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			ice and snow \citep{campin08}; ice does not float on top of the
85			ocean, but within the ocean according to its buoyancy.
86
87			Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
88			the tangential velocities points lie on the boundary. For C-grids, the
89			lateral boundary condition for tangential velocities is realized via
90			``ghost points'', allowing alternatively no-slip or free-slip
91			conditions. In ocean models free-slip boundary conditions in
92			conjunction with piecewise-constant (``castellated'') coastlines have
93			been shown to reduce to no-slip boundary conditions
94			\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
95			lateral boundary conditions have not yet been studied (as far as we
96			know).
97
98			Moving sea ice exerts a surface stress on the ocean. In coupling the
99			sea-ice model to the ocean model, this stess is applied directly to
100			the surface layer of the ocean model. An alternative ocean stress
101			formulation is given by \citet{hibler87}. Rather than applying the
102			interfacial stress directly, the stress is derived from integrating
103			over the ice thickness to the bottom of the oceanic surface layer. In
104			the resulting equation for the \emph{combined} ocean-ice momentum, the
105			interfacial stress cancels and the total stress appears as the sum of
106	dimitri	1.6	wind stress and divergence of internal ice stresses \citep[see also
107	mlosch	1.2	Eq.\,2 of][]{hibler87}. While this formulation tightly embeds the
108			sea-ice into the surface layer of the ocean, its disadvantage is that
109			now the velocity in the surface layer of the ocean that is used to
110			advect ocean tracers, is really an average over the ocean surface
111			velocity and the ice velocity leading to an inconsistency as the ice
112			temperature and salinity are different from the oceanic variables.
113			Both stress coupling options are available for a direct comparison of
114			the their effects on the sea-ice solution.
115
116	dimitri	1.6	The discretization of the momentum equation is straightforward. It is
117	mlosch	1.2	similar to that of \citet{zhang98, zhang03}, but differs fundamentally
118			in the underlying grid, namely the Arakawa C-grid. The EVP model, in
119			particular, is discretized naturally on the C-grid with the diagonal
120			components of the stress tensor on the center points and the
121			off-diagonal term on the corner (or vorticity) points of the grid.
122			With this choice all derivatives are discretized as central
123			differences and very little averaging is involved. Apart from the
124			standard C-grid implementation, the original B-grid implementation of
125			\citet{zhang97} is also available as an option in the code.
126
127			\subsection{Thermodynamics}
128			\label{sec:thermodynamics}
129
130			In its original formulation the sea ice model \citep{menemenlis05}
131			uses simple thermodynamics following the appendix of
132			\citet{semtner76}. This formulation does not allow storage of heat,
133			that is, the heat capacity of ice is zero. Upward conductive heat flux
134			is parameterized assuming a linear temperature profile and a constant
135			ice conductivity. This type of model is often refered to as a
136			``zero-layer'' model. The surface heat flux is computed in a similar
137			way to that of \citet{parkinson79} and \citet{manabe79}.
138
139			The conductive heat flux depends strongly on the ice thickness $h$.
140			However, the ice thickness in the model represents a mean over a
141			potentially very heterogeneous thickness distribution. In order to
142			parameterize a sub-grid scale distribution for heat flux computations,
143			the mean ice thickness $h$ is split into seven thickness categories
144			$H_{n}$ that are equally distributed between $2h$ and a minimum
145			imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
146			for $n\in[1,7]$. The heat fluxes computed for each thickness category
147			is area-averaged to give the total heat flux \citep{hibler84}.
148
149			The atmospheric heat flux is balanced by an oceanic heat flux from
150			below. The oceanic flux is proportional to the difference between
151	mlosch	1.3	ocean surface temperature and the freezing point temperature of sea
152			water, which is a function of salinity. Contrary to
153	mlosch	1.2	\citet{menemenlis05}, this flux is not assumed to instantaneously melt
154	mlosch	1.3	or create ice, but a time scale of three days is used to relax the
155	dimitri	1.7	ocean temperature to the freezing point \citep{fen09}. While this
156	mlosch	1.3	parameterization is not new \citep[it follows the ideas of
157			e.g.,][]{mellor86, mcphee92, lohmann98, notz03}, it avoids a
158			discontinuity in the functional relationship between model variables,
159	dimitri	1.6	which is crucial for making the code differentiable for adjoint code
160	dimitri	1.7	generation \citep{fen09}.
161	mlosch	1.2	%
162			The parameterization of lateral and vertical growth of sea ice follows
163			that of \citet{hibler79, hibler80}.
164
165			On top of the ice there is a layer of snow that modifies the heat flux
166			and the albedo as in \citet{zhang98}. If enough snow accumulates so
167			that its weight submerges the ice and the snow is flooded, a simple
168			mass conserving parameterization of snowice formation (a flood-freeze
169			algorithm following Archimedes' principle) turns snow into ice until
170			the ice surface is back at $z=0$ \citep{leppaeranta83}.
171
172			The concentration $c$, effective ice thickness (ice volume per unit
173			area, $c\cdot{h}$), and effective snow thickness ($c\cdot{h}_{s}$) are
174			advected by ice velocities.
175			%
176			From the various advection scheme that are available in the MITgcm
177			\citep{mitgcm02}, we choose flux-limited schemes
178			\citep[multidimensional 2nd and 3rd-order advection scheme with flux
179			limiter,][]{roe85, hundsdorfer94} to preserve sharp gradients and
180			edges that are typical of sea ice distributions and to rule out
181			unphysical over- and undershoots (negative thickness or
182	dimitri	1.6	concentration). These schemes conserve volume and horizontal area and
183	mlosch	1.2	are unconditionally stable, so that no extra diffusion is required.
184			The original 2nd order central differences scheme requires additional
185			diffusion and is used here for comparison.
186
187			Finally, there is considerable doubt about the reliability of a
188			``zero-layer'' thermodynamic model---\citet{semtner84} found
189			significant errors in phase (one month lead) and amplitude
190			($\approx$50\%\,overestimate) in such models---, so that today many
191	dimitri	1.6	sea ice models employ more complex thermodynamics. The thermodynamics model of
192			\citet{winton00} is implemented in the MITsim.
193	mlosch	1.4	It is based on the 3-layer model of \citet{semtner76} and treats brine
194			content by means of enthalphy conservation. This model requires
195			additional state variables, namely the enthalphy of the two ice
196			layers, to be advected by ice velocities. \ml{[Jean-Michel, your
197	mlosch	1.2	turf: ]Care must be taken in advecting these quantities in order to
198			ensure conservation of enthalphy. Currently this can only be
199			accomplished with a 2nd-order advection scheme with flux limiter
200			\citep{roe85}.}
201	heimbach	1.1
202			%%% Local Variables:
203			%%% mode: latex
204	mlosch	1.2	%%% TeX-master: "ceaice_part1"
205	heimbach	1.1	%%% End: