1 |
|
% $Header$ |
2 |
|
% $Name$ |
3 |
\documentclass[12pt]{article} |
\documentclass[12pt]{article} |
4 |
|
|
5 |
\usepackage[]{graphicx} |
\usepackage[]{graphicx} |
52 |
\maketitle |
\maketitle |
53 |
|
|
54 |
\begin{abstract} |
\begin{abstract} |
55 |
Some blabla |
As part of ongoing efforts to obtain a best possible synthesis of most |
56 |
|
available, global-scale, ocean and sea ice data, a dynamic and thermodynamic |
57 |
|
sea-ice model has been coupled to the Massachusetts Institute of Technology |
58 |
|
general circulation model (MITgcm). Ice mechanics follow a viscous plastic |
59 |
|
rheology and the ice momentum equations are solved numerically using either |
60 |
|
line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic |
61 |
|
models. Ice thermodynamics are represented using either a zero-heat-capacity |
62 |
|
formulation or a two-layer formulation that conserves enthalpy. The model |
63 |
|
includes prognostic variables for snow and for sea-ice salinity. The above |
64 |
|
sea ice model components were borrowed from current-generation climate models |
65 |
|
but they were reformulated on an Arakawa C-grid in order to match the MITgcm |
66 |
|
oceanic grid and they were modified in many ways to permit efficient and |
67 |
|
accurate automatic differentiation. This paper describes the MITgcm sea ice |
68 |
|
model; it presents example Arctic and Antarctic results from a realistic, |
69 |
|
eddy-permitting, global ocean and sea-ice configuration; it compares B-grid |
70 |
|
and C-grid dynamic solvers in a regional Arctic configuration; and it presents |
71 |
|
example results from coupled ocean and sea-ice adjoint-model integrations. |
72 |
|
|
73 |
\end{abstract} |
\end{abstract} |
74 |
|
|
75 |
\section{Introduction} |
\section{Introduction} |
76 |
\label{sec:intro} |
\label{sec:intro} |
77 |
|
|
78 |
more blabla |
The availability of an adjoint model as a powerful research tool |
79 |
|
complementary to an ocean model was a major design requirement early |
80 |
\section{Model} |
on in the development of the MIT general circulation model (MITgcm) |
81 |
\label{sec:model} |
[Marshall et al. 1997a, Marotzke et al. 1999, Adcroft et al. 2002]. It |
82 |
|
was recognized that the adjoint model permitted computing the |
83 |
|
gradients of various scalar-valued model diagnostics, norms or, |
84 |
|
generally, objective functions with respect to external or independent |
85 |
|
parameters very efficiently. The information associtated with these |
86 |
|
gradients is useful in at least two major contexts. First, for state |
87 |
|
estimation problems, the objective function is the sum of squared |
88 |
|
differences between observations and model results weighted by the |
89 |
|
inverse error covariances. The gradient of such an objective function |
90 |
|
can be used to reduce this measure of model-data misfit to find the |
91 |
|
optimal model solution in a least-squares sense. Second, the |
92 |
|
objective function can be a key oceanographic quantity such as |
93 |
|
meridional heat or volume transport, ocean heat content or mean |
94 |
|
surface temperature index. In this case the gradient provides a |
95 |
|
complete set of sensitivities of this quantity to all independent |
96 |
|
variables simultaneously. These sensitivities can be used to address |
97 |
|
the cause of, say, changing net transports accurately. |
98 |
|
|
99 |
|
References to existing sea-ice adjoint models, explaining that they are either |
100 |
|
for simplified configurations, for ice-only studies, or for short-duration |
101 |
|
studies to motivate the present work. |
102 |
|
|
103 |
Traditionally, probably for historical reasons and the ease of |
Traditionally, probably for historical reasons and the ease of |
104 |
treating the Coriolis term, most standard sea-ice models are |
treating the Coriolis term, most standard sea-ice models are |
108 |
and fresh-water pose no difficulty for a B-grid sea-ice model |
and fresh-water pose no difficulty for a B-grid sea-ice model |
109 |
\citep[e.g.,][]{timmermann02a}. However, surface stress is defined at |
\citep[e.g.,][]{timmermann02a}. However, surface stress is defined at |
110 |
velocities points and thus needs to be interpolated between a B-grid |
velocities points and thus needs to be interpolated between a B-grid |
111 |
sea-ice model and a C-grid ocean model. While the smoothing implicitly |
sea-ice model and a C-grid ocean model. Smoothing implicitly |
112 |
associated with this interpolation may mask grid scale noise, it may |
associated with this interpolation may mask grid scale noise and may |
113 |
in two-way coupling lead to a computational mode as will be shown. By |
contribute to stabilizing the solution. On the other hand, by |
114 |
choosing a C-grid for the sea-ice model, we circumvene this difficulty |
smoothing the stress signals are damped which could lead to reduced |
115 |
altogether and render the stress coupling as consistent as the |
variability of the system. By choosing a C-grid for the sea-ice model, |
116 |
buoyancy coupling. |
we circumvent this difficulty altogether and render the stress |
117 |
|
coupling as consistent as the buoyancy coupling. |
118 |
|
|
119 |
A further advantage of the C-grid formulation is apparent in narrow |
A further advantage of the C-grid formulation is apparent in narrow |
120 |
straits. In the limit of only one grid cell between coasts there is no |
straits. In the limit of only one grid cell between coasts there is no |
121 |
flux allowed for a B-grid (with no-slip lateral boundary counditions), |
flux allowed for a B-grid (with no-slip lateral boundary counditions), |
122 |
whereas the C-grid formulation allows a flux of sea-ice through this |
and models have used topographies artificially widened straits to |
123 |
passage for all types of lateral boundary conditions. We (will) |
avoid this problem \citep{holloway07}. The C-grid formulation on the |
124 |
demonstrate this effect in the Candian archipelago. |
other hand allows a flux of sea-ice through narrow passages if |
125 |
|
free-slip along the boundaries is allowed. We demonstrate this effect |
126 |
|
in the Candian archipelago. |
127 |
|
|
128 |
|
Talk about problems that make the sea-ice-ocean code very sensitive and |
129 |
|
changes in the code that reduce these sensitivities. |
130 |
|
|
131 |
|
This paper describes the MITgcm sea ice |
132 |
|
model; it presents example Arctic and Antarctic results from a realistic, |
133 |
|
eddy-permitting, global ocean and sea-ice configuration; it compares B-grid |
134 |
|
and C-grid dynamic solvers in a regional Arctic configuration; and it presents |
135 |
|
example results from coupled ocean and sea-ice adjoint-model integrations. |
136 |
|
|
137 |
|
\section{Model} |
138 |
|
\label{sec:model} |
139 |
|
|
140 |
\subsection{Dynamics} |
\subsection{Dynamics} |
141 |
\label{sec:dynamics} |
\label{sec:dynamics} |
142 |
|
|
143 |
The momentum equations of the sea-ice model are standard with |
The momentum equation of the sea-ice model is |
144 |
\begin{equation} |
\begin{equation} |
145 |
\label{eq:momseaice} |
\label{eq:momseaice} |
146 |
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} + |
147 |
\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, |
\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, |
148 |
\end{equation} |
\end{equation} |
149 |
where $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$ |
where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area; |
150 |
the ice mass per unit area, $f$ the Coriolis parameter, $g$ is the |
$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector; |
151 |
gravity accelation, $\nabla\phi$ is the gradient (tilt) of the sea |
$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$ |
152 |
surface height potential beneath the ice. $\phi$ is the sum of |
directions, respectively; |
153 |
atmpheric pressure $p_{a}$ and loading due to ice and snow |
$f$ is the Coriolis parameter; |
154 |
$(m_{i}+m_{s})g$. $\vtau_{air}$ and $\vtau_{ocean}$ are the wind and |
$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses, |
155 |
ice-ocean stresses, respectively. $\vek{F}$ is the interaction force |
respectively; |
156 |
and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the |
$g$ is the gravity accelation; |
157 |
$x$, $y$, and $z$ directions. Advection of sea-ice momentum is |
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; |
158 |
neglected. The wind and ice-ocean stress terms are given by |
$\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential |
159 |
|
in response to ocean dynamics ($g\eta$) and to atmospheric pressure |
160 |
|
loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density); |
161 |
|
and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress |
162 |
|
tensor $\sigma_{ij}$. |
163 |
|
When using the rescaled vertical coordinate system, z$^\ast$, of |
164 |
|
\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice |
165 |
|
loading, $mg/\rho_{0}$. |
166 |
|
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress |
167 |
|
terms are given by |
168 |
\begin{align*} |
\begin{align*} |
169 |
\vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\ |
\vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}| |
170 |
\vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}| |
R_{air} (\vek{U}_{air} -\vek{u}), \\ |
171 |
|
\vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| |
172 |
R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ |
R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ |
173 |
\end{align*} |
\end{align*} |
174 |
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere |
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere |
175 |
and surface currents of the ocean, respectively. $C_{air/ocean}$ are |
and surface currents of the ocean, respectively; $C_{air/ocean}$ are |
176 |
air and ocean drag coefficients, $\rho_{air/ocean}$ reference |
air and ocean drag coefficients; $\rho_{air/ocean}$ are reference |
177 |
densities, and $R_{air/ocean}$ rotation matrices that act on the |
densities; and $R_{air/ocean}$ are rotation matrices that act on the |
178 |
wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence |
wind/current vectors. |
179 |
of the interal stress tensor $\sigma_{ij}$. |
|
180 |
|
For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can |
181 |
For an isotropic system this stress tensor can be related to the ice |
be related to the ice strain rate and strength by a nonlinear |
182 |
strain rate and strength by a nonlinear viscous-plastic (VP) |
viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}: |
|
constitutive law \citep{hibler79, zhang98}: |
|
183 |
\begin{equation} |
\begin{equation} |
184 |
\label{eq:vpequation} |
\label{eq:vpequation} |
185 |
\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} |
\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} |
197 |
both thickness $h$ and compactness (concentration) $c$: |
both thickness $h$ and compactness (concentration) $c$: |
198 |
\begin{equation} |
\begin{equation} |
199 |
P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]}, |
P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]}, |
200 |
\label{icestrength} |
\label{eq:icestrength} |
201 |
\end{equation} |
\end{equation} |
202 |
with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear |
with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear |
203 |
viscosities $\eta$ and $\zeta$ are functions of ice strain rate |
viscosities $\eta$ and $\zeta$ are functions of ice strain rate |
219 |
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a |
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a |
220 |
maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where |
maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where |
221 |
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress |
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress |
222 |
tensor compuation the replacement pressure $P = 2\,\Delta\zeta$ |
tensor computation the replacement pressure $P = 2\,\Delta\zeta$ |
223 |
\citep{hibler95} is used so that the stress state always lies on the |
\citep{hibler95} is used so that the stress state always lies on the |
224 |
elliptic yield curve by definition. |
elliptic yield curve by definition. |
225 |
|
|
227 |
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
228 |
\begin{equation} |
\begin{equation} |
229 |
\label{eq:etatem} |
\label{eq:etatem} |
230 |
\eta = \min(\frac{\zeta}{e^2} |
\eta = \min\left(\frac{\zeta}{e^2}, |
231 |
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
232 |
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
233 |
+4\dot{\epsilon}_{12}^2}} |
+4\dot{\epsilon}_{12}^2}}\right). |
234 |
\end{equation} |
\end{equation} |
235 |
|
|
236 |
In the current implementation, the VP-model is integrated with the |
In the current implementation, the VP-model is integrated with the |
237 |
semi-implicit line successive over relaxation (LSOR)-solver of |
semi-implicit line successive over relaxation (LSOR)-solver of |
238 |
\citet{zhang98}, which allows for long time steps that, in our case, |
\citet{zhang98}, which allows for long time steps that, in our case, |
239 |
is limited by the explicit treatment of the Coriolis term. The |
are limited by the explicit treatment of the Coriolis term. The |
240 |
explicit treatment of the Coriolis term does not represent a severe |
explicit treatment of the Coriolis term does not represent a severe |
241 |
limitation because it restricts the time step to approximately the |
limitation because it restricts the time step to approximately the |
242 |
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 |
243 |
treated explicitly. |
treated explicitly. |
244 |
|
|
245 |
\citet{hunke97}'s introduced an elastic contribution to the strain |
\citet{hunke97}'s introduced an elastic contribution to the strain |
246 |
rate elatic-viscous-plastic in order to regularize |
rate in order to regularize Eq.\refeq{vpequation} in such a way that |
247 |
Eq.\refeq{vpequation} in such a way that the resulting |
the resulting elastic-viscous-plastic (EVP) and VP models are |
248 |
elatic-viscous-plastic (EVP) and VP models are identical at steady |
identical at steady state, |
|
state, |
|
249 |
\begin{equation} |
\begin{equation} |
250 |
\label{eq:evpequation} |
\label{eq:evpequation} |
251 |
\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + |
\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + |
270 |
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension |
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension |
271 |
and shearing strain rates, $D_T = |
and shearing strain rates, $D_T = |
272 |
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = |
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = |
273 |
2\dot{\epsilon}_{12}$, respectively and using the above abbreviations, |
2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations, |
274 |
the equations can be written as: |
the equations can be written as: |
275 |
\begin{align} |
\begin{align} |
276 |
\label{eq:evpstresstensor1} |
\label{eq:evpstresstensor1} |
292 |
For details of the spatial discretization, the reader is referred to |
For details of the spatial discretization, the reader is referred to |
293 |
\citet{zhang98, zhang03}. Our discretization differs only (but |
\citet{zhang98, zhang03}. Our discretization differs only (but |
294 |
importantly) in the underlying grid, namely the Arakawa C-grid, but is |
importantly) in the underlying grid, namely the Arakawa C-grid, but is |
295 |
otherwise straightforward. The EVP model in particular is discretized |
otherwise straightforward. The EVP model, in particular, is discretized |
296 |
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 |
297 |
center points and $\sigma_{12}$ on the corner (or vorticity) points of |
center points and $\sigma_{12}$ on the corner (or vorticity) points of |
298 |
the grid. With this choice all derivatives are discretized as central |
the grid. With this choice all derivatives are discretized as central |
300 |
$P$ at vorticity points. |
$P$ at vorticity points. |
301 |
|
|
302 |
For a general curvilinear grid, one needs in principle to take metric |
For a general curvilinear grid, one needs in principle to take metric |
303 |
terms into account that arise in the transformation a curvilinear grid |
terms into account that arise in the transformation of a curvilinear |
304 |
on the sphere. However, for now we can neglect these metric terms |
grid on the sphere. For now, however, we can neglect these metric |
305 |
because they are very small on the cubed sphere grids used in this |
terms because they are very small on the \ml{[modify following |
306 |
paper; in particular, only near the edges of the cubed sphere grid, we |
section3:] cubed sphere grids used in this paper; in particular, |
307 |
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 |
308 |
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\ |
309 |
simulations. Everywhere else the coordinate system is locally nearly |
which are never covered by sea-ice in our simulations. Everywhere |
310 |
cartesian. However, for last-glacial-maximum or snowball-earth-like |
else the coordinate system is locally nearly cartesian.} However, for |
311 |
simulations the question of metric terms needs to be reconsidered. |
last-glacial-maximum or snowball-earth-like simulations the question |
312 |
Either, one includes these terms as in \citet{zhang03}, or one finds a |
of metric terms needs to be reconsidered. Either, one includes these |
313 |
vector-invariant formulation fo the sea-ice internal stress term that |
terms as in \citet{zhang03}, or one finds a vector-invariant |
314 |
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 |
315 |
of the MITgcm \citep{adcroft04:_cubed_sphere}. |
any metric terms, as it is done in the ocean dynamics of the MITgcm |
316 |
|
\citep{adcroft04:_cubed_sphere}. |
317 |
|
|
318 |
|
Lateral boundary conditions are naturally ``no-slip'' for B-grids, as |
319 |
|
the tangential velocities points lie on the boundary. For C-grids, the |
320 |
|
lateral boundary condition for tangential velocities is realized via |
321 |
|
``ghost points'', allowing alternatively no-slip or free-slip |
322 |
|
conditions. In ocean models free-slip boundary conditions in |
323 |
|
conjunction with piecewise-constant (``castellated'') coastlines have |
324 |
|
been shown to reduce in effect to no-slip boundary conditions |
325 |
|
\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of |
326 |
|
lateral boundary conditions have not yet been studied. |
327 |
|
|
328 |
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 |
329 |
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is |
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is |
342 |
temperature and salinity are different from the oceanic variables. |
temperature and salinity are different from the oceanic variables. |
343 |
|
|
344 |
Sea ice distributions are characterized by sharp gradients and edges. |
Sea ice distributions are characterized by sharp gradients and edges. |
345 |
For this reason, we employ a positive 3rd-order advection scheme |
For this reason, we employ positive, multidimensional 2nd and 3rd-order |
346 |
\citep{hundsdorfer94} for the thermodynamic variables discussed in the |
advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the |
347 |
next section. |
thermodynamic variables discussed in the next section. |
348 |
|
|
349 |
\subparagraph{boundary conditions: no-slip, free-slip, half-slip} |
\subparagraph{boundary conditions: no-slip, free-slip, half-slip} |
350 |
|
|
381 |
addition to ice-thickness and compactness (fractional area) additional |
addition to ice-thickness and compactness (fractional area) additional |
382 |
state variables to be advected by ice velocities, namely enthalphy of |
state variables to be advected by ice velocities, namely enthalphy of |
383 |
the two ice layers and the thickness of the overlying snow layer. |
the two ice layers and the thickness of the overlying snow layer. |
384 |
|
\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these |
385 |
|
quantities in order to ensure conservation of enthalphy. Currently |
386 |
|
this can only be accomplished with a 2nd-order advection scheme with |
387 |
|
flux limiter \citep{roe85}.} |
388 |
|
|
|
\section{Funnel Experiments} |
|
|
\label{sec:funnel} |
|
|
|
|
|
For a first/detailed comparison between the different variants of the |
|
|
MIT sea ice model an idealized geometry of a periodic channel, |
|
|
1000\,km long and 500\,m wide on a non-rotating plane, with converging |
|
|
walls forming a symmetric funnel and a narrow strait of 40\,km width |
|
|
is used. The horizontal resolution is 5\,km throughout the domain |
|
|
making the narrow strait 8 grid points wide. The ice model is |
|
|
initialized with a complete ice cover of 50\,cm uniform thickness. The |
|
|
ice model is driven by a constant along channel eastward ocean current |
|
|
of 25\,cm/s that does not see the walls in the domain. All other |
|
|
ice-ocean-atmosphere interactions are turned off, in particular there |
|
|
is no feedback of ice dynamics on the ocean current. All thermodynamic |
|
|
processes are turned off so that ice thickness variations are only |
|
|
caused by convergent or divergent ice flow. Ice volume (effective |
|
|
thickness) and concentration are advected with a third-order scheme |
|
|
with a flux limiter \citep{hundsdorfer94} to avoid undershoots. This |
|
|
scheme is unconditionally stable and does not require additional |
|
|
diffusion. The time step used here is 1\,h. |
|
|
|
|
|
\reffig{funnelf0} compares the dynamic fields ice concentration $c$, |
|
|
effective thickness $h_{eff} = h\cdot{c}$, and velocities $(u,v)$ for |
|
|
five different cases at steady state (after 10\,years of integration): |
|
|
\begin{description} |
|
|
\item[B-LSRns:] LSR solver with no-slip boundary conditions on a B-grid, |
|
|
\item[C-LSRns:] LSR solver with no-slip boundary conditions on a C-grid, |
|
|
\item[C-LSRfs:] LSR solver with free-slip boundary conditions on a C-grid, |
|
|
\item[C-EVPns:] EVP solver with no-slip boundary conditions on a C-grid, |
|
|
\item[C-EVPfs:] EVP solver with free-slip boundary conditions on a C-grid, |
|
|
\end{description} |
|
|
\ml{[We have not implemented the EVP solver on a B-grid.]} |
|
|
\begin{figure*}[htbp] |
|
|
%GET \includegraphics[width=\widefigwidth]{\fpath/all_086280} |
|
|
\caption{Ice concentration, effective thickness [m], and ice |
|
|
velocities [m/s] |
|
|
for 5 different numerical solutions.} |
|
|
\label{fig:funnelf0} |
|
|
\end{figure*} |
|
|
At a first glance, the solutions look similar. This is encouraging as |
|
|
the details of discretization and numerics should not affect the |
|
|
solutions to first order. In all cases the ice-ocean stress pushes the |
|
|
ice cover eastwards, where it converges in the funnel. In the narrow |
|
|
channel the ice moves quickly (nearly free drift) and leaves the |
|
|
channel as narrow band. |
|
|
|
|
|
A close look reveals interesting differences between the B- and C-grid |
|
|
results. The zonal velocity in the narrow channel is nearly the free |
|
|
drift velocity ( = ocean velocity) of 25\,cm/s for the C-grid |
|
|
solutions, regardless of the boundary conditions, while it is just |
|
|
above 20\,cm/s for the B-grid solution. The ice accelerates to |
|
|
25\,cm/s after it exits the channel. Concentrating on the solutions |
|
|
B-LSRns and C-LSRns, the ice volume (effective thickness) along the |
|
|
boundaries in the narrow channel is larger in the B-grid case although |
|
|
the ice concentration is reduces in the C-grid case. The combined |
|
|
effect leads to a larger actual ice thickness at smaller |
|
|
concentrations in the C-grid case. However, since the effective |
|
|
thickness determines the ice strength $P$ in Eq\refeq{icestrength}, |
|
|
the ice strength and thus the bulk and shear viscosities are larger in |
|
|
the B-grid case leading to more horizontal friction. This circumstance |
|
|
might explain why the no-slip boundary conditions in the B-grid case |
|
|
appear to be more effective in reducing the flow within the narrow |
|
|
channel, than in the C-grid case. Further, the viscosities are also |
|
|
sensitive to details of the velocity gradients. Via $\Delta$, these |
|
|
gradients enter the viscosities in the denominator so that large |
|
|
gradients tend to reduce the viscosities. This again favors more flow |
|
|
along the boundaries in the C-grid case: larger velocities |
|
|
(\reffig{funnelf0}) on grid points that are closer to the boundary by |
|
|
a factor $\frac{1}{2}$ than in the B-grid case because of the stagger |
|
|
nature of the C-grid lead numerically to larger tangential gradients |
|
|
across the boundary; these in turn make the viscosities smaller for |
|
|
less tangential friction and allow more tangential flow along the |
|
|
boundaries. |
|
|
|
|
|
The above argument can also be invoked to explain the small |
|
|
differences between the free-slip and no-slip solutions on the C-grid. |
|
|
Because of the non-linearities in the ice viscosities, in particular |
|
|
along the boundaries, the no-slip boundary conditions have only a small |
|
|
impact on the solution. |
|
|
|
|
|
The difference between LSR and EVP solutions is largest in the |
|
|
effective thickness and meridional velocity fields. The EVP velocity |
|
|
fields appears to be a little noisy. This noise has been address by |
|
|
\citet{hunke01}. It can be dealt with by reducing EVP's internal time |
|
|
step (increasing the number of iterations along with the computational |
|
|
cost) or by regularizing the bulk and shear viscosities. We revisit |
|
|
the latter option by reproducing some of the results of |
|
|
\citet{hunke01}, namely the experiment described in her section~4, for |
|
|
our C-grid no-slip cases: in a square domain with a few islands the |
|
|
ice model is initialized with constant ice thickness and linearly |
|
|
increasing ice concentration to the east. The model dynamics are |
|
|
forced with a constant anticyclonic ocean gyre and by variable |
|
|
atmospheric wind whose mean direction is diagnonal to the north-east |
|
|
corner of the domain; ice volume and concentration are held constant |
|
|
(no thermodynamics and no advection by ice velocity). |
|
|
\reffig{hunke01} shows the ice velocity field, its divergence, and the |
|
|
bulk viscosity $\zeta$ for the cases C-LRSns and C-EVPns, and for a |
|
|
C-EVPns case, where \citet{hunke01}'s regularization has been |
|
|
implemented; compare to Fig.\,4 in \citet{hunke01}. The regularization |
|
|
contraint limits ice strength and viscosities as a function of damping |
|
|
time scale, resolution and EVP-time step, effectively allowing the |
|
|
elastic waves to damp out more quickly \citep{hunke01}. |
|
|
\begin{figure*}[htbp] |
|
|
%GET \includegraphics[width=\widefigwidth]{\fpath/hun12days} |
|
|
\caption{Ice flow, divergence and bulk viscosities of three |
|
|
experiments with \citet{hunke01}'s test case: C-LSRns (top), |
|
|
C-EVPns (middle), and C-EVPns with damping described in |
|
|
\citet{hunke01} (bottom).} |
|
|
\label{fig:hunke01} |
|
|
\end{figure*} |
|
|
|
|
|
In the far right (``east'') side of the domain the ice concentration |
|
|
is close to one and the ice should be nearly rigid. The applied wind |
|
|
tends to push ice toward the upper right corner. Because the highly |
|
|
compact ice is confined by the boundary, it resists any further |
|
|
compression and exhibits little motion in the rigid region on the |
|
|
right hand side. The C-LSRns solution (top row) allows high |
|
|
viscosities in the rigid region suppressing nearly all flow. |
|
|
\citet{hunke01}'s regularization for the C-EVPns solution (bottom row) |
|
|
clearly suppresses the noise present in $\nabla\cdot\vek{u}$ and |
|
|
$\log_{10}\zeta$ in the |
|
|
unregularized case (middle row), at the cost of reduced viscosities. |
|
|
These reduced viscosities lead to small but finite ice velocities |
|
|
which in turn can have a strong effect on solutions in the limit of |
|
|
nearly rigid regimes (arching and blocking, not shown). |
|
389 |
|
|
390 |
\subsection{C-grid} |
\subsection{C-grid} |
391 |
\begin{itemize} |
\begin{itemize} |
432 |
\subsection{Arctic Domain with Open Boundaries} |
\subsection{Arctic Domain with Open Boundaries} |
433 |
\label{sec:arctic} |
\label{sec:arctic} |
434 |
|
|
435 |
The Arctic domain of integration is illustrated in Fig.~\ref{???}. It |
The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}. It |
436 |
is carved out from, and obtains open boundary conditions from, the |
is carved out from, and obtains open boundary conditions from, the |
437 |
global cubed-sphere configuration of the Estimating the Circulation |
global cubed-sphere configuration of the Estimating the Circulation |
438 |
and Climate of the Ocean, Phase II (ECCO2) project |
and Climate of the Ocean, Phase II (ECCO2) project |
439 |
\citet{menemenlis05}. The domain size is 420 by 384 grid boxes |
\citet{menemenlis05}. The domain size is 420 by 384 grid boxes |
440 |
horizontally with mean horizontal grid spacing of 18 km. |
horizontally with mean horizontal grid spacing of 18 km. |
441 |
|
|
442 |
|
\begin{figure} |
443 |
|
%\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}} |
444 |
|
\caption{Bathymetry of Arctic Domain.\label{fig:arctic1}} |
445 |
|
\end{figure} |
446 |
|
|
447 |
There are 50 vertical levels ranging in thickness from 10 m near the surface |
There are 50 vertical levels ranging in thickness from 10 m near the surface |
448 |
to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from |
to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from |
449 |
the National Geophysical Data Center (NGDC) 2-minute gridded global relief |
the National Geophysical Data Center (NGDC) 2-minute gridded global relief |
752 |
helpful discussions. ML thanks Elizabeth Hunke for multiple explanations. |
helpful discussions. ML thanks Elizabeth Hunke for multiple explanations. |
753 |
|
|
754 |
\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram} |
\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram} |
|
%\bibliography{journal_abrvs,seaice,genocean,maths,mixing,mitgcmuv,bib/fram} |
|
755 |
|
|
756 |
\end{document} |
\end{document} |
757 |
|
|