/[MITgcm]/MITgcm_contrib/articles/ceaice_split_version/ceaice_part1/ceaice_model.tex
ViewVC logotype

Annotation of /MITgcm_contrib/articles/ceaice_split_version/ceaice_part1/ceaice_model.tex

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


Revision 1.11 - (hide annotations) (download) (as text)
Sat Dec 6 01:46:38 2008 UTC (16 years, 8 months ago) by dimitri
Branch: MAIN
Changes since 1.10: +38 -30 lines
File MIME type: application/x-tex
Suggested reorganization of results section:
ceaice_forward section split into ceaice_global and ceaice_arctic

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 jmc 1.9 %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     \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 dimitri 1.11 and free-slip lateral boundary conditions,
43 heimbach 1.1 \item two different solution methods for solving the nonlinear
44 dimitri 1.11 momentum equations, LSOR \citep{zhang97} and EVP
45     \citep{hunke97, hunke01}, have been adopted,
46     \item ice-ocean stress can be formulated as in \citet{hibler87},
47     \item ice concentration and thickness, snow, and ice salinity or enthalpy can
48     be advected by sophisticated, conservative
49     advection schemes with flux limiters, and
50 heimbach 1.1 \item growth and melt parameterizations have been refined and extended
51     in order to allow for automatic differentiation of the code.
52     \end{itemize}
53    
54 jmc 1.9 The sea ice model is tightly coupled to the ocean component of the
55 heimbach 1.1 MITgcm \citep{marshall97:_finit_volum_incom_navier_stokes, mitgcm02}.
56     Heat, fresh water fluxes and surface stresses are computed from the
57     atmospheric state and modified by the ice model at every time step.
58 mlosch 1.2 The remainder of this section describes the model equations and
59     details of their numerical realization. Further documentation and
60     model code can be found at \url{http://mitgcm.org}.
61    
62     \subsection{Dynamics}
63     \label{sec:dynamics}
64    
65     Sea-ice motion is driven by ice-atmosphere, ice-ocean and internal
66     stresses. The internal stresses are evaluated following a
67     viscous-plastic (VP) constitutive law with an elliptic yield curve as
68     in \citet{hibler79}. The full momentum equations for the sea-ice model
69     and the solution by line successive over-relaxation (LSOR) are
70     described in \citet{zhang97}. Alternatively, the momentum equation
71     can be solved with an elastic-viscous-plastic (EVP) solver following
72     \citet{hunke01}. In this technique, evolution equations for the
73 dimitri 1.6 internal stress tensor components are solved by sub-cycling the sea ice
74     momentum solver within one ocean model time step.
75 mlosch 1.2
76     In both cases, the bulk viscosities can be bounded from above (if
77     required for numerical reasons). For stress tensor computations the
78     replacement pressure \citep{hibler95} is used so that the stress state
79     always lies on the elliptic yield curve by definition. Alternatively,
80     in the so-called truncated ellipse method (TEM) the shear viscosity is
81     capped to suppress any tensile stress \citep{hibler97, geiger98}.
82    
83     The horizontal gradient of the ocean's surface is estimated directly
84     from ocean sea surface height and pressure loading from atmosphere,
85 dimitri 1.11 ice and snow \citep{campin08}. Ice does not float on top of the
86     ocean, instead it depresses the ocean surface according to its thickness and
87     buoyancy.
88 mlosch 1.2
89     Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
90     the tangential velocities points lie on the boundary. For C-grids, the
91     lateral boundary condition for tangential velocities is realized via
92     ``ghost points'', allowing alternatively no-slip or free-slip
93     conditions. In ocean models free-slip boundary conditions in
94     conjunction with piecewise-constant (``castellated'') coastlines have
95     been shown to reduce to no-slip boundary conditions
96     \citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
97     lateral boundary conditions have not yet been studied (as far as we
98     know).
99    
100     Moving sea ice exerts a surface stress on the ocean. In coupling the
101 jmc 1.9 sea-ice model to the ocean model, this stress is applied directly to
102 mlosch 1.2 the surface layer of the ocean model. An alternative ocean stress
103     formulation is given by \citet{hibler87}. Rather than applying the
104     interfacial stress directly, the stress is derived from integrating
105     over the ice thickness to the bottom of the oceanic surface layer. In
106     the resulting equation for the \emph{combined} ocean-ice momentum, the
107     interfacial stress cancels and the total stress appears as the sum of
108 dimitri 1.6 wind stress and divergence of internal ice stresses \citep[see also
109 mlosch 1.2 Eq.\,2 of][]{hibler87}. While this formulation tightly embeds the
110     sea-ice into the surface layer of the ocean, its disadvantage is that
111 dimitri 1.11 the velocity in the surface layer of the ocean that is used to
112     advect ocean tracers is really an average over the ocean surface
113     velocity and the ice velocity, leading to an inconsistency as the ice
114 mlosch 1.2 temperature and salinity are different from the oceanic variables.
115     Both stress coupling options are available for a direct comparison of
116     the their effects on the sea-ice solution.
117    
118 dimitri 1.6 The discretization of the momentum equation is straightforward. It is
119 mlosch 1.2 similar to that of \citet{zhang98, zhang03}, but differs fundamentally
120     in the underlying grid, namely the Arakawa C-grid. The EVP model, in
121     particular, is discretized naturally on the C-grid with the diagonal
122     components of the stress tensor on the center points and the
123     off-diagonal term on the corner (or vorticity) points of the grid.
124     With this choice all derivatives are discretized as central
125     differences and very little averaging is involved. Apart from the
126     standard C-grid implementation, the original B-grid implementation of
127     \citet{zhang97} is also available as an option in the code.
128    
129     \subsection{Thermodynamics}
130     \label{sec:thermodynamics}
131    
132     In its original formulation the sea ice model \citep{menemenlis05}
133     uses simple thermodynamics following the appendix of
134     \citet{semtner76}. This formulation does not allow storage of heat,
135     that is, the heat capacity of ice is zero. Upward conductive heat flux
136     is parameterized assuming a linear temperature profile and a constant
137 jmc 1.9 ice conductivity. This type of model is often referred to as a
138 mlosch 1.2 ``zero-layer'' model. The surface heat flux is computed in a similar
139     way to that of \citet{parkinson79} and \citet{manabe79}.
140    
141     The conductive heat flux depends strongly on the ice thickness $h$.
142     However, the ice thickness in the model represents a mean over a
143     potentially very heterogeneous thickness distribution. In order to
144     parameterize a sub-grid scale distribution for heat flux computations,
145     the mean ice thickness $h$ is split into seven thickness categories
146     $H_{n}$ that are equally distributed between $2h$ and a minimum
147     imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
148     for $n\in[1,7]$. The heat fluxes computed for each thickness category
149     is area-averaged to give the total heat flux \citep{hibler84}.
150    
151     The atmospheric heat flux is balanced by an oceanic heat flux from
152     below. The oceanic flux is proportional to the difference between
153 mlosch 1.3 ocean surface temperature and the freezing point temperature of sea
154     water, which is a function of salinity. Contrary to
155 mlosch 1.2 \citet{menemenlis05}, this flux is not assumed to instantaneously melt
156 mlosch 1.3 or create ice, but a time scale of three days is used to relax the
157 mlosch 1.8 ocean temperature to the freezing point. While this
158 dimitri 1.11 parameterization is not new \citep[it follows the ideas of,
159 mlosch 1.3 e.g.,][]{mellor86, mcphee92, lohmann98, notz03}, it avoids a
160     discontinuity in the functional relationship between model variables,
161 dimitri 1.6 which is crucial for making the code differentiable for adjoint code
162 dimitri 1.7 generation \citep{fen09}.
163 mlosch 1.2 %
164     The parameterization of lateral and vertical growth of sea ice follows
165     that of \citet{hibler79, hibler80}.
166    
167     On top of the ice there is a layer of snow that modifies the heat flux
168     and the albedo as in \citet{zhang98}. If enough snow accumulates so
169     that its weight submerges the ice and the snow is flooded, a simple
170 dimitri 1.11 mass conserving parameterization of snow ice formation (a flood-freeze
171 mlosch 1.2 algorithm following Archimedes' principle) turns snow into ice until
172     the ice surface is back at $z=0$ \citep{leppaeranta83}.
173    
174     The concentration $c$, effective ice thickness (ice volume per unit
175 dimitri 1.11 area, $c\cdot{h}$), effective snow thickness ($c\cdot{h}_{s}$), and effective
176     ice salinity (in g\,m$^{-2}$) are advected by ice velocities.
177 mlosch 1.2 %
178     From the various advection scheme that are available in the MITgcm
179 dimitri 1.11 \citep{mitgcm02}, we choose flux-limited schemes, i.e., multidimensional 2nd
180     and 3rd-order advection schemes with flux limiters
181     \citep{roe85, hundsdorfer94}, to preserve sharp gradients and
182 mlosch 1.2 edges that are typical of sea ice distributions and to rule out
183     unphysical over- and undershoots (negative thickness or
184 dimitri 1.6 concentration). These schemes conserve volume and horizontal area and
185 mlosch 1.2 are unconditionally stable, so that no extra diffusion is required.
186 dimitri 1.11 The original 2nd order central differences scheme, which requires additional
187     diffusion, has been retained for legacy tests and for comparison
188     with the newer schemes.
189 mlosch 1.2
190     Finally, there is considerable doubt about the reliability of a
191 dimitri 1.11 ``zero-layer'' thermodynamic model --- \citet{semtner84} found
192 mlosch 1.2 significant errors in phase (one month lead) and amplitude
193 dimitri 1.11 ($\approx$50\%\,overestimate) in such models --- so that today many
194     sea ice models employ more complex thermodynamics. The MITgcm
195     sea ice model provides the option to use the thermodynamics model of
196     \citet{winton00}, which in turn
197     is based on the 3-layer model of \citet{semtner76} and which treats brine
198 jmc 1.9 content by means of enthalpy conservation. This model requires
199     additional state variables, namely the enthalpy of the two ice
200 dimitri 1.11 layers (instead of effective ice salinity), to be advected by ice velocities.
201 jmc 1.9 % \ml{[Jean-Michel, your
202     % turf: ]Care must be taken in advecting these quantities in order to
203     % ensure conservation of enthalpy. Currently this can only be
204     % accomplished with a 2nd-order advection scheme with flux limiter
205     % \citep{roe85}.}
206 dimitri 1.11 The internal sea ice temperature is inferred from ice enthalpy.
207 jmc 1.9 Since negative ice concentration or negative ice thickness is
208     not allowed, only positive advection schemes are suitable;
209     the 2nd order advection scheme with SuperBee limiter falls into this
210 dimitri 1.11 category and has been used in this study to advect all sea-ice-related
211     quantities of the \citet{winton00} thermodynamic model.
212 jmc 1.9 Because of the non-linearity of the advection scheme,
213     care must be taken in advecting these quantities: when simply using
214     ice velocity to advect enthalpy, the total energy (i.e., the volume
215     integral of enthalpy) is not conserved. Alternatively, one can advect
216     the energy content (i.e., product of ice-volume and enthalpy)
217     but then false enthalpy extrema can occur,
218     which then leads to unrealistic ice temperature.
219     The solution currently implemented consists in
220 dimitri 1.11 using the sea ice mass flux to advect the enthalpy: this
221 jmc 1.9 ensures conservation of enthalphy and prevents
222     the occurrence of false enthalpy extrema.
223 heimbach 1.1
224 dimitri 1.11 In \refsec{globalmodel} and \ref{sec:arcticmodel} we provide example
225     applications of the MITgcm sea ice model and we exercise and compare several
226     of the options, which have been discussed hereiabove.
227    
228 heimbach 1.1 %%% Local Variables:
229     %%% mode: latex
230 mlosch 1.2 %%% TeX-master: "ceaice_part1"
231 heimbach 1.1 %%% End:

  ViewVC Help
Powered by ViewVC 1.1.22