| 3 |
\documentclass[12pt]{article} |
\documentclass[12pt]{article} |
| 4 |
|
|
| 5 |
\usepackage[]{graphicx} |
\usepackage[]{graphicx} |
| 6 |
|
%\usepackage[draft]{graphicx} |
| 7 |
\usepackage{subfigure} |
\usepackage{subfigure} |
| 8 |
|
|
| 9 |
\usepackage[round,comma]{natbib} |
\usepackage[round,comma]{natbib} |
| 13 |
\newcommand\bmmax{10} \newcommand\hmmax{10} |
\newcommand\bmmax{10} \newcommand\hmmax{10} |
| 14 |
\usepackage{bm} |
\usepackage{bm} |
| 15 |
|
|
| 16 |
|
\usepackage{url} |
| 17 |
|
|
| 18 |
% math abbreviations |
% math abbreviations |
| 19 |
\newcommand{\vek}[1]{\ensuremath{\mathbf{#1}}} |
\newcommand{\vek}[1]{\ensuremath{\mathbf{#1}}} |
| 20 |
\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}} |
\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}} |
| 45 |
% commenting scheme |
% commenting scheme |
| 46 |
\newcommand{\ml}[1]{\textsf{\slshape #1}} |
\newcommand{\ml}[1]{\textsf{\slshape #1}} |
| 47 |
|
|
| 48 |
\title{A Dynamic-Thermodynamic Sea ice Model for Ocean Climate |
\title{A Dynamic-Thermodynamic Sea Ice Model on an Arakawa C-Grid |
| 49 |
Estimation on an Arakawa C-Grid} |
for Ocean Climate Estimation and Sensitivity Studies} |
| 50 |
|
|
| 51 |
|
%Alternative title suggested by Chris Hill: |
| 52 |
|
%\title{A Sea Ice Model Designed for Ocean State Estimation and its |
| 53 |
|
% Application to Studying Sea Ice Model Dynamics in the Canadian Arctic |
| 54 |
|
% Archipelago} |
| 55 |
|
|
| 56 |
\author{Martin Losch, Dimitris Menemenlis, Patrick Heimbach, \\ |
\author{Martin Losch, Dimitris Menemenlis, Patrick Heimbach, \\ |
| 57 |
Jean-Michel Campin, and Chris Hill} |
Jean-Michel Campin, and Chris Hill} |
| 59 |
|
|
| 60 |
\maketitle |
\maketitle |
| 61 |
|
|
| 62 |
\begin{abstract} |
\input{ceaice_abstract.tex} |
| 63 |
As part of ongoing efforts to obtain a best possible synthesis of most |
|
| 64 |
available, global-scale, ocean and sea ice data, a dynamic and thermodynamic |
\input{ceaice_intro.tex} |
| 65 |
sea-ice model has been coupled to the Massachusetts Institute of Technology |
|
| 66 |
general circulation model (MITgcm). Ice mechanics follow a viscous plastic |
\input{ceaice_model.tex} |
| 67 |
rheology and the ice momentum equations are solved numerically using either |
|
| 68 |
line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic |
\input{ceaice_forward.tex} |
| 69 |
models. Ice thermodynamics are represented using either a zero-heat-capacity |
|
| 70 |
formulation or a two-layer formulation that conserves enthalpy. The model |
\input{ceaice_adjoint.tex} |
| 71 |
includes prognostic variables for snow and for sea-ice salinity. The above |
|
| 72 |
sea ice model components were borrowed from current-generation climate models |
\input{ceaice_concl.tex} |
| 73 |
but they were reformulated on an Arakawa C-grid in order to match the MITgcm |
|
| 74 |
oceanic grid and they were modified in many ways to permit efficient and |
\appendix |
| 75 |
accurate automatic differentiation. This paper describes the MITgcm sea ice |
\input{ceaice_appendix.tex} |
|
model; it presents example Arctic and Antarctic results from a realistic, |
|
|
eddy-permitting, global ocean and sea-ice configuration; it compares B-grid |
|
|
and C-grid dynamic solvers in a regional Arctic configuration; and it presents |
|
|
example results from coupled ocean and sea-ice adjoint-model integrations. |
|
|
|
|
|
\end{abstract} |
|
|
|
|
|
\section{Introduction} |
|
|
\label{sec:intro} |
|
|
|
|
|
The availability of an adjoint model as a powerful research tool |
|
|
complementary to an ocean model was a major design requirement early |
|
|
on in the development of the MIT general circulation model (MITgcm) |
|
|
[Marshall et al. 1997a, Marotzke et al. 1999, Adcroft et al. 2002]. It |
|
|
was recognized that the adjoint model permitted computing the |
|
|
gradients of various scalar-valued model diagnostics, norms or, |
|
|
generally, objective functions with respect to external or independent |
|
|
parameters very efficiently. The information associtated with these |
|
|
gradients is useful in at least two major contexts. First, for state |
|
|
estimation problems, the objective function is the sum of squared |
|
|
differences between observations and model results weighted by the |
|
|
inverse error covariances. The gradient of such an objective function |
|
|
can be used to reduce this measure of model-data misfit to find the |
|
|
optimal model solution in a least-squares sense. Second, the |
|
|
objective function can be a key oceanographic quantity such as |
|
|
meridional heat or volume transport, ocean heat content or mean |
|
|
surface temperature index. In this case the gradient provides a |
|
|
complete set of sensitivities of this quantity to all independent |
|
|
variables simultaneously. These sensitivities can be used to address |
|
|
the cause of, say, changing net transports accurately. |
|
|
|
|
|
References to existing sea-ice adjoint models, explaining that they are either |
|
|
for simplified configurations, for ice-only studies, or for short-duration |
|
|
studies to motivate the present work. |
|
|
|
|
|
Traditionally, probably for historical reasons and the ease of |
|
|
treating the Coriolis term, most standard sea-ice models are |
|
|
discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99, |
|
|
kreyscher00, zhang98, hunke97}. From the perspective of coupling a |
|
|
sea ice-model to a C-grid ocean model, the exchange of fluxes of heat |
|
|
and fresh-water pose no difficulty for a B-grid sea-ice model |
|
|
\citep[e.g.,][]{timmermann02a}. However, surface stress is defined at |
|
|
velocities points and thus needs to be interpolated between a B-grid |
|
|
sea-ice model and a C-grid ocean model. Smoothing implicitly |
|
|
associated with this interpolation may mask grid scale noise and may |
|
|
contribute to stabilizing the solution. On the other hand, by |
|
|
smoothing the stress signals are damped which could lead to reduced |
|
|
variability of the system. By choosing a C-grid for the sea-ice model, |
|
|
we circumvent this difficulty altogether and render the stress |
|
|
coupling as consistent as the buoyancy coupling. |
|
|
|
|
|
A further advantage of the C-grid formulation is apparent in narrow |
|
|
straits. In the limit of only one grid cell between coasts there is no |
|
|
flux allowed for a B-grid (with no-slip lateral boundary counditions), |
|
|
and models have used topographies artificially widened straits to |
|
|
avoid this problem \citep{holloway07}. The C-grid formulation on the |
|
|
other hand allows a flux of sea-ice through narrow passages if |
|
|
free-slip along the boundaries is allowed. We demonstrate this effect |
|
|
in the Candian archipelago. |
|
|
|
|
|
Talk about problems that make the sea-ice-ocean code very sensitive and |
|
|
changes in the code that reduce these sensitivities. |
|
|
|
|
|
This paper describes the MITgcm sea ice |
|
|
model; it presents example Arctic and Antarctic results from a realistic, |
|
|
eddy-permitting, global ocean and sea-ice configuration; it compares B-grid |
|
|
and C-grid dynamic solvers in a regional Arctic configuration; and it presents |
|
|
example results from coupled ocean and sea-ice adjoint-model integrations. |
|
|
|
|
|
\section{Model} |
|
|
\label{sec:model} |
|
|
|
|
|
\subsection{Dynamics} |
|
|
\label{sec:dynamics} |
|
|
|
|
|
The momentum equation of the sea-ice model is |
|
|
\begin{equation} |
|
|
\label{eq:momseaice} |
|
|
m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} + |
|
|
\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, |
|
|
\end{equation} |
|
|
where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area; |
|
|
$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector; |
|
|
$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$ |
|
|
directions, respectively; |
|
|
$f$ is the Coriolis parameter; |
|
|
$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses, |
|
|
respectively; |
|
|
$g$ is the gravity accelation; |
|
|
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; |
|
|
$\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential |
|
|
in response to ocean dynamics ($g\eta$) and to atmospheric pressure |
|
|
loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density); |
|
|
and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress |
|
|
tensor $\sigma_{ij}$. |
|
|
When using the rescaled vertical coordinate system, z$^\ast$, of |
|
|
\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice |
|
|
loading, $mg/\rho_{0}$. |
|
|
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress |
|
|
terms are given by |
|
|
\begin{align*} |
|
|
\vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}| |
|
|
R_{air} (\vek{U}_{air} -\vek{u}), \\ |
|
|
\vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| |
|
|
R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ |
|
|
\end{align*} |
|
|
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere |
|
|
and surface currents of the ocean, respectively; $C_{air/ocean}$ are |
|
|
air and ocean drag coefficients; $\rho_{air/ocean}$ are reference |
|
|
densities; and $R_{air/ocean}$ are rotation matrices that act on the |
|
|
wind/current vectors. |
|
|
|
|
|
For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can |
|
|
be related to the ice strain rate and strength by a nonlinear |
|
|
viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}: |
|
|
\begin{equation} |
|
|
\label{eq:vpequation} |
|
|
\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} |
|
|
+ \left[\zeta(\dot{\epsilon}_{ij},P) - |
|
|
\eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij} |
|
|
- \frac{P}{2}\delta_{ij}. |
|
|
\end{equation} |
|
|
The ice strain rate is given by |
|
|
\begin{equation*} |
|
|
\dot{\epsilon}_{ij} = \frac{1}{2}\left( |
|
|
\frac{\partial{u_{i}}}{\partial{x_{j}}} + |
|
|
\frac{\partial{u_{j}}}{\partial{x_{i}}}\right). |
|
|
\end{equation*} |
|
|
The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on |
|
|
both thickness $h$ and compactness (concentration) $c$: |
|
|
\begin{equation} |
|
|
P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]}, |
|
|
\label{eq:icestrength} |
|
|
\end{equation} |
|
|
with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear |
|
|
viscosities $\eta$ and $\zeta$ are functions of ice strain rate |
|
|
invariants and ice strength such that the principal components of the |
|
|
stress lie on an elliptical yield curve with the ratio of major to |
|
|
minor axis $e$ equal to $2$; they are given by: |
|
|
\begin{align*} |
|
|
\zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})}, |
|
|
\zeta_{\max}\right) \\ |
|
|
\eta =& \frac{\zeta}{e^2} \\ |
|
|
\intertext{with the abbreviation} |
|
|
\Delta = & \left[ |
|
|
\left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right) |
|
|
(1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 + |
|
|
2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2}) |
|
|
\right]^{-\frac{1}{2}} |
|
|
\end{align*} |
|
|
The bulk viscosities are bounded above by imposing both a minimum |
|
|
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a |
|
|
maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where |
|
|
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress |
|
|
tensor computation the replacement pressure $P = 2\,\Delta\zeta$ |
|
|
\citep{hibler95} is used so that the stress state always lies on the |
|
|
elliptic yield curve by definition. |
|
|
|
|
|
In the so-called truncated ellipse method the shear viscosity $\eta$ |
|
|
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
|
|
\begin{equation} |
|
|
\label{eq:etatem} |
|
|
\eta = \min\left(\frac{\zeta}{e^2}, |
|
|
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
|
|
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
|
|
+4\dot{\epsilon}_{12}^2}}\right). |
|
|
\end{equation} |
|
|
|
|
|
In the current implementation, the VP-model is integrated with the |
|
|
semi-implicit line successive over relaxation (LSOR)-solver of |
|
|
\citet{zhang98}, which allows for long time steps that, in our case, |
|
|
are limited by the explicit treatment of the Coriolis term. The |
|
|
explicit treatment of the Coriolis term does not represent a severe |
|
|
limitation because it restricts the time step to approximately the |
|
|
same length as in the ocean model where the Coriolis term is also |
|
|
treated explicitly. |
|
|
|
|
|
\citet{hunke97}'s introduced an elastic contribution to the strain |
|
|
rate in order to regularize Eq.\refeq{vpequation} in such a way that |
|
|
the resulting elastic-viscous-plastic (EVP) and VP models are |
|
|
identical at steady state, |
|
|
\begin{equation} |
|
|
\label{eq:evpequation} |
|
|
\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + |
|
|
\frac{1}{2\eta}\sigma_{ij} |
|
|
+ \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij} |
|
|
+ \frac{P}{4\zeta}\delta_{ij} |
|
|
= \dot{\epsilon}_{ij}. |
|
|
\end{equation} |
|
|
%In the EVP model, equations for the components of the stress tensor |
|
|
%$\sigma_{ij}$ are solved explicitly. Both model formulations will be |
|
|
%used and compared the present sea-ice model study. |
|
|
The EVP-model uses an explicit time stepping scheme with a short |
|
|
timestep. According to the recommendation of \citet{hunke97}, the |
|
|
EVP-model is stepped forward in time 120 times within the physical |
|
|
ocean model time step (although this parameter is under debate), to |
|
|
allow for elastic waves to disappear. Because the scheme does not |
|
|
require a matrix inversion it is fast in spite of the small timestep |
|
|
\citep{hunke97}. For completeness, we repeat the equations for the |
|
|
components of the stress tensor $\sigma_{1} = |
|
|
\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and |
|
|
$\sigma_{12}$. Introducing the divergence $D_D = |
|
|
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension |
|
|
and shearing strain rates, $D_T = |
|
|
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = |
|
|
2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations, |
|
|
the equations can be written as: |
|
|
\begin{align} |
|
|
\label{eq:evpstresstensor1} |
|
|
\frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} + |
|
|
\frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\ |
|
|
\label{eq:evpstresstensor2} |
|
|
\frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T} |
|
|
&= \frac{P}{2T\Delta} D_T \\ |
|
|
\label{eq:evpstresstensor12} |
|
|
\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} |
|
|
&= \frac{P}{4T\Delta} D_S |
|
|
\end{align} |
|
|
Here, the elastic parameter $E$ is redefined in terms of a damping timescale |
|
|
$T$ for elastic waves \[E=\frac{\zeta}{T}.\] |
|
|
$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and |
|
|
the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend |
|
|
$E_{0} = \frac{1}{3}$. |
|
|
|
|
|
For details of the spatial discretization, the reader is referred to |
|
|
\citet{zhang98, zhang03}. Our discretization differs only (but |
|
|
importantly) in the underlying grid, namely the Arakawa C-grid, but is |
|
|
otherwise straightforward. The EVP model, in particular, is discretized |
|
|
naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the |
|
|
center points and $\sigma_{12}$ on the corner (or vorticity) points of |
|
|
the grid. With this choice all derivatives are discretized as central |
|
|
differences and averaging is only involved in computing $\Delta$ and |
|
|
$P$ at vorticity points. |
|
|
|
|
|
For a general curvilinear grid, one needs in principle to take metric |
|
|
terms into account that arise in the transformation of a curvilinear |
|
|
grid on the sphere. For now, however, we can neglect these metric |
|
|
terms because they are very small on the \ml{[modify following |
|
|
section3:] cubed sphere grids used in this paper; in particular, |
|
|
only near the edges of the cubed sphere grid, we expect them to be |
|
|
non-zero, but these edges are at approximately 35\degS\ or 35\degN\ |
|
|
which are never covered by sea-ice in our simulations. Everywhere |
|
|
else the coordinate system is locally nearly cartesian.} However, for |
|
|
last-glacial-maximum or snowball-earth-like simulations the question |
|
|
of metric terms needs to be reconsidered. Either, one includes these |
|
|
terms as in \citet{zhang03}, or one finds a vector-invariant |
|
|
formulation for the sea-ice internal stress term that does not require |
|
|
any metric terms, as it is done in the ocean dynamics of the MITgcm |
|
|
\citep{adcroft04:_cubed_sphere}. |
|
|
|
|
|
Lateral boundary conditions are naturally ``no-slip'' for B-grids, as |
|
|
the tangential velocities points lie on the boundary. For C-grids, the |
|
|
lateral boundary condition for tangential velocities is realized via |
|
|
``ghost points'', allowing alternatively no-slip or free-slip |
|
|
conditions. In ocean models free-slip boundary conditions in |
|
|
conjunction with piecewise-constant (``castellated'') coastlines have |
|
|
been shown to reduce in effect to no-slip boundary conditions |
|
|
\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of |
|
|
lateral boundary conditions have not yet been studied. |
|
|
|
|
|
Moving sea ice exerts a stress on the ocean which is the opposite of |
|
|
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is |
|
|
applied directly to the surface layer of the ocean model. An |
|
|
alternative ocean stress formulation is given by \citet{hibler87}. |
|
|
Rather than applying $\vtau_{ocean}$ directly, the stress is derived |
|
|
from integrating over the ice thickness to the bottom of the oceanic |
|
|
surface layer. In the resulting equation for the \emph{combined} |
|
|
ocean-ice momentum, the interfacial stress cancels and the total |
|
|
stress appears as the sum of windstress and divergence of internal ice |
|
|
stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also |
|
|
Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that |
|
|
now the velocity in the surface layer of the ocean that is used to |
|
|
advect tracers, is really an average over the ocean surface |
|
|
velocity and the ice velocity leading to an inconsistency as the ice |
|
|
temperature and salinity are different from the oceanic variables. |
|
|
|
|
|
Sea ice distributions are characterized by sharp gradients and edges. |
|
|
For this reason, we employ positive, multidimensional 2nd and 3rd-order |
|
|
advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the |
|
|
thermodynamic variables discussed in the next section. |
|
|
|
|
|
\subparagraph{boundary conditions: no-slip, free-slip, half-slip} |
|
|
|
|
|
\begin{itemize} |
|
|
\item transition from existing B-Grid to C-Grid |
|
|
\item boundary conditions: no-slip, free-slip, half-slip |
|
|
\item fancy (multi dimensional) advection schemes |
|
|
\item VP vs.\ EVP \citep{hunke97} |
|
|
\item ocean stress formulation \citep{hibler87} |
|
|
\end{itemize} |
|
|
|
|
|
\subsection{Thermodynamics} |
|
|
\label{sec:thermodynamics} |
|
|
|
|
|
In the original formulation the sea ice model \citep{menemenlis05} |
|
|
uses simple thermodynamics following the appendix of |
|
|
\citet{semtner76}. This formulation does not allow storage of heat |
|
|
(heat capacity of ice is zero, and this type of model is often refered |
|
|
to as a ``zero-layer'' model). Upward heat flux is parameterized |
|
|
assuming a linear temperature profile and together with a constant ice |
|
|
conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is |
|
|
the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the |
|
|
difference between water and ice surface temperatures. The surface |
|
|
heat budget is computed in a similar way to that of |
|
|
\citet{parkinson79} and \citet{manabe79}. |
|
|
|
|
|
There is considerable doubt about the reliability of such a simple |
|
|
thermodynamic model---\citet{semtner84} found significant errors in |
|
|
phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in |
|
|
such models---, so that today many sea ice models employ more complex |
|
|
thermodynamics. A popular thermodynamics model of \citet{winton00} is |
|
|
based on the 3-layer model of \citet{semtner76} and treats brine |
|
|
content by means of enthalphy conservation. This model requires in |
|
|
addition to ice-thickness and compactness (fractional area) additional |
|
|
state variables to be advected by ice velocities, namely enthalphy of |
|
|
the two ice layers and the thickness of the overlying snow layer. |
|
|
\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these |
|
|
quantities in order to ensure conservation of enthalphy. Currently |
|
|
this can only be accomplished with a 2nd-order advection scheme with |
|
|
flux limiter \citep{roe85}.} |
|
|
|
|
|
|
|
|
\subsection{C-grid} |
|
|
\begin{itemize} |
|
|
\item no-slip vs. free-slip for both lsr and evp; |
|
|
"diagnostics" to look at and use for comparison |
|
|
\begin{itemize} |
|
|
\item ice transport through Fram Strait/Denmark Strait/Davis |
|
|
Strait/Bering strait (these are general) |
|
|
\item ice transport through narrow passages, e.g.\ Nares-Strait |
|
|
\end{itemize} |
|
|
\item compare different advection schemes (if lsr turns out to be more |
|
|
effective, then with lsr otherwise I prefer evp), eg. |
|
|
\begin{itemize} |
|
|
\item default 2nd-order with diff1=0.002 |
|
|
\item 1st-order upwind with diff1=0. |
|
|
\item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me) |
|
|
\item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.) |
|
|
\end{itemize} |
|
|
That should be enough. Here, total ice mass and location of ice edge |
|
|
is interesting. However, this comparison can be done in an idealized |
|
|
domain, may not require full Arctic Domain? |
|
|
\item |
|
|
Do a little study on the parameters of LSR and EVP |
|
|
\begin{enumerate} |
|
|
\item convergence of LSR, how many iterations do you need to get a |
|
|
true elliptic yield curve |
|
|
\item vary deltaTevp and the relaxation parameter for EVP and see when |
|
|
the EVP solution breaks down (relative to the forcing time scale). |
|
|
For this, it is essential that the evp solver gives use "stripeless" |
|
|
solutions, that is your dtevp = 1sec solutions/or 10sec solutions |
|
|
with SEAICE\_evpDampC = 615. |
|
|
\end{enumerate} |
|
|
\end{itemize} |
|
|
|
|
|
\section{Forward sensitivity experiments} |
|
|
\label{sec:forward} |
|
|
|
|
|
A second series of forward sensitivity experiments have been carried out on an |
|
|
Arctic Ocean domain with open boundaries. Once again the objective is to |
|
|
compare the old B-grid LSR dynamic solver with the new C-grid LSR and EVP |
|
|
solvers. One additional experiment is carried out to illustrate the |
|
|
differences between the two main options for sea ice thermodynamics in the MITgcm. |
|
|
|
|
|
\subsection{Arctic Domain with Open Boundaries} |
|
|
\label{sec:arctic} |
|
|
|
|
|
The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}. It |
|
|
is carved out from, and obtains open boundary conditions from, the |
|
|
global cubed-sphere configuration of the Estimating the Circulation |
|
|
and Climate of the Ocean, Phase II (ECCO2) project |
|
|
\citet{menemenlis05}. The domain size is 420 by 384 grid boxes |
|
|
horizontally with mean horizontal grid spacing of 18 km. |
|
|
|
|
|
\begin{figure} |
|
|
%\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}} |
|
|
\caption{Bathymetry of Arctic Domain.\label{fig:arctic1}} |
|
|
\end{figure} |
|
|
|
|
|
There are 50 vertical levels ranging in thickness from 10 m near the surface |
|
|
to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from |
|
|
the National Geophysical Data Center (NGDC) 2-minute gridded global relief |
|
|
data (ETOPO2) and the model employs the partial-cell formulation of |
|
|
\citet{adcroft97:_shaved_cells}, which permits accurate representation of the bathymetry. The |
|
|
model is integrated in a volume-conserving configuration using a finite volume |
|
|
discretization with C-grid staggering of the prognostic variables. In the |
|
|
ocean, the non-linear equation of state of \citet{jackett95}. The ocean model is |
|
|
coupled to a sea-ice model described hereinabove. |
|
|
|
|
|
This particular ECCO2 simulation is initialized from rest using the |
|
|
January temperature and salinity distribution from the World Ocean |
|
|
Atlas 2001 (WOA01) [Conkright et al., 2002] and it is integrated for |
|
|
32 years prior to the 1996--2001 period discussed in the study. Surface |
|
|
boundary conditions are from the National Centers for Environmental |
|
|
Prediction and the National Center for Atmospheric Research |
|
|
(NCEP/NCAR) atmospheric reanalysis [Kistler et al., 2001]. Six-hourly |
|
|
surface winds, temperature, humidity, downward short- and long-wave |
|
|
radiations, and precipitation are converted to heat, freshwater, and |
|
|
wind stress fluxes using the \citet{large81, large82} bulk formulae. |
|
|
Shortwave radiation decays exponentially as per Paulson and Simpson |
|
|
[1977]. Additionally the time-mean river run-off from Large and Nurser |
|
|
[2001] is applied and there is a relaxation to the monthly-mean |
|
|
climatological sea surface salinity values from WOA01 with a |
|
|
relaxation time scale of 3 months. Vertical mixing follows |
|
|
\citet{large94} with background vertical diffusivity of |
|
|
$1.5\times10^{-5}\text{\,m$^{2}$\,s$^{-1}$}$ and viscosity of |
|
|
$10^{-3}\text{\,m$^{2}$\,s$^{-1}$}$. A third order, direct-space-time |
|
|
advection scheme with flux limiter is employed \citep{hundsdorfer94} |
|
|
and there is no explicit horizontal diffusivity. Horizontal viscosity |
|
|
follows \citet{lei96} but |
|
|
modified to sense the divergent flow as per Fox-Kemper and Menemenlis |
|
|
[in press]. Shortwave radiation decays exponentially as per Paulson |
|
|
and Simpson [1977]. Additionally, the time-mean runoff of Large and |
|
|
Nurser [2001] is applied near the coastline and, where there is open |
|
|
water, there is a relaxation to monthly-mean WOA01 sea surface |
|
|
salinity with a time constant of 45 days. |
|
|
|
|
|
Open water, dry |
|
|
ice, wet ice, dry snow, and wet snow albedo are, respectively, 0.15, 0.85, |
|
|
0.76, 0.94, and 0.8. |
|
|
|
|
|
\begin{itemize} |
|
|
\item Configuration |
|
|
\item OBCS from cube |
|
|
\item forcing |
|
|
\item 1/2 and full resolution |
|
|
\item with a few JFM figs from C-grid LSR no slip |
|
|
ice transport through Canadian Archipelago |
|
|
thickness distribution |
|
|
ice velocity and transport |
|
|
\end{itemize} |
|
|
|
|
|
\begin{itemize} |
|
|
\item Arctic configuration |
|
|
\item ice transport through straits and near boundaries |
|
|
\item focus on narrow straits in the Canadian Archipelago |
|
|
\end{itemize} |
|
|
|
|
|
\begin{itemize} |
|
|
\item B-grid LSR no-slip |
|
|
\item C-grid LSR no-slip |
|
|
\item C-grid LSR slip |
|
|
\item C-grid EVP no-slip |
|
|
\item C-grid EVP slip |
|
|
\item C-grid LSR + TEM (truncated ellipse method, no tensile stress, new flag) |
|
|
\item C-grid LSR no-slip + Winton |
|
|
\item speed-performance-accuracy (small) |
|
|
ice transport through Canadian Archipelago differences |
|
|
thickness distribution differences |
|
|
ice velocity and transport differences |
|
|
\end{itemize} |
|
|
|
|
|
We anticipate small differences between the different models due to: |
|
|
\begin{itemize} |
|
|
\item advection schemes: along the ice-edge and regions with large |
|
|
gradients |
|
|
\item C-grid: less transport through narrow straits for no slip |
|
|
conditons, more for free slip |
|
|
\item VP vs.\ EVP: speed performance, accuracy? |
|
|
\item ocean stress: different water mass properties beneath the ice |
|
|
\end{itemize} |
|
|
|
|
|
\section{Adjoint sensiivities of the MITsim} |
|
|
|
|
|
\subsection{The adjoint of MITsim} |
|
|
|
|
|
The ability to generate tangent linear and adjoint model components |
|
|
of the MITsim has been a main design task. |
|
|
For the ocean the adjoint capability has proven to be an |
|
|
invaluable tool for sensitivity analysis as well as state estimation. |
|
|
In short, the adjoint enables very efficient computation of the gradient |
|
|
of scalar-valued model diagnostics (called cost function or objective function) |
|
|
with respect to many model "variables". |
|
|
These variables can be two- or three-dimensional fields of initial |
|
|
conditions, model parameters such as mixing coefficients, or |
|
|
time-varying surface or lateral (open) boundary conditions. |
|
|
When combined, these variables span a potentially high-dimensional |
|
|
(e.g. O(10$^8$)) so-called control space. Performing parameter perturbations |
|
|
to assess model sensitivities quickly becomes prohibitive at these scales. |
|
|
Alternatively, (time-varying) sensitivities of the objective function |
|
|
to any element of the control space can be computed very efficiently in |
|
|
one single adjoint |
|
|
model integration, provided an efficient adjoint model is available. |
|
|
|
|
|
[REFERENCES] |
|
|
|
|
|
|
|
|
The adjoint operator (ADM) is the transpose of the tangent linear operator (TLM) |
|
|
of the full (in general nonlinear) forward model, i.e. the MITsim. |
|
|
The TLM maps perturbations of elements of the control space |
|
|
(e.g. initial ice thickness distribution) |
|
|
via the model Jacobian |
|
|
to a perturbation in the objective function |
|
|
(e.g. sea-ice export at the end of the integration interval). |
|
|
\textit{Tangent} linearity ensures that the derivatives are evaluated |
|
|
with respect to the underlying model trajectory at each point in time. |
|
|
This is crucial for nonlinear trajectories and the presence of different |
|
|
regimes (e.g. effect of the seaice growth term at or away from the |
|
|
freezing point of the ocean surface). |
|
|
Ensuring tangent linearity can be easily achieved by integrating |
|
|
the full model in sync with the TLM to provide the underlying model state. |
|
|
Ensuring \textit{tangent} adjoints is equally crucial, but much more |
|
|
difficult to achieve because of the reverse nature of the integration: |
|
|
the adjoint accumulates sensitivities backward in time, |
|
|
starting from a unit perturbation of the objective function. |
|
|
The adjoint model requires the model state in reverse order. |
|
|
This presents one of the major complications in deriving an |
|
|
exact, i.e. \textit{tangent} adjoint model. |
|
|
|
|
|
Following closely the development and maintenance of TLM and ADM |
|
|
components of the MITgcm we have relied heavily on the |
|
|
autmomatic differentiation (AD) tool |
|
|
"Transformation of Algorithms in Fortran" (TAF) |
|
|
developed by Fastopt (Giering and Kaminski, 1998) |
|
|
to derive TLM and ADM code of the MITsim. |
|
|
Briefly, the nonlinear parent model is fed to the AD tool which produces |
|
|
derivative code for the specified control space and objective function. |
|
|
Following this approach has (apart from its evident success) |
|
|
several advantages: |
|
|
(1) the adjoint model is the exact adjoint operator of the parent model, |
|
|
(2) the adjoint model can be kept up to date with respect to ongoing |
|
|
development of the parent model, and adjustments to the parent model |
|
|
to extend the automatically generated adjoint are incremental changes |
|
|
only, rather than extensive re-developments, |
|
|
(3) the parallel structure of the parent model is preserved |
|
|
by the adjoint model, ensuring efficient use in high performance |
|
|
computing environments. |
|
|
|
|
|
Some initial code adjustments are required to support dependency analysis |
|
|
of the flow reversal and certain language limitations which may lead |
|
|
to irreducible flow graphs (e.g. GOTO statements). |
|
|
The problem of providing the required model state in reverse order |
|
|
at the time of evaluating nonlinear or conditional |
|
|
derivatives is solved via balancing |
|
|
storing vs. recomputation of the model state in a multi-level |
|
|
checkpointing loop. |
|
|
Again, an initial code adjustment is required to support TAFs |
|
|
checkpointing capability. |
|
|
The code adjustments are sufficiently simple so as not to cause |
|
|
major limitations to the full nonlinear parent model. |
|
|
Once in place, an adjoint model of a new model configuration |
|
|
may be derived in about 10 minutes. |
|
|
|
|
|
[HIGHLIGHT COUPLED NATURE OF THE ADJOINT!] |
|
|
|
|
|
\subsection{Special considerations} |
|
|
|
|
|
* growth term(?) |
|
|
|
|
|
* small active denominators |
|
|
|
|
|
* dynamic solver (implicit function theorem) |
|
|
|
|
|
* approximate adjoints |
|
|
|
|
|
|
|
|
\subsection{An example: sensitivities of sea-ice export through Fram Strait} |
|
|
|
|
|
We demonstrate the power of the adjoint method |
|
|
in the context of investigating sea-ice export sensitivities through Fram Strait |
|
|
(for details of this study see Heimbach et al., 2007). |
|
|
%\citep[for details of this study see][]{heimbach07}. %Heimbach et al., 2007). |
|
|
The domain chosen is a coarsened version of the Arctic face of the |
|
|
high-resolution cubed-sphere configuration of the ECCO2 project |
|
|
\citep[see][]{menemenlis05}. It covers the entire Arctic, |
|
|
extends into the North Pacific such as to cover the entire |
|
|
ice-covered regions, and comprises parts of the North Atlantic |
|
|
down to XXN to enable analysis of remote influences of the |
|
|
North Atlantic current to sea-ice variability and export. |
|
|
The horizontal resolution varies between XX and YY km |
|
|
with 50 unevenly spaced vertical levels. |
|
|
The adjoint models run efficiently on 80 processors |
|
|
(benchmarks have been performed both on an SGI Altix as well as an |
|
|
IBM SP5 at NASA/ARC). |
|
|
|
|
|
Following a 1-year spinup, the model has been integrated for four |
|
|
years between 1992 and 1995. It is forced using realistic 6-hourly |
|
|
NCEP/NCAR atmospheric state variables. Over the open ocean these are |
|
|
converted into air-sea fluxes via the bulk formulae of |
|
|
\citet{large04}. Derivation of air-sea fluxes in the presence of |
|
|
sea-ice is handled by the ice model as described in \refsec{model}. |
|
|
The objective function chosen is sea-ice export through Fram Strait |
|
|
computed for December 1995. The adjoint model computes sensitivities |
|
|
to sea-ice export back in time from 1995 to 1992 along this |
|
|
trajectory. In principle all adjoint model variable (i.e., Lagrange |
|
|
multipliers) of the coupled ocean/sea-ice model are available to |
|
|
analyze the transient sensitivity behaviour of the ocean and sea-ice |
|
|
state. Over the open ocean, the adjoint of the bulk formula scheme |
|
|
computes sensitivities to the time-varying atmospheric state. Over |
|
|
ice-covered parts, the sea-ice adjoint converts surface ocean |
|
|
sensitivities to atmospheric sensitivities. |
|
|
|
|
|
\reffig{4yradjheff}(a--d) depict sensitivities of sea-ice export |
|
|
through Fram Strait in December 1995 to changes in sea-ice thickness |
|
|
12, 24, 36, 48 months back in time. Corresponding sensitivities to |
|
|
ocean surface temperature are depicted in |
|
|
\reffig{4yradjthetalev1}(a--d). The main characteristics is |
|
|
consistency with expected advection of sea-ice over the relevant time |
|
|
scales considered. The general positive pattern means that an |
|
|
increase in sea-ice thickness at location $(x,y)$ and time $t$ will |
|
|
increase sea-ice export through Fram Strait at time $T_e$. Largest |
|
|
distances from Fram Strait indicate fastest sea-ice advection over the |
|
|
time span considered. The ice thickness sensitivities are in close |
|
|
correspondence to ocean surface sentivitites, but of opposite sign. |
|
|
An increase in temperature will incur ice melting, decrease in ice |
|
|
thickness, and therefore decrease in sea-ice export at time $T_e$. |
|
|
|
|
|
The picture is fundamentally different and much more complex |
|
|
for sensitivities to ocean temperatures away from the surface. |
|
|
\reffig{4yradjthetalev10??}(a--d) depicts ice export sensitivities to |
|
|
temperatures at roughly 400 m depth. |
|
|
Primary features are the effect of the heat transport of the North |
|
|
Atlantic current which feeds into the West Spitsbergen current, |
|
|
the circulation around Svalbard, and ... |
|
|
|
|
|
\begin{figure}[t!] |
|
|
\centerline{ |
|
|
\subfigure[{\footnotesize -12 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim072_cmax2.0E+02.eps}} |
|
|
%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf} |
|
|
% |
|
|
\subfigure[{\footnotesize -24 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim145_cmax2.0E+02.eps}} |
|
|
} |
|
|
|
|
|
\centerline{ |
|
|
\subfigure[{\footnotesize |
|
|
-36 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim218_cmax2.0E+02.eps}} |
|
|
% |
|
|
\subfigure[{\footnotesize |
|
|
-48 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim292_cmax2.0E+02.eps}} |
|
|
} |
|
|
\caption{Sensitivity of sea-ice export through Fram Strait in December 2005 to |
|
|
sea-ice thickness at various prior times. |
|
|
\label{fig:4yradjheff}} |
|
|
\end{figure} |
|
|
|
|
|
|
|
|
\begin{figure}[t!] |
|
|
\centerline{ |
|
|
\subfigure[{\footnotesize -12 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim072_cmax5.0E+01.eps}} |
|
|
%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf} |
|
|
% |
|
|
\subfigure[{\footnotesize -24 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim145_cmax5.0E+01.eps}} |
|
|
} |
|
|
|
|
|
\centerline{ |
|
|
\subfigure[{\footnotesize |
|
|
-36 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim218_cmax5.0E+01.eps}} |
|
|
% |
|
|
\subfigure[{\footnotesize |
|
|
-48 months}] |
|
|
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim292_cmax5.0E+01.eps}} |
|
|
} |
|
|
\caption{Same as \reffig{4yradjheff} but for sea surface temperature |
|
|
\label{fig:4yradjthetalev1}} |
|
|
\end{figure} |
|
|
|
|
|
|
|
|
|
|
|
\section{Discussion and conclusion} |
|
|
\label{sec:concl} |
|
|
|
|
|
The story of the paper could be: |
|
|
B-grid ice model + C-grid ocean model does not work properly for us, |
|
|
therefore C-grid ice model with advantages: |
|
|
\begin{enumerate} |
|
|
\item stress coupling simpler (no interpolation required) |
|
|
\item different boundary conditions |
|
|
\item advection schemes carry over trivially from main code |
|
|
\end{enumerate} |
|
|
LSR/EVP solutions are similar with appropriate bcs, evp parameters as |
|
|
a function of forcing time scale (when does VP solution break |
|
|
down). Same for LSR solver, provided that it works (o: |
|
|
Which scheme is more efficient for the resolution/time stepping |
|
|
parameters that we use here. What about other resolutions? |
|
| 76 |
|
|
| 77 |
\paragraph{Acknowledgements} |
\paragraph{Acknowledgements} |
| 78 |
We thank Jinlun Zhang for providing the original B-grid code and many |
We thank Jinlun Zhang for providing the original B-grid code and many |
| 79 |
helpful discussions. ML thanks Elizabeth Hunke for multiple explanations. |
helpful discussions. ML thanks Elizabeth Hunke for multiple explanations. |
| 80 |
|
|
| 81 |
\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram} |
This work is a contribution to Estimating the Circulation and Climate of the |
| 82 |
|
Ocean, Phase II (ECCO2). The ECCO2 project (http://ecco2.org/) is sponsored |
| 83 |
|
by the NASA Modeling Analysis and Prediction (MAP) program. D. Menemenlis |
| 84 |
|
carried out this work at the Jet Propulsion Laboratory, California Institute |
| 85 |
|
of Technology under contract with the National Aeronautics and Space |
| 86 |
|
Administration. |
| 87 |
|
|
| 88 |
|
\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram,bib/mit_biblio} |
| 89 |
|
|
| 90 |
\end{document} |
\end{document} |
| 91 |
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch