ceaice_split_version/ceaice_part1/ceaice_model.tex

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

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; these are the most important ones:
\begin{itemize}
\item the model has been rewritten for an Arakawa~C grid, both B- and
  C-grid variants are available; the finite-volume 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{hunke01, hunke02}, have been adopted,
\item ice-ocean stress can be formulated as in \citet{hibler87} as an
  alternative to the standard method of applying ice-ocean stress
  directly,
\item ice concentration and thickness, snow, and ice salinity or enthalpy can
  be advected by sophisticated, conservative
  advection schemes with flux limiters.%, and
%ph(
% This is Ian's work, and will be reported elsewhere, I insist;
% or Ian should be co-author
%\item growth and melt parameterizations have been refined and extended
%  in order to allow for automatic differentiation of the code.
%ph)
\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, freshwater 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; and by the horizontal surface elevation gradient
of the ocean. 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}. %
Implicit solvers such as LSOR usually require capping very high
  viscosities for numerical stability reasons. Alternatively, the
  elastic-viscous-plastic (EVP) technique following \citet{hunke01}
  regularizes large viscosities by adding an extra term in the
  constitutive law that introduces  damped elastic waves. The
  EVP-solver relaxes the ice state towards the VP rheology by
  sub-cycling the evolution equations for the internal stress tensor
  components and the sea ice momentum solver within one ocean model
  time step. Neither solver requires limiting the viscosities from below (see
\refapp{dynamics} for details).

%\ml{[Reviewer 2 asks for more explanation of VP/EVP/LSOR, do we want
%  to give more?]} DM: add the darn thing to appendix and let editor cut it out.

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.  In an alternative to the
elliptic yield curve, 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).
Free-slip boundary conditions are not implemented for the B~grid.

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