141 |
\begin{equation} |
\begin{equation} |
142 |
\label{eq:momseaice} |
\label{eq:momseaice} |
143 |
m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} + |
m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} + |
144 |
\vtau_{ocean} - mg \nabla{\phi(0)} + \vek{F}, |
\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, |
145 |
\end{equation} |
\end{equation} |
146 |
where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area; |
where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area; |
147 |
$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector; |
$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector; |
152 |
respectively; |
respectively; |
153 |
$g$ is the gravity accelation; |
$g$ is the gravity accelation; |
154 |
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; |
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; |
155 |
$\phi(0)$ is the sea surface height potential in response to ocean dynamics |
$\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential |
156 |
and to atmospheric pressure loading; |
in response to ocean dynamics ($g\eta$) and to atmospheric pressure |
157 |
|
loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density); |
158 |
and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress |
and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress |
159 |
tensor $\sigma_{ij}$. |
tensor $\sigma_{ij}$. |
160 |
When using the rescaled vertical coordinate system, z$^\ast$, of |
When using the rescaled vertical coordinate system, z$^\ast$, of |
161 |
\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice loading, $mg$. |
\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice |
162 |
|
loading, $mg/\rho_{0}$. |
163 |
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress |
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress |
164 |
terms are given by |
terms are given by |
165 |
\begin{align*} |
\begin{align*} |
174 |
densities; and $R_{air/ocean}$ are rotation matrices that act on the |
densities; and $R_{air/ocean}$ are rotation matrices that act on the |
175 |
wind/current vectors. |
wind/current vectors. |
176 |
|
|
177 |
For an isotropic system this stress tensor can be related to the ice |
For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can |
178 |
strain rate and strength by a nonlinear viscous-plastic (VP) |
be related to the ice strain rate and strength by a nonlinear |
179 |
constitutive law \citep{hibler79, zhang98}: |
viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}: |
180 |
\begin{equation} |
\begin{equation} |
181 |
\label{eq:vpequation} |
\label{eq:vpequation} |
182 |
\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} |
\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} |
224 |
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
225 |
\begin{equation} |
\begin{equation} |
226 |
\label{eq:etatem} |
\label{eq:etatem} |
227 |
\eta = \min(\frac{\zeta}{e^2} |
\eta = \min\left(\frac{\zeta}{e^2}, |
228 |
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
229 |
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
230 |
+4\dot{\epsilon}_{12}^2}} |
+4\dot{\epsilon}_{12}^2}}\right). |
231 |
\end{equation} |
\end{equation} |
232 |
|
|
233 |
In the current implementation, the VP-model is integrated with the |
In the current implementation, the VP-model is integrated with the |
234 |
semi-implicit line successive over relaxation (LSOR)-solver of |
semi-implicit line successive over relaxation (LSOR)-solver of |
235 |
\citet{zhang98}, which allows for long time steps that, in our case, |
\citet{zhang98}, which allows for long time steps that, in our case, |
236 |
is limited by the explicit treatment of the Coriolis term. The |
are limited by the explicit treatment of the Coriolis term. The |
237 |
explicit treatment of the Coriolis term does not represent a severe |
explicit treatment of the Coriolis term does not represent a severe |
238 |
limitation because it restricts the time step to approximately the |
limitation because it restricts the time step to approximately the |
239 |
same length as in the ocean model where the Coriolis term is also |
same length as in the ocean model where the Coriolis term is also |
240 |
treated explicitly. |
treated explicitly. |
241 |
|
|
242 |
\citet{hunke97}'s introduced an elastic contribution to the strain |
\citet{hunke97}'s introduced an elastic contribution to the strain |
243 |
rate elastic-viscous-plastic in order to regularize |
rate in order to regularize Eq.\refeq{vpequation} in such a way that |
244 |
Eq.\refeq{vpequation} in such a way that the resulting |
the resulting elastic-viscous-plastic (EVP) and VP models are |
245 |
elastic-viscous-plastic (EVP) and VP models are identical at steady |
identical at steady state, |
|
state, |
|
246 |
\begin{equation} |
\begin{equation} |
247 |
\label{eq:evpequation} |
\label{eq:evpequation} |
248 |
\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + |
\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + |
289 |
For details of the spatial discretization, the reader is referred to |
For details of the spatial discretization, the reader is referred to |
290 |
\citet{zhang98, zhang03}. Our discretization differs only (but |
\citet{zhang98, zhang03}. Our discretization differs only (but |
291 |
importantly) in the underlying grid, namely the Arakawa C-grid, but is |
importantly) in the underlying grid, namely the Arakawa C-grid, but is |
292 |
otherwise straightforward. The EVP model in particular is discretized |
otherwise straightforward. The EVP model, in particular, is discretized |
293 |
naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the |
naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the |
294 |
center points and $\sigma_{12}$ on the corner (or vorticity) points of |
center points and $\sigma_{12}$ on the corner (or vorticity) points of |
295 |
the grid. With this choice all derivatives are discretized as central |
the grid. With this choice all derivatives are discretized as central |
297 |
$P$ at vorticity points. |
$P$ at vorticity points. |
298 |
|
|
299 |
For a general curvilinear grid, one needs in principle to take metric |
For a general curvilinear grid, one needs in principle to take metric |
300 |
terms into account that arise in the transformation of a curvilinear grid |
terms into account that arise in the transformation of a curvilinear |
301 |
on the sphere. For now, however, we can neglect these metric terms |
grid on the sphere. For now, however, we can neglect these metric |
302 |
because they are very small on the cubed sphere grids used in this |
terms because they are very small on the \ml{[modify following |
303 |
paper; in particular, only near the edges of the cubed sphere grid, we |
section3:] cubed sphere grids used in this paper; in particular, |
304 |
expect them to be non-zero, but these edges are at approximately |
only near the edges of the cubed sphere grid, we expect them to be |
305 |
35\degS\ or 35\degN\ which are never covered by sea-ice in our |
non-zero, but these edges are at approximately 35\degS\ or 35\degN\ |
306 |
simulations. Everywhere else the coordinate system is locally nearly |
which are never covered by sea-ice in our simulations. Everywhere |
307 |
cartesian. However, for last-glacial-maximum or snowball-earth-like |
else the coordinate system is locally nearly cartesian.} However, for |
308 |
simulations the question of metric terms needs to be reconsidered. |
last-glacial-maximum or snowball-earth-like simulations the question |
309 |
Either, one includes these terms as in \citet{zhang03}, or one finds a |
of metric terms needs to be reconsidered. Either, one includes these |
310 |
vector-invariant formulation for the sea-ice internal stress term that |
terms as in \citet{zhang03}, or one finds a vector-invariant |
311 |
does not require any metric terms, as it is done in the ocean dynamics |
formulation for the sea-ice internal stress term that does not require |
312 |
of the MITgcm \citep{adcroft04:_cubed_sphere}. |
any metric terms, as it is done in the ocean dynamics of the MITgcm |
313 |
|
\citep{adcroft04:_cubed_sphere}. |
314 |
|
|
315 |
|
Lateral boundary conditions are naturally ``no-slip'' for B-grids, as |
316 |
|
the tangential velocities points lie on the boundary. For C-grids, the |
317 |
|
lateral boundary condition for tangential velocities is realized via |
318 |
|
``ghost points'', allowing alternatively no-slip or free-slip |
319 |
|
conditions. In ocean models free-slip boundary conditions in |
320 |
|
conjunction with piecewise-constant (``castellated'') coastlines have |
321 |
|
been shown to reduce in effect to no-slip boundary conditions |
322 |
|
\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of |
323 |
|
lateral boundary conditions have not yet been studied. |
324 |
|
|
325 |
Moving sea ice exerts a stress on the ocean which is the opposite of |
Moving sea ice exerts a stress on the ocean which is the opposite of |
326 |
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is |
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is |
339 |
temperature and salinity are different from the oceanic variables. |
temperature and salinity are different from the oceanic variables. |
340 |
|
|
341 |
Sea ice distributions are characterized by sharp gradients and edges. |
Sea ice distributions are characterized by sharp gradients and edges. |
342 |
For this reason, we employ a positive 3rd-order advection scheme |
For this reason, we employ positive, multidimensional 2nd and 3rd-order |
343 |
\citep{hundsdorfer94} for the thermodynamic variables discussed in the |
advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the |
344 |
next section. |
thermodynamic variables discussed in the next section. |
345 |
|
|
346 |
\subparagraph{boundary conditions: no-slip, free-slip, half-slip} |
\subparagraph{boundary conditions: no-slip, free-slip, half-slip} |
347 |
|
|
378 |
addition to ice-thickness and compactness (fractional area) additional |
addition to ice-thickness and compactness (fractional area) additional |
379 |
state variables to be advected by ice velocities, namely enthalphy of |
state variables to be advected by ice velocities, namely enthalphy of |
380 |
the two ice layers and the thickness of the overlying snow layer. |
the two ice layers and the thickness of the overlying snow layer. |
381 |
|
\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these |
382 |
|
quantities in order to ensure conservation of enthalphy. Currently |
383 |
|
this can only be accomplished with a 2nd-order advection scheme with |
384 |
|
flux limiter \citep{roe85}.} |
385 |
|
|
386 |
|
|
387 |
\subsection{C-grid} |
\subsection{C-grid} |