| 52 |
\maketitle |
\maketitle |
| 53 |
|
|
| 54 |
\begin{abstract} |
\begin{abstract} |
|
|
|
| 55 |
As part of ongoing efforts to obtain a best possible synthesis of most |
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 |
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 |
sea-ice model has been coupled to the Massachusetts Institute of Technology |
| 75 |
\section{Introduction} |
\section{Introduction} |
| 76 |
\label{sec:intro} |
\label{sec:intro} |
| 77 |
|
|
| 78 |
\section{Model} |
The availability of an adjoint model as a powerful research |
| 79 |
\label{sec:model} |
tool complementary to an ocean model was a major design |
| 80 |
|
requirement early on in the development of the MIT general |
| 81 |
|
circulation model (MITgcm) [Marshall et al. 1997a, |
| 82 |
|
Marotzke et al. 1999, Adcroft et al. 2002]. It was recognized |
| 83 |
|
that the adjoint permitted very efficient computation |
| 84 |
|
of gradients of various scalar-valued model diagnostics, |
| 85 |
|
norms or, generally, objective functions with respect |
| 86 |
|
to external or independent parameters. Such gradients |
| 87 |
|
arise in at least two major contexts. If the objective function |
| 88 |
|
is the sum of squared model vs. obervation differences |
| 89 |
|
weighted by e.g. the inverse error covariances, the gradient |
| 90 |
|
of the objective function can be used to optimize this measure |
| 91 |
|
of model vs. data misfit in a least-squares sense. One |
| 92 |
|
is then solving a problem of statistical state estimation. |
| 93 |
|
If the objective function is a key oceanographic quantity |
| 94 |
|
such as meridional heat or volume transport, ocean heat |
| 95 |
|
content or mean surface temperature index, the gradient |
| 96 |
|
provides a complete set of sensitivities of this quantity |
| 97 |
|
with respect to all independent variables simultaneously. |
| 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 |
| 105 |
discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99, |
discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99, |
| 106 |
kreyscher00, zhang98, hunke97}. From the perspective of coupling a |
kreyscher00, zhang98, hunke97}. From the perspective of coupling a |
| 107 |
sea ice-model to a C-grid ocean model, the exchange of fluxes of heat |
sea ice-model to a C-grid ocean model, the exchange of fluxes of heat |
| 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 |
| 111 |
sea-ice model and a C-grid ocean model. While the smoothing implicitly |
sea-ice model and a C-grid ocean model. While the smoothing implicitly |
| 112 |
associated with this interpolation may mask grid scale noise, it may |
associated with this interpolation may mask grid scale noise, it may |
| 113 |
in two-way coupling lead to a computational mode as will be shown. By |
in two-way coupling lead to a computational mode as will be shown. By |
| 114 |
choosing a C-grid for the sea-ice model, we circumvene this difficulty |
choosing a C-grid for the sea-ice model, we circumvent this difficulty |
| 115 |
altogether and render the stress coupling as consistent as the |
altogether and render the stress coupling as consistent as the |
| 116 |
buoyancy coupling. |
buoyancy coupling. |
| 117 |
|
|
| 119 |
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 |
| 120 |
flux allowed for a B-grid (with no-slip lateral boundary counditions), |
flux allowed for a B-grid (with no-slip lateral boundary counditions), |
| 121 |
whereas the C-grid formulation allows a flux of sea-ice through this |
whereas the C-grid formulation allows a flux of sea-ice through this |
| 122 |
passage for all types of lateral boundary conditions. We (will) |
passage for all types of lateral boundary conditions. We |
| 123 |
demonstrate this effect in the Candian archipelago. |
demonstrate this effect in the Candian archipelago. |
| 124 |
|
|
| 125 |
|
Talk about problems that make the sea-ice-ocean code very sensitive and |
| 126 |
|
changes in the code that reduce these sensitivities. |
| 127 |
|
|
| 128 |
|
This paper describes the MITgcm sea ice |
| 129 |
|
model; it presents example Arctic and Antarctic results from a realistic, |
| 130 |
|
eddy-permitting, global ocean and sea-ice configuration; it compares B-grid |
| 131 |
|
and C-grid dynamic solvers in a regional Arctic configuration; and it presents |
| 132 |
|
example results from coupled ocean and sea-ice adjoint-model integrations. |
| 133 |
|
|
| 134 |
|
\section{Model} |
| 135 |
|
\label{sec:model} |
| 136 |
|
|
| 137 |
\subsection{Dynamics} |
\subsection{Dynamics} |
| 138 |
\label{sec:dynamics} |
\label{sec:dynamics} |
| 139 |
|
|
| 140 |
The momentum equations of the sea-ice model are standard with |
The momentum equation of the sea-ice model is |
| 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} - m \nabla{\phi(0)} + \vek{F}, |
\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, |
| 145 |
\end{equation} |
\end{equation} |
| 146 |
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; |
| 147 |
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; |
| 148 |
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$ |
| 149 |
surface height potential beneath the ice. $\phi$ is the sum of |
directions, respectively; |
| 150 |
atmpheric pressure $p_{a}$ and loading due to ice and snow |
$f$ is the Coriolis parameter; |
| 151 |
$(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, |
| 152 |
ice-ocean stresses, respectively. $\vek{F}$ is the interaction force |
respectively; |
| 153 |
and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the |
$g$ is the gravity accelation; |
| 154 |
$x$, $y$, and $z$ directions. Advection of sea-ice momentum is |
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; |
| 155 |
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 |
| 156 |
|
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 |
| 159 |
|
tensor $\sigma_{ij}$. |
| 160 |
|
When using the rescaled vertical coordinate system, z$^\ast$, of |
| 161 |
|
\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 |
| 164 |
|
terms are given by |
| 165 |
\begin{align*} |
\begin{align*} |
| 166 |
\vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\ |
\vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}| |
| 167 |
\vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}| |
R_{air} (\vek{U}_{air} -\vek{u}), \\ |
| 168 |
|
\vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| |
| 169 |
R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ |
R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ |
| 170 |
\end{align*} |
\end{align*} |
| 171 |
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere |
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere |
| 172 |
and surface currents of the ocean, respectively. $C_{air/ocean}$ are |
and surface currents of the ocean, respectively; $C_{air/ocean}$ are |
| 173 |
air and ocean drag coefficients, $\rho_{air/ocean}$ reference |
air and ocean drag coefficients; $\rho_{air/ocean}$ are reference |
| 174 |
densities, and $R_{air/ocean}$ rotation matrices that act on the |
densities; and $R_{air/ocean}$ are rotation matrices that act on the |
| 175 |
wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence |
wind/current vectors. |
| 176 |
of the interal stress tensor $\sigma_{ij}$. |
|
| 177 |
|
For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can |
| 178 |
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 |
| 179 |
strain rate and strength by a nonlinear viscous-plastic (VP) |
viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}: |
|
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} |
| 216 |
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a |
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a |
| 217 |
maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where |
maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where |
| 218 |
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress |
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress |
| 219 |
tensor compuation the replacement pressure $P = 2\,\Delta\zeta$ |
tensor computation the replacement pressure $P = 2\,\Delta\zeta$ |
| 220 |
\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 |
| 221 |
elliptic yield curve by definition. |
elliptic yield curve by definition. |
| 222 |
|
|
| 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 elatic-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 |
elatic-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}} + |
| 267 |
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension |
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension |
| 268 |
and shearing strain rates, $D_T = |
and shearing strain rates, $D_T = |
| 269 |
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = |
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = |
| 270 |
2\dot{\epsilon}_{12}$, respectively and using the above abbreviations, |
2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations, |
| 271 |
the equations can be written as: |
the equations can be written as: |
| 272 |
\begin{align} |
\begin{align} |
| 273 |
\label{eq:evpstresstensor1} |
\label{eq:evpstresstensor1} |
| 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 a curvilinear grid |
terms into account that arise in the transformation of a curvilinear |
| 301 |
on the sphere. However, for now 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 fo 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} |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch