/[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.17 - (hide annotations) (download) (as text)
Tue Feb 10 22:53:58 2009 UTC (16 years, 5 months ago) by cnh
Branch: MAIN
Changes since 1.16: +3 -2 lines
File MIME type: application/x-tex
Reached the end. Still need to track down a few refs.

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     Sea-ice motion is driven by ice-atmosphere, ice-ocean and internal
64     stresses. The internal stresses are evaluated following a
65     viscous-plastic (VP) constitutive law with an elliptic yield curve as
66     in \citet{hibler79}. The full momentum equations for the sea-ice model
67     and the solution by line successive over-relaxation (LSOR) are
68     described in \citet{zhang97}. Alternatively, the momentum equation
69     can be solved with an elastic-viscous-plastic (EVP) solver following
70 mlosch 1.12 \citet{hunke01}. In this technique, the evolution equations for the
71 dimitri 1.6 internal stress tensor components are solved by sub-cycling the sea ice
72     momentum solver within one ocean model time step.
73 mlosch 1.2
74     In both cases, the bulk viscosities can be bounded from above (if
75     required for numerical reasons). For stress tensor computations the
76     replacement pressure \citep{hibler95} is used so that the stress state
77 cnh 1.16 always lies within the elliptic yield curve by definition. Alternatively,
78 mlosch 1.2 in the so-called truncated ellipse method (TEM) the shear viscosity is
79     capped to suppress any tensile stress \citep{hibler97, geiger98}.
80    
81     The horizontal gradient of the ocean's surface is estimated directly
82     from ocean sea surface height and pressure loading from atmosphere,
83 dimitri 1.11 ice and snow \citep{campin08}. Ice does not float on top of the
84     ocean, instead it depresses the ocean surface according to its thickness and
85     buoyancy.
86 mlosch 1.2
87 mlosch 1.15 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 cnh 1.16 lateral boundary condition for tangential velocities
90     allow alternatively no-slip or free-slip
91 mlosch 1.2 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 cnh 1.16 \citep{adcroft98:_slippery_coast}; for coupled ocean sea-ice models the effects of
95 mlosch 1.2 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 jmc 1.9 sea-ice model to the ocean model, this stress is applied directly to
100 mlosch 1.2 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 mlosch 1.15 the resulting equation for the combined ocean-ice momentum, the
105 mlosch 1.2 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 mlosch 1.15 sea ice into the surface layer of the ocean, its disadvantage is that
109 dimitri 1.11 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 mlosch 1.2 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 mlosch 1.15 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 mlosch 1.2 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 cnh 1.16 Upward conductive heat flux
131 mlosch 1.2 is parameterized assuming a linear temperature profile and a constant
132 jmc 1.9 ice conductivity. This type of model is often referred to as a
133 cnh 1.16 ``zero-layer'' model \citet{semtner76}. The surface heat flux is computed in a similar
134 mlosch 1.2 way to that of \citet{parkinson79} and \citet{manabe79}.
135    
136     The conductive heat flux depends strongly on the ice thickness $h$.
137     However, the ice thickness in the model represents a mean over a
138     potentially very heterogeneous thickness distribution. In order to
139     parameterize a sub-grid scale distribution for heat flux computations,
140     the mean ice thickness $h$ is split into seven thickness categories
141     $H_{n}$ that are equally distributed between $2h$ and a minimum
142     imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
143     for $n\in[1,7]$. The heat fluxes computed for each thickness category
144     is area-averaged to give the total heat flux \citep{hibler84}.
145    
146 mlosch 1.12 The atmospheric heat flux is balanced by an oceanic heat flux.
147     The oceanic flux is proportional to the difference between
148 mlosch 1.3 ocean surface temperature and the freezing point temperature of sea
149 cnh 1.16 water, which is a function of salinity.
150     This flux is not assumed to instantaneously melt
151 mlosch 1.3 or create ice, but a time scale of three days is used to relax the
152 mlosch 1.8 ocean temperature to the freezing point. While this
153 dimitri 1.11 parameterization is not new \citep[it follows the ideas of,
154 mlosch 1.3 e.g.,][]{mellor86, mcphee92, lohmann98, notz03}, it avoids a
155     discontinuity in the functional relationship between model variables,
156 dimitri 1.6 which is crucial for making the code differentiable for adjoint code
157 cnh 1.16 generation (see companion, part 2, paper).
158 dimitri 1.13 %\ml{[ONCE IT IS SUBMITTED, otherwise pers. communcations:]}\citep{fen09}
159 mlosch 1.2 The parameterization of lateral and vertical growth of sea ice follows
160     that of \citet{hibler79, hibler80}.
161    
162     On top of the ice there is a layer of snow that modifies the heat flux
163     and the albedo as in \citet{zhang98}. If enough snow accumulates so
164     that its weight submerges the ice and the snow is flooded, a simple
165 dimitri 1.11 mass conserving parameterization of snow ice formation (a flood-freeze
166 mlosch 1.2 algorithm following Archimedes' principle) turns snow into ice until
167     the ice surface is back at $z=0$ \citep{leppaeranta83}.
168    
169     The concentration $c$, effective ice thickness (ice volume per unit
170 dimitri 1.11 area, $c\cdot{h}$), effective snow thickness ($c\cdot{h}_{s}$), and effective
171     ice salinity (in g\,m$^{-2}$) are advected by ice velocities.
172 mlosch 1.2 %
173     From the various advection scheme that are available in the MITgcm
174 dimitri 1.11 \citep{mitgcm02}, we choose flux-limited schemes, i.e., multidimensional 2nd
175     and 3rd-order advection schemes with flux limiters
176     \citep{roe85, hundsdorfer94}, to preserve sharp gradients and
177 mlosch 1.2 edges that are typical of sea ice distributions and to rule out
178     unphysical over- and undershoots (negative thickness or
179 dimitri 1.6 concentration). These schemes conserve volume and horizontal area and
180 mlosch 1.2 are unconditionally stable, so that no extra diffusion is required.
181    
182 cnh 1.16 There is considerable doubt about the reliability of a
183 dimitri 1.11 ``zero-layer'' thermodynamic model --- \citet{semtner84} found
184 mlosch 1.2 significant errors in phase (one month lead) and amplitude
185 dimitri 1.11 ($\approx$50\%\,overestimate) in such models --- so that today many
186     sea ice models employ more complex thermodynamics. The MITgcm
187     sea ice model provides the option to use the thermodynamics model of
188     \citet{winton00}, which in turn
189     is based on the 3-layer model of \citet{semtner76} and which treats brine
190 cnh 1.16 content by means of enthalpy conservation. This scheme requires
191 jmc 1.9 additional state variables, namely the enthalpy of the two ice
192 dimitri 1.11 layers (instead of effective ice salinity), to be advected by ice velocities.
193 jmc 1.9 % \ml{[Jean-Michel, your
194     % turf: ]Care must be taken in advecting these quantities in order to
195     % ensure conservation of enthalpy. Currently this can only be
196     % accomplished with a 2nd-order advection scheme with flux limiter
197     % \citep{roe85}.}
198 dimitri 1.11 The internal sea ice temperature is inferred from ice enthalpy.
199 cnh 1.16 To avoid unphysical (negative) values for ice thickness and
200 mlosch 1.12 concentration, a positive 2nd-order advection scheme with a SuperBee
201     flux limiter \citep{roe85}
202     is used in this study to advect all sea-ice-related
203 dimitri 1.11 quantities of the \citet{winton00} thermodynamic model.
204 jmc 1.9 Because of the non-linearity of the advection scheme,
205     care must be taken in advecting these quantities: when simply using
206     ice velocity to advect enthalpy, the total energy (i.e., the volume
207     integral of enthalpy) is not conserved. Alternatively, one can advect
208     the energy content (i.e., product of ice-volume and enthalpy)
209     but then false enthalpy extrema can occur,
210     which then leads to unrealistic ice temperature.
211 mlosch 1.15 In the currently implemented solution, the sea-ice mass flux is used
212 cnh 1.17 to advect the enthalpy in order to ensure conservation of enthalpy
213 mlosch 1.15 and to prevent false enthalpy extrema.
214 heimbach 1.1
215 cnh 1.16 In \refsec{globalmodel} and \ref{sec:arcticmodel}
216     we exercise and compare several
217     of the options, which have been discussed above, intercompare
218     the impact of the different formulations (all of which are widely
219 cnh 1.17 used in sea ice modeling today) on Arctic sea ice simulation
220     \citep{prosh07:_aomipspecial}.
221 cnh 1.16 %% Got to here..... more later
222     %% Add reference to JGR special issue here.....
223 dimitri 1.11
224 heimbach 1.1 %%% Local Variables:
225     %%% mode: latex
226 mlosch 1.2 %%% TeX-master: "ceaice_part1"
227 heimbach 1.1 %%% End:

  ViewVC Help
Powered by ViewVC 1.1.22