ceaice_split_version/ceaice_part2/ceaice_adjoint.tex

%---------------------------------------------------------------------------
\section{MITgcm adjoint code generation}
\label{sec:adjoint}
%---------------------------------------------------------------------------

There is now a growing body of literature on adjoint applications
in oceanography and adjoint code generation via automatic differentiation (AD).
We therefore limit the description of the method to a brief summary.
For discrete problems as considered here,
the adjoint model operator (ADM) is the transpose of the 
Jacobian or tangent
linear model operator (TLM) of the full (in general nonlinear) forward
model, in this case, the MITgcm coupled ocean and sea ice model. 
Consider a scalar-valued model diagnostics, referred to as
objective function,
and an $m$-dimensional control space (referred to as independent variables)
whose elements we may wish to perturb to assess their impact on the 
objective function.
In the context of data assimilation the objective function may be the
least-square model vs. data misfit, whereas here, we may choose almost
any function that is (at least piece-wise) differentiable with respect to 
the control variables. Here, we shall be focusing on the 
solid freshwater export through Lancaster Sound.
Two- and three-dimensional control variables used in the present
study are listed in Table \ref{tab:controlvars}.

\begin{table}[t!]
\caption{List of control variables used.
The controls are either part of the oceanic (O) or sea-ice (I) state,
or time-varying elements of the atmospheric (A) boundary conditions.}
\label{tab:controlvars}
\begin{tabular}{cccc}
\hline
component & variable & dim. & time  \\
\hline \hline
O & temperature & 3-D & init. \\
O & salinity & 3-D & init. \\
O & vertical diffusivity & 3-D & const. \\
I & concentration & 2-D & init. \\
I & thickness & 2-D & init. \\
A & air temperature & 2-D & 2-day \\
A & specific humidity & 2-D & 2-day \\
A & shortwave radiation & 2-D & 2-day \\
A & precipitation & 2-D & 2-day \\\
A & zonal windspeed  & 2-D & 2-day \\
A & merid. windspeed  & 2-D & 2-day \\
\hline
\end{tabular}
\end{table}

The controls consists, in this study, of two- or
three-dimensional fields of initial conditions of the ocean or sea-ice
state, ocean vertical mixing coefficients, and time-varying
surface boundary conditions (surface air temperature,
specific humidity, shortwave radiation, precipitation, 
zonal and meridional wind speed).
The TLM computes objective functions's directional derivatives 
for a given perturbation direction.
In contrast, the ADM computes the the full gradient
of the objective function with respect to all control variables.
When combined, the control
variables may span a potentially high-dimensional, e.g., O(10$^8$),
control space. At this problem dimension, perturbing
individual parameters to assess model sensitivities is
prohibitive. By contrast, transient sensitivities of the objective
function to any element of the control and model state space can be
computed very efficiently in one single adjoint model integration,
provided an adjoint model is available.

Conventionally, adjoint models are developed ``by hand'' through
implementing code which solves the adjoint equations 
\citep[e.g.][]{marc:95,wuns:96} of the given forward equations.
The burden of developing ``by hand''
an adjoint model in general matches that of
the forward model development. The substantial extra investment
often prevents serious attempts at making available adjoint 
components of sophisticated models. Furthermore, the work of keeping
the adjoint model up-to-date with its forward parent model matches the
work of forward model development.
The alternative route of rigorous application of AD tools has proven
very successful in the context of MITgcm ocean modeling applications.

Certain limitations regarding coding standards apply.
Although they vary from tool to tool, they are similar across various
tools and are related to the ability to efficiently reverse the flow 
through the model. 
Work is thus required initially to make the model amenable to
efficient adjoint code generation for a given AD tool.
This part of the adjoint code generation is not automatic
(we sometimes refer to it as semi-automatic)
and can be substantial for legacy code, in particular if the code
is badly modularized and contains many irreducible control flows
(e.g. GO TO statements, which are considered bad coding practice anyways).

It is important to note, nevertheless, that once the tailoring of the
model code to the AD code is in place, any further forward model 
development can be easily incorporated in the adjoint model via AD.
Furthermore, the notion of \textit{the adjoint} is misleading, since the
structure of the adjoint depends critically on the control problem posed
(a passive tracer sensitivity yields a very different Jacobian
to an active tracer sensitivity). A clear example of the dependence
of the structure of the adjoint model on the control problem
is the extension of an adjoint model which uses bottom topography
as a control variable \citep{losc-heim:07}.
The AD approach enables a much more thorough and smoother 
adjoint model extension than would be possible via hand-coding.
Further comments on advantages and limitations on AD are given
in Appendix  \ref{app:adissues}.

The adjoint model of the MITgcm has become an invaluable
tool for sensitivity analysis as well as for state estimation \citep[for a
recent overview and summary, see][]{heim:08}. 
AD also enables a large variety of configurations 
and studies to be conducted with adjoint methods without the onerous task of
modifying the adjoint of each new configuration by hand.
\cite{gier-kami:98} discuss in detail the advantages of the AD approach.

The AD route was also taken in developing and adapting the sea-ice 
component of the MITgcm, so that tangent linear and adjoint components can be
obtained and kept up to date without excessive effort.
As for the TLM and ADM components of the MITgcm ocean model, we rely on the
automatic differentiation (AD) tool ``Transformation of Algorithms in
Fortran'' (TAF) developed by Fastopt \citep{gier-kami:98} to generate
TLM and ADM code of the MITgcm sea ice model \citep[for details
  see][]{maro-etal:99,heim-etal:05}.
Note that for the ocean component, we are now also able to generate
  efficient derivative code using the new open-source tool OpenAD
  \citep{utke-etal:08}.
Appendix \ref{app:adissues} provides details of
adjoint code generation for the coupled ocean and sea ice MITgcm
configuration.

Since conducting this study, further changes to the
thermodynamic formulation have been implemented, which improve certain
aspects of forward and adjoint model behavior.
These changes are discussed in detail in \cite{fent:09} along with application
of the coupled ocean and sea ice MITgcm adjoint to estimating the state of the
Labrador Sea during 1996--1997.

To conclude this section, we emphasize the coupled nature of the MITgcm 
ocean and sea ice adjoint.
\reffigure{couplingschematic} 
illustrates the relationship between control variables and the
objective function $J$ when using the tangent linear model 
(TLM, left diagram), or the adjoint model (ADM, right diagram).
The left diagram depicts how
each perturbation of an element of the control space
which consists atmospheric perturbations
(surface air temperature $\delta T_a$, precipitation $\delta p$),
sea-ice perturbations
(e.g. ice concentration $\delta c$, ice thickness $\delta h$, 
and oceanic perturbations
(e.g. potential temperature $\delta \Theta$, salinity $\delta S$)
leads to a perturbed objective function $\delta J$).
In contrast, the right diagram shows the reverse propagation of
\textit{adjoint variables} or
\textit{sensitivities} labeled with an asterisk $^{\ast}$.
The notation reflect the fact that adjoint variables are formally
Lagrange multipliers or elements of the model's 
co-tangent space (as opposed to perturbations which are formally
elements of the model's tangent space).
For example, $\delta^{\ast} c$ refers to the gradient
$ \partial J / \partial c$.
The aim of the diagram is to show (in a simplified way) two things:
(1) It depicts how sensitivities of a
sea ice export objective function (to be defined below)
which is purely a function of the sea-ice state,
to changes, e.g., in ice concentration $ \partial J / \partial c$
is affected by changes, e.g. in ocean temperature, i.e.
$ \frac{\partial J}{\partial \theta} \, = \,
  \frac{\partial J}{\partial c} \cdot 
  \frac{\partial c}{\partial \theta} $.
The adjoint model thus maps the adjoint objective function state
to the adjoint sea-ice state, and from there to the coupled
adjoint oceanic and surface atmospheric state.
(2) It can be seen that the ADM maps from a
1-dimensional state ($\delta^{\ast} J$) to a multi-dimensional state
($\delta^{\ast} c, \delta^{\ast} h, \delta^{\ast} T_a, 
\delta^{\ast} p, \delta^{\ast} \Theta, \delta^{\ast} S$), 
whereas the TLM maps from a multi-dimensional state
($\delta c, \delta h, \delta T_a, 
\delta p, \delta \Theta, \delta S$) to a 1-dimensional state
($\delta J$). This is at the heart of the reason, why only
only 1 adjoint integration is needed to assemble the full 
objective function gradient,
but as many tangent linear integrations as dimensions of the control space 
are needed to assemble the same gradient.
Rigorous derivations can be found, e.g., in Chapter 5 of the MITgcm
documentation \citep{adcr-etal:02}, \cite{wuns:06}, or 
\cite{gier-kami:98}.

\begin{figure}[t]
\newcommand{\textinfigure}[1]{{\normalsize\textbf{\textsf{#1}}}} 
\newcommand{\mathinfigure}[1]{\normalsize\ensuremath{{#1}}}
\psfrag{delS}{\mathinfigure{\delta S}}
\psfrag{delT}{\mathinfigure{\delta \Theta}}
\psfrag{delc}{\mathinfigure{\delta c}}
\psfrag{delh}{\mathinfigure{\delta h}}
\psfrag{delAT}{\mathinfigure{\delta T_a}}
\psfrag{delP}{\mathinfigure{\delta p}}
\psfrag{delJ}{\mathinfigure{\delta J}}
%
\psfrag{addS}{\mathinfigure{\delta^{\ast} S}}
\psfrag{addT}{\mathinfigure{\delta^{\ast} \Theta}}
\psfrag{addc}{\mathinfigure{\delta^{\ast} c}}
\psfrag{addh}{\mathinfigure{\delta^{\ast} h}}
\psfrag{addAT}{\mathinfigure{\delta^{\ast} T_a}}
\psfrag{addP}{\mathinfigure{\delta^{\ast} p}}
\psfrag{addJ}{\mathinfigure{\delta^{\ast} J}}
  \centerline{
  \includegraphics*[width=.95\textwidth]{\fpath/coupling_schematic}
  }
  \caption{
    This schematic diagram illustrates how, 
    the tangent linear model (TLM, left panel) maps perturbations in 
    the oceanic, atmospheric, or sea-ice state into a perturbation
    of the objective function $\delta J$, 
    whereas the adjoint model (ADM, right panel) maps the adjoint
    objective function $\delta^{\ast} J$
    into the adjoint sea-ice state
    which is a sensitivity or gradient, e.g., 
    $\delta^{\ast} c \, = \,  \partial J / \partial c$,
    and into the coupled ocean and atmospheric adjoint states.
    The TLM thus computes how a perturbation in \textit{one} input
    affects \textit{all} outputs (but considered here as output is just $J$),
    whereas the adjoint model computes how \textit{one} output 
    (again, just $J$) is affected by \textit{all} inputs.
    \label{fig:couplingschematic}}
\end{figure}

%---------------------------------------------------------------------------
\section{A case study: Sensitivities of sea-ice export through
Lancaster Sound}
%---------------------------------------------------------------------------

We demonstrate the power of the adjoint method in the context of
investigating sea-ice export sensitivities through Lancaster Sound.
The rationale for this choice is to complement the analysis of sea-ice
dynamics in the presence of narrow straits of Part 1. 
Lancaster Sound is one of
the main paths of sea ice export through the Canadian Arctic
Archipelago (CAA)
\citep{mell:02,prin-hami:05,mich-etal:06,muen-etal:06,kwok:06}.
\reffigure{sverdrupbasin} %taken from \cite{mell:02} 
reflects the intricate local geography of CAA
straits, sounds, and islands.
Export sensitivities reflect dominant pathways
through the CAA, as resolved by the model.  Sensitivity maps provide a very
detailed view of
%shed a very detailed light on
various quantities affecting the sea-ice export
(and thus the underlying propagation pathways).  
A caveat of the present study is the limited resolution, which
is not adequate to realistically simulate the CAA.
For example, while the dominant
circulation through Lancaster Sound is toward the East, there is a
small Westward flow to the North, hugging the coast of Devon Island, 
which is not resolved in our simulation.
Nevertheless, the focus here is on elucidating model sensitivities in a 
general way. For any given simulation, whether deemed
``realistic'' or not, the adjoint provides exact model sensitivities, which
help test whether hypothesized processes are actually
borne out by the model dynamics.
Note that the resolution used in this study is at least as good or better 
than the resolution used for IPCC-type calculations.

\begin{figure}[t]
  \centering
  \includegraphics*[width=0.95\textwidth]{\fpath/map_part2_2}
%  \includegraphics*[width=0.9\textwidth]{\fpath/map_sverdrup_basin_melling_2002}
  \caption{Map of the Canadian Arctic Archipelago with model
    coastlines and grid (filled grey boxes are land). The black
    contours are the true coastlines as taken from the GSHHS data base
    \citep{wessel96}. The gate at 82$^{\circ}W$ 
    across which the solid freshwater export is computed 
    is indicated as black line.
    \label{fig:sverdrupbasin}}
\end{figure}

%---------------------------------------------------------------------------
\subsection{The model configuration}
%---------------------------------------------------------------------------

The model domain is similar to the one described in Part 1.
It is carved out from the Arctic face of a global, eddy-admitting,
cubed-sphere simulation \citep{menemenlis05}
but with 36-km instead of 18-km grid cell width, 
i.e., coarsened horizontal resolution compared to 
the configuration described in Part 1.
%, now amounting to roughly 36 km..
The vertical discretization is the same as in Part 1, i.e. the model has
50 vertical depth levels, which are unevenly spaced, ranging from 10 m
layer thicknesses in the top 100 m to a maximum of 456 m layer thickness
at depth.
The adjoint model for this configuration runs efficiently on 80 processors,
inferred from benchmarks on both an SGI Altix and on an IBM SP5 at NASA/ARC
and at NCAR/CSL, respectively.
Following a 4-year spinup (1985 to 1988), the model is integrated for an
additional four
years and nine months between January 1, 1989 and September 30, 1993.
It is forced at the surface 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}.  The air-sea fluxes in the presence of
%sea-ice are handled by the ice model as described in \refsec{model}.
The objective function $J$ is chosen as the ``solid'' freshwater
export through Lancaster Sound (LC), at approximately 74\degN, 82\degW\ in
\reffig{sverdrupbasin}, integrated over the final 12-month period, i.e.,
October 1, 1992 to September 30, 1993.
That is,

\begin{linenomath*}
\begin{equation}
\label{eq:costls}
J \, = 
\frac{1}{\rho_{fresh}}
\int_{\mathrm{Oct\,92}}^{\mathrm{Sep\,93}}
\int_{\mathrm{LC}}
 \, (\rho \, h \, c \, + \, \rho_{s} h_{s}c)\,u \,ds \,dt,
\end{equation} 
\end{linenomath*}

%\ml{[ML: shouldn't $J$ be normalized by $\rho_{\mathrm{fresh}}$ to
%  give the units that we use in the figures?]}
is the mass export of ice and snow converted to units of freshwater.
Furthermore, for each grid cell $(i,j)$ of the section, along which the
integral $\int \ldots ds$ is taken,
$c(i,j)$ is the fractional ice cover, $u(i,j)$ is the along-channel ice drift
velocity, $h(i,j)$ and $h_s(i,j)$ are the ice and snow
thicknesses, and $\rho$, $\rho_s$, and $\rho_{fresh}$
are the ice, snow and freshwater densities, respectively.
At the given resolution, the section amounts to three grid points.
The forward trajectory of the model integration resembles broadly that
of the model in Part~1 but some details are different due
to the different resolution and integration period.
For example, the differences in annual solid
freshwater export through Lancaster Sound as defined in eqn. (\refeq{costls})
are smaller between no-slip and
free-slip lateral boundary conditions at higher resolution,
as shown in Part 1, Section 4.3
($91\pm85\text{\,km$^{3}$\,y$^{-1}$}$ and 
$77\pm110\text{\,km$^{3}$\,y$^{-1}$}$ for free-slip and no-slip, respectively,
and for the C-grid LSR solver; $\pm$ values refer to standard deviations 
of the annual mean) than at lower resolution
($116\pm101\text{\,km$^{3}$\,y$^{-1}$}$ and 
$39\pm64\text{\,km$^{3}$\,y$^{-1}$}$ for free-slip and no-slip, respectively).
The large range of these estimates emphasizes the need to
  better understand the model sensitivities to lateral boundary
  conditions and to different configuration details. We aim to explore
  these sensitivities across the entire model state space in a
  comprehensive manner by means of the adjoint model.
%The large discrepancy between all these numbers underlines the need to
%better understand the model sensitivities across the entire model state space
%resulting from different lateral boundary conditions and different
%configurations, and which we aim to explore in a more
%comprehensive manner through the adjoint.

The adjoint model is the transpose of the tangent linear model
operator.  It runs backwards in time from September 1993 to
January 1989.  During this integration period, the Lagrange multipliers
of the model subject to objective function \refeq{costls} are
accumulated. These Langrange multipliers
are the sensitivities (or derivatives) of the objective function with respect
%ML which can be interpreted as sensitivities of the objective function
to each control variable and to each element of the intermediate 
coupled ocean and sea ice model state variables.
Thus, all sensitivity elements of the model state and of the surface
atmospheric state are
available for analysis of the transient sensitivity behavior.  Over the
open ocean, the adjoint of the \cite{larg-yeag:04} bulk formula scheme computes
sensitivities to the time-varying atmospheric state. 
Specifically, ocean sensitivities propagate to air-sea flux sensitivities,
which are mapped to atmospheric state sensitivities via the
bulk formula adjoint.
Similarly, over ice-covered areas, the sea-ice model adjoint, 
rather than the bulk formula adjoint converts surface ocean sensitivities to
atmospheric sensitivities. 


%---------------------------------------------------------------------------
\subsection{Adjoint sensitivities}
%---------------------------------------------------------------------------

\begin{figure*}[t]
  \includegraphics*[width=\textwidth]{\fpath/adj_canarch_freeslip_ADJheff}
  \caption{Sensitivity $\partial{J}/\partial{(hc)}$ in
    m$^3$\,s$^{-1}$/m for four different times using free-slip
    lateral sea ice boundary conditions. The color scale is chosen
    to illustrate the patterns of the sensitivities.
    The objective function (\refeq{costls}) was evaluated between
    October 1992 and September 1993.
    Sensitivity patterns extend backward in time upstream of the 
    Lancaster Sound section.
    \label{fig:adjhefffreeslip}}
\end{figure*}

\begin{figure*}[t]
  \includegraphics*[width=\textwidth]{\fpath/adj_canarch_noslip_ADJheff}
  \caption{Same as in \reffig{adjhefffreeslip} but for no-slip
  lateral sea ice boundary conditions.
    \label{fig:adjheffnoslip}}
\end{figure*}

The most readily interpretable ice-export sensitivity is that to
ice thickness, $\partial{J} / \partial{(hc)}$.
Maps of transient sensitivities $\partial{J} / \partial{(hc)}$ are shown for
free-slip (\reffig{adjhefffreeslip}) and for no-slip
(\reffig{adjheffnoslip}) boundary conditions.
Each figure depicts four sensitivity snapshots for the objective function $J$, 
starting October 1, 1992, i.e., at the beginning of the 12-month averaging
period, and going back in time to October 2, 1989.
As a reminder, the full period over which the adjoint sensitivities
are calculated is between January 1, 1989 and September 30, 1993.

The sensitivity patterns for ice thickness are predominantly positive.
An increase in ice volume in most places west, i.e., ``upstream'', of
%``upstream'' of
Lancaster Sound increases the solid freshwater export at the exit section.
The transient nature of the sensitivity patterns is evident:
the area upstream of Lancaster Sound that
contributes to the export sensitivity is larger in the earlier snapshot.
In the free slip case, the sensivity follows (backwards in time) the dominant
pathway through Barrow Strait
into Viscount Melville Sound, and from there trough M'Clure Strait
into the Arctic Ocean 
%
\footnote{
(the branch of the ``Northwest Passage'' apparently 
discovered by Robert McClure during his 1850 to 1854 expedition;
McClure lost his vessel in the Viscount Melville Sound)
}. 
%
Secondary paths are northward from
Viscount Melville Sound through Byam Martin Channel into
Prince Gustav Adolf Sea and through Penny Strait into MacLean Strait. 

There are large differences between the free-slip and no-slip
solutions.  By the end of the adjoint integration in January 1989, the
no-slip sensitivities (\reffig{adjheffnoslip}) are generally weaker than the
free slip sensitivities and hardly reach beyond the western end of
Barrow Strait. In contrast, the free-slip sensitivities
(\reffig{adjhefffreeslip})
extend through most of the CAA and into the Arctic interior, both to
the West (M'Clure Strait)  and to the North (Ballantyne Strait, Prince
Gustav Adolf Sea, Massey Sound). In this case the ice can
drift more easily through narrow straits and a positive ice
volume anomaly anywhere upstream in the CAA increases ice export
through Lancaster Sound within the simulated 4-year period.

One peculiar feature in the October 1992 sensitivity maps
are the negative sensivities to the East and, albeit much weaker,
to the West of Lancaster Sound.
The former can be explained by indirect effects: less ice eastward
  of Lancaster Sound results in
less resistance to eastward drift and thus more export.
A similar mechanism might account for the latter,
albeit more speculative: less ice to
the West means that more ice can be moved eastward from Barrow Strait
into Lancaster Sound leading to more ice export. 
%\\ \ml{[ML: This
%  paragraph is very weak, need to think of something else, longer
%  fetch maybe? PH: Not sure what you mean. ML: I cannot remember,
%  either, so maybe we should just leave it as is it, but the paragraph
%  is weak, maybe we can drop it altogether and if reviewer comment on
%  these negative sensitivies we put something back in?]}

\begin{figure*}
  \centerline{
  \includegraphics*[height=.7\textheight]{\fpath/lancaster_adj-line}
  }
  \caption{Time vs. latitude diagrams along the axis of Viscount Melville
    Sound, Barrow Strait, and Lancaster Sound.  The diagrams show the
    sensitivities (derivatives) of the ``solid'' freshwater, i.e.,
    ice and snow, export $J$ through Lancaster Sound
    (\reffig{sverdrupbasin}) with respect to
    ice thickness ($hc$, top panels), ice and ocean surface temperature 
    (in short SST, middle panels), and
    precipitation ($p$, bottom panels) for free slip (left panels) and no slip
    (right panels) sea ice boundary conditions. 
    The objective function was integrated over the last year (period above
    green line). A precipitation perturbation during Oct. 1st 1990 (dashed line)
    or Apr. 1st. 1991 (dash-dottel line) leads to a negative or positive
    ice export anomaly, respectively (see also Fig. \ref{fig:lancpert}).
    Each plot is overlaid with contours 1 and 3 of the normalized ice strength
    $P/P^*=(hc)\,\exp[-C\,(1-c)]$.
    \label{fig:lancasteradj}}
\end{figure*}

The temporal evolution of several ice export sensitivities 
along a zonal axis through
Lancaster Sound, Barrow Strait, and Melville Sound (115\degW\ to
80\degW, averaged across the passages) are depicted in \reffig{lancasteradj}
as Hovmoeller-type diagrams, that is, two-dimensional plots of sensitivities 
as function of longitude and time. 
Serving as examples for the ocean, sea-ice, and atmospheric forcing components
%In order to represent sensitivities to elements of the state of
of the model,  we depict, from top to bottom, the
sensitivities to ice thickness ($hc$), 
ice and ocean surface temperature (in short SST) 
defined as the temperature of the
ocean model component's top grid cell, and precipitation ($p$) for free slip
(left column) and no slip (right column) ice drift boundary conditions.
The green line marks the time from which onward the ice export
objective function was integrated (1 Oct. 1992 to 30 Sep. 1993).
Also indicated are times when a perturbation in precipitation
leads to a negative (Oct. 1990) or to a positive (Apr. 1991) ice export
anomaly (see also Fig. \ref{fig:lancpert}).

The Hovmoeller-type diagrams of ice thickness (top row) and SST
(second row) sensitivities are coherent: 
more ice in Lancaster Sound leads
to more export and one way to form more ice is by colder surface
temperatures. In the free-slip case the
sensitivities spread out in ``pulses'' following a seasonal cycle:
ice can propagate eastward (forward in time) and thus sensitivities
propagate westward (backwards in time) when the ice strength is low
in late summer to early autumn
(\reffig{lancasterfwd1}, bottom panels).
In contrast, during winter, the sensitivities show little to no
westward propagation as the ice is frozen solid and does not move.
In the no-slip case the (normalized)
ice strength does not fall below 1 during the winters of 1991 to 1993
(mainly because the ice concentrations remain near 100\%, not
shown). Ice is therefore blocked and cannot drift eastwards 
(forward in time) through the Viscount
Melville Sound, Barrow Strait, and Lancaster Sound channel system.
Consequently, the sensitivities do not propagate westward (backwards in
time) and the export through Lancaster Sound is only affected by
local ice formation and melting for the entire integration period.

\begin{figure*}
  \centerline{
  \includegraphics*[height=.8\textheight]{\fpath/lancaster_fwd_1-line}
  }
  \caption{Hovmoeller-type diagrams along the axis of Viscount Melville
    Sound, Barrow Strait, and Lancaster Sound.  The diagrams show ice
    thickness ($hc$, top panels), snow thickness ($h_{s}c$, middle panels) and
    normalized ice strength ($P/P^*=(hc)\,\exp[-C\,(1-c)]$, bottom panels) for
    free slip (left panels) and no slip (right panels) sea ice boundary
    conditions.  For orientation, each plot is overlaid with contours 1 and 3
    of the normalized ice strength. 
    Green line is as in Fig. \ref{fig:lancasteradj}.
    \label{fig:lancasterfwd1}}
\end{figure*}

It is worth contrasting the sensitivity 
diagrams of \reffig{lancasteradj}
with the Hovmoeller-type diagrams of the corresponding state variables
(Figs.~\ref{fig:lancasterfwd1} and \ref{fig:lancasterfwd2}).
The sensitivities show clear causal connections of ice motion
over the years, that is, they expose the winter arrest and the summer
evolution of the ice.  These causal connections cannot
easily be inferred from the Hovmoeller-type diagrams of ice and snow
thickness.  This example illustrates the usefulness and complementary nature
of the adjoint variables for investigating dynamical linkages in the
ocean/sea-ice system.

\begin{figure*}
  \centerline{
  \includegraphics*[height=.8\textheight]{\fpath/lancaster_fwd_2-line}
  }
  \caption{Same as in \reffig{lancasterfwd1} but for SST (top panels), SSS
    (middle panels), and precipitation minus evaporation minus runoff, E-P-R
    (bottom panels).
    \label{fig:lancasterfwd2}}
\end{figure*}

The sensitivities to precipitation are more complex.
%exhibit a more complex behaviour.
To first order, they have an oscillatory pattern
with negative sensitivity (more precipitation leads to less export) 
between roughly September and December and mostly positive sensitivity
from January through June (sensitivities are negligible during the summer).
%A fairly accurate description would note an oscillatory behaviour:
%they are negative (more precipitation leads to less export) 
%before January (more precisely, between roughly August and December) 
%and mostly positive after January
%(more precisely, January through July). 
Times of positive sensitivities coincide with times of
normalized ice strengths exceeding values of~3.
This pattern is broken only immediatly preceding the evaluation
period of the ice export objective function in 1992.  In contrast to previous
years, the sensitivity is negative between January and August~1992
and east of 95\degW.
  
We shall elucidate the mechanisms underlying
these precipitation sensitivities
in Section \ref{sec:oscillprecip} 
in the context of forward perturbation experiments.


%---------------------------------------------------------------------------
\subsection{Forward perturbation experiments}
\label{sec:forwardpert}
%---------------------------------------------------------------------------

Applying an automatically generated adjoint model
%Using an adjoint model obtained via automatic differentiation
%and applied 
under potentially highly nonlinear conditions
%, and one generated automatically, relying on AD tools
incites the question 
to what extent the adjoint sensitivities are ``reliable''
in the sense of accurately representing forward model sensitivities.
Adjoint sensitivities that are physically interpretable provide
%Obtaining adjoint fields that are physically interpretable provides
a partial answer but an independent, quantitative test is needed to
gain confidence in the calculations.
%credence to the calculations.
Such a verification can be achieved by comparing adjoint-derived gradients
with ones obtained from finite-difference perturbation experiments.
Specifically, for a control variable $\mathbf{u}$ of interest,
we can readily calculate an expected change $\delta J$ in the objective
function for an applied perturbation $\mathbf{\delta u}$ over domain $A$
based on adjoint sensitivities $\partial J / \partial \mathbf{u}$:

\begin{linenomath*}
\begin{equation}
\delta J \, = \, \int_A \frac{\partial J}{\partial \mathbf{u}} \,
\mathbf{\delta u} \, dA
\label{eqn:adjpert}
\end{equation}
\end{linenomath*}

Alternatively we can infer the magnitude of the objective perturbation $\delta J$
without use of the adjoint.  Instead we apply the same perturbation
$\mathbf{\delta u}$ to the control space over the same domain $A$ and
integrate the forward model.  The perturbed objective function is

\begin{linenomath*}
\begin{equation}
\delta J \, = \, 
J(\mathbf{u}+\mathbf{\delta u}) - J(\mathbf{u}).
\label{eqn:fdpert}
\end{equation}
\end{linenomath*}

The degree to which Eqns.~(\ref{eqn:adjpert}) and (\ref{eqn:fdpert}) agree
depends both on the magnitude of perturbation $\mathbf{\delta u}$
and on the length of the integration period. 
%(note that forward and adjoint models are evaluated over the same period).

We distinguish two types of adjoint-model tests.  First there are finite
difference tests performed over short time intervals,
over which the assumption of linearity is expected to hold,
and where individual elements of the control vector are perturbed.
We refer to these tests as gradient checks.  Gradient checks are performed
on a routine, automated basis for various MITgcm verification setups,
including verification setups that exercise coupled ocean and sea ice model
configurations.  These automated tests insure that updates to the MITgcm
repository do not break the differentiability of the code.

\begin{table*}
\caption{Summary of forward perturbation experiments
and comparison of adjoint-based and finite-difference-based objective function
sensitivities.  All perturbations were applied to a region centered at
101.24$^{\circ}$W, 75.76$^{\circ}$N. The reference value for ice and snow
export through Lancaster Sound is $J_0$ = 69.6 km$^3/yr$.
For perturbations to the time-varying precipitation $p$ the perturbation
interval is indicated by $ \Delta t$. 
}
\label{tab:pertexp}
\centering
\begin{tabular}{ccc@{\hspace{2ex}}c@{\hspace{2ex}}cr@{\hspace{2ex}}r@{\hspace{2ex}}r}
\hline
\textsf{exp.} & variable & time & $\Delta t$ & $\mathbf{\delta u}$ & 
$\frac{\delta J(adj.)}{km^3/yr}$ & $\frac{\delta J(fwd.)}{km^3/yr}$ &
\% diff. \\
\hline \hline
\textsf{ICE1} & $hc$ & 1-Jan-89 & init. & 0.5 m & 0.98 & 1.1 & 11 \\
\textsf{OCE1} & SST & 1-Jan-89 & init. & 0.5$^{\circ}$C & -0.125 & -0.108 & 16 \\
\textsf{ATM1} & $p$ & 1-Apr-91 & 10 dy & 1.6$\cdot10^{-7}$ m/s & 0.185 & 0.191 & 3 \\
\textsf{ATM2} & $p$ & 1-Nov-90 & 10 dy & 1.6$\cdot10^{-7}$ m/s & -0.435 & -1.016 & 57 \\
\textsf{ATM3} & $p$ & 1-Apr-91 & 10 dy & -1.6$\cdot10^{-7}$ m/s & -0.185 & -0.071 & 62 \\
\textsf{ATM4} & $p$ & 1-Nov-90 & 10 dy & -1.6$\cdot10^{-7}$ m/s & 0.435 & 0.259 & 40 \\
\hline
\end{tabular}
\end{table*}

A second type of adjoint-model tests is
finite difference tests performed over longer time intervals
% comparable to the ones used for actual sensitivity studies such as this one, 
and where a whole area is perturbed, guided by the adjoint sensitivity maps,
in order to investigate physical mechanisms.
The examples discussed herein and summarized in Table \ref{tab:pertexp}
are of this second type of sensitivity experiments.
For nonlinear models, the deviations between Eqns.~(\ref{eqn:adjpert}) and
(\ref{eqn:fdpert}) are expected to increase both with
perturbation magnitude as well as with integration time.

\begin{figure}
  %\centerline{
  \subfigure %[$hc$]
  {\includegraphics*[width=.49\textwidth]{\fpath/lanc_pert_heff-box}}

  \subfigure %[SST]
  {\includegraphics*[width=.49\textwidth]{\fpath/lanc_pert_theta-box}}

  \subfigure %[$p$]
  {\includegraphics*[width=.49\textwidth]{\fpath/lanc_pert_precip-box}}
  %}
  \caption{
    Difference in solid freshwater export at 82$^{\circ}$W between perturbed
    and unperturbed forward integrations.  From top to bottom, perturbations
    are initial ice thickness (\textsf{ICE1} in Table \ref{tab:pertexp}),
    initial sea-surface temperature (\textsf{OCE1}), and precipitation
    (\textsf{ATM1} and \textsf{ATM2}). The grey box indicates the period
    during which the ice export objective function $J$ is integrated.
    \label{fig:lancpert}}
\end{figure}

Comparison between finite-difference and adjoint-derived ice-export
perturbations show remarkable agreement for initial value perturbations of
ice thickness (\textsf{ICE1}) or sea surface temperature (\textsf{OCE1}).
Deviations between perturbed objective function values remain below 16\% (see Table
\ref{tab:pertexp}).
\reffigure{lancpert} depicts the temporal evolution of
perturbed minus unperturbed ice export through Lancaster Sound for initial ice
thickness
(top panel) and SST (middle panel) perturbations.
In both cases, differences are confined to the melting season, during which
the ice unlocks and which can lead to significant export.
Large differences are seen during (but are not confined to) the period
during which the ice export objective function $J$ is integrated (grey box).
As ``predicted'' by the adjoint, the two curves are of opposite sign
and scales differ by almost an order of magnitude.

%---------------------------------------------------------------------------
\subsection{Sign change of precipitation sensitivities}
\label{sec:oscillprecip}
%---------------------------------------------------------------------------

Our next goal is to explain the sign and magnitude changes through time 
of the transient precipitation sensitivities.
To investigate this, we have carried out the following two perturbation
experiments: (i) an experiment labeled \textsf{ATM1}, in which we perturb
precipitation over a 10-day period between April 1 and 10, 1991, coincident
with a period of positive adjoint sensitivities, and (ii) an experiment
labeled \textsf{ATM2}, in which we apply the same perturbation over a 10-day
period between November 1 and 10, 1990, coincident with a period of negative
adjoint sensitivities.
The perturbation magnitude chosen is
$\mathbf{\delta u} = 1.6 \times 10^{-7}$ m/s, which is
of comparable magnitude with the standard deviation of precipitation.
%as a measure of spatial mean standard deviation of precipitation
%variability.  The results are as follows: First
The perturbation experiments confirm the sign change
when perturbing in different seasons.
We observe good quantitative agreement for the April 1991 case
and a 50\% deviation for the November 1990 case.
%
The discrepancy between the finite-difference and adjoint-based sensitivity
estimates results from model nonlinearities and from the multi-year
integration period.
To support this statement, we repeated perturbation experiments \textsf{ATM1}
and \textsf{ATM2} but applied a perturbation with opposite sign, i.e.,
$\mathbf{\delta u} = -1.6 \times 10^{-7}$ m/s (experiments \textsf{ATM3} and
\textsf{ATM4} in Table \ref{tab:pertexp}).
For negative $\mathbf{\delta u}$, both perturbation periods lead to about
50\% discrepancies between finite-difference and adjoint-derived
ice export sensitivities.
%
The finite-difference export changes are different in amplitude for positive
and for negative perturbations, confirming that model nonlinearities impact
these calculations.

These experiments constitute severe tests of the adjoint model in the sense
that they push the limit of the linearity assumption.  Nevertheless, the
results confirm that adjoint sensitivities provide useful qualitative, and
within certain limits quantitative,
information of comprehensive model sensitivities that
cannot realistically be computed otherwise.

\begin{figure*}
  \centerline{
  \includegraphics*[width=.9\textwidth]{\fpath/lancaster_pert_hov}
  }
  \caption{
    Same as in \reffig{lancasterfwd1} but restricted to the period
    1991--1993 and for the differences
    in (from top to bottom)
    ice thickness $(hc)$, snow thickness $(h_\mathrm{snow}c)$, sea-surface
    temperature (SST), and shortwave radiation (for completeness)
    between a perturbed and unperturbed run in precipitation of
    $1.6\times10^{-1}\text{\,m\,s$^{-1}$}$ on November 1, 1990 (left panels)
    and on April 1, 1991 (right panels).  For orientation, each plot is
    overlaid with contours 1 and 3 of the normalized ice strength
    $P/P^*=(hc)\,\exp[-C\,(1-c)]$.
    \label{fig:lancasterperthov}}
\end{figure*}

To investigate in more detail the oscillatory behavior of precipitation
sensitivities
we have plotted differences in ice thickness, snow thicknesses, and SST,
between perturbed and unperturbed simulations 
along the Lancaster Sound axis as a function of time.
\reffigure{lancasterperthov} shows how the
small localized perturbations of precipitation are propagated,
depending on whether applied during \textit{early} winter (left column)
or \textit{late} winter (right column).
More precipation
leads to more snow on the ice in all cases. 
However, the same perturbation in different
seasons has an opposite effect on the solid freshwater export
through Lancaster Sound.
Both the adjoint and the perturbation results suggest the following
mechanism to be at play:
%ML: why not let LaTeX do it? Elsevier might have it's own layout
\begin{itemize}
\item
More snow in November (on thin ice) insulates the ice by reducing 
the effective conductivity and thus the heat flux through the ice. 
This insulating effect slows down the cooling of the surface water 
underneath the ice. In summary, more snow early in the winter limits the ice growth 
from above and below (negative sensitivity).
\item 
More snow in April (on thick ice) insulates the
ice against melting.
Shortwave radiation cannot penetrate the snow cover and snow has
a higher albedo than ice (0.85 for dry snow and 0.75 for dry ice in our
simulations); thus it protects the ice against melting in the spring,
more specifically, after January, and it leads to more ice in the 
following growing season.
\end{itemize}
% \\ $\bullet$
% More snow in November (on thin ice) insulates the ice by reducing 
% the effective conductivity and thus the heat flux through the ice. 
% This insulating effect slows down the cooling of the surface water 
% underneath the ice. In summary, more snow early in the winter limits the ice growth 
% from above and below (negative sensitivity).
% \\ $\bullet$
% More snow in April (on thick ice) insulates the
% ice against melting.
% Short wave radiation cannot penetrate the snow cover and has
% a higher albedo than ice (0.85 for dry snow and 0.75 for dry ice in our
% case); thus it protects the ice against melting in spring 
% (more specifically, after January), and leads to more ice in the 
% following growing season.

A secondary effect is the
accumulation of snow, which increases the exported volume.
The feedback from SST appears to be negligible because
there is little connection of anomalies beyond a full seasonal cycle.

We note that the effect of snow vs rain seems to be irrelevant 
in explaining positive vs negative sensitivity patterns.
  In the current implementation, the model differentiates between
  snow and rain depending on the thermodynamic growth rate of sea ice; when it
  is cold enough for ice to grow, all precipitation is assumed to be
  snow. The surface atmospheric conditions most of the year in the Lancaster 
  Sound region are such that almost all precipitation is treated as snow,
  except for a short period in July and August; even then, air
  temperatures are only slightly above freezing. 
  
Finally, the negative sensitivities to precipitation between 95\degW\ and
85\degW\ during the spring of 1992, which break the oscillatory pattern,
may also be explained by the presence of
  snow: in an area of large snow accumulation 
  (almost 50\,cm: see \reffig{lancasterfwd1}, middle panel), 
  ice cannot melt and it
tends to block the channel so that ice coming from the West cannot
pass, thus leading to less ice export in the next season.
%
%\ml{PH: Why is this true for 1992 but not 1991?}
The reason why this is true for the spring of 1992 but not for the spring of
1991 is that by then the high
  sensitivites have propagated westward out of the area of thick 
  snow and ice around 90\degW.

%(*)
%The sensitivity in Baffin Bay are more complex.
%The pattern evolves along the Western boundary, connecting
%the Lancaster Sound Polynya, the Coburg Island Polynya, and the
%North Water Polynya, and reaches into Nares Strait and the Kennedy Channel.
%The sign of sensitivities has an oscillatory character
%[AT FREQUENCY OF SEASONAL CYCLE?].
%First, we need to establish whether forward perturbation runs
%corroborate the oscillatory behaviour.
%Then, several possible explanations:
%(i) connection established through Nares Strait throughflow
%which extends into Western boundary current in Northern Baffin Bay.
%(ii) sea-ice concentration there is seasonal, i.e. partly
%ice-free during the year. Seasonal cycle in sensitivity likely
%connected to ice-free vs ice-covered parts of the year.
%Negative sensitivities can potentially be attributed
%to blocking of Lancaster Sound ice export by Western boundary ice
%in Baffin Bay.
%(iii) Alternatively to (ii), flow reversal in Lancaster Sound is a possibility
%(in reality there's a Northern counter current hugging the coast of
%Devon Island which we probably don't resolve).

%Remote control of Kennedy Channel on Lancaster Sound ice export
%seems a nice test for appropriateness of free-slip vs no-slip BCs.

%\paragraph{Sensitivities to the sea-ice area}

%\refigure{XXX} depicts transient sea-ice export sensitivities
%to changes in sea-ice concentration
% $\partial J / \partial area$ using free-slip
%(left column) and no-slip (right column) boundary conditions.
%Sensitivity snapshots are depicted for (from top to bottom) 
%12, 24, 36, and 48 months prior to May 2003.
%Contrary to the steady patterns seen for thickness sensitivities,
%the ice-concentration sensitivities exhibit a strong seasonal cycle
%in large parts of the domain (but synchronized on large scale).
%The following discussion is w.r.t. free-slip run.

%(*)
%Months, during which sensitivities are negative:
%\\
%0 to 5   Db=N/A, Dr=5 (May-Jan) \\
%10 to 17 Db=7, Dr=5 (Jul-Jan) \\
%22 to 29 Db=7, Dr=5 (Jul-Jan) \\
%34 to 41 Db=7, Dr=5 (Jul-Jan) \\
%46 to 49 D=N/A \\
%%
%These negative sensitivities seem to be connected to months
%during which main parts of the CAA are essentially entirely ice-covered.
%This means that increase in ice concentration during this period
%will likely reduce ice export due to blocking
%[NEED TO EXPLAIN WHY THIS IS NOT THE CASE FOR dJ/dHEFF].
%Only during periods where substantial parts of the CAA are
%ice free (i.e. sea-ice concentration is less than one in larger parts of
%the CAA) will an increase in ice-concentration increase ice export.

%(*)
%Sensitivities peak about 2-3 months before sign reversal, i.e.
%max. negative sensitivities are expected end of July
%[DOUBLE CHECK THIS].

%(*) 
%Peaks/bursts of sensitivities for months
%14-17, 19-21, 27-29, 30-33, 38-40, 42-45

%(*) 
%Spatial ``anti-correlation'' (in sign) between main sensitivity branch
%(essentially Northwest Passage and immediate connecting channels),
%and remote places.
%For example: month 20, 28, 31.5, 40, 43.
%The timings of max. sensitivity extent are similar between
%free-slip and no-slip run; and patterns are similar within CAA,
%but differ in the Arctic Ocean interior.

%(*) 
%Interesting (but real?) patterns in Arctic Ocean interior.

%\paragraph{Sensitivities to the sea-ice velocity}

%(*)
%Patterns of ADJuice at almost any point in time are rather complicated
%(in particular with respect to spatial structure of signs).
%Might warrant perturbation tests.
%Patterns of ADJvice, on the other hand, are more spatially coherent,
%but still hard to interpret (or even counter-intuitive
%in many places).

%(*)
%``Growth in extent of sensitivities'' goes in clear pulses:
%almost no change between months: 0-5, 10-20, 24-32, 36-44
%These essentially correspond to months of


%\subsection{Sensitivities to the oceanic state}

%\paragraph{Sensitivities to theta}

%\textit{Sensitivities at the surface (z = 5 m)}

%(*)
%mabye redo with caxmax=0.02 or even 0.05

%(*)
%Core of negative sensitivities spreading through the CAA as 
%one might expect [TEST]:
%Increase in SST will decrease ice thickness and therefore ice export.

%(*)
%What's maybe unexpected is patterns of positive sensitivities
%at the fringes of the ``core'', e.g. in the Southern channels
%(Bellot St., Peel Sound, M'Clintock Channel), and to the North
%(initially MacLean St., Prince Gustav Adolf Sea, Hazen St.,
%then shifting Northward into the Arctic interior).

%(*)
%Marked sensitivity from the Arctic interior roughly along 60$^{\circ}$W
%propagating into Lincoln Sea, then
%entering Nares Strait and Smith Sound, periodically
%warming or cooling[???] the Lancaster Sound exit.

%\textit{Sensitivities at depth (z = 200 m)}

%(*)
%Negative sensitivities almost everywhere, as might be expected.

%(*)
%Sensitivity patterns between free-slip and no-slip BCs
%are quite similar, except in Lincoln Sea (North of Nares St),
%where the sign is reversed (but pattern remains similar).

%\paragraph{Sensitivities to salt}

%T.B.D.

%\paragraph{Sensitivities to velocity}

%T.B.D.

%\subsection{Sensitivities to the atmospheric state}

%\begin{itemize}
%%
%\item
%plot of ATEMP for 12, 24, 36, 48 months
%%
%\item
%plot of HEFF for 12, 24, 36, 48 months
%%
%\end{itemize}


%\reffigure{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.
%\reffigure{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}}
%%}
%%
%%\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}


%\ml{[based on the movie series
%  zzz\_run\_export\_canarch\_freeslip\_4yr\_1989\_ADJ*:]} The ice
%export through the Canadian Archipelag is highly sensitive to the
%previous state of the ocean-ice system in the Archipelago and the
%Western Arctic. According to the \ml{(adjoint)} senstivities of the
%eastward ice transport through Lancaster Sound (\reffig{sverdrupbasin})
%with respect to ice volume (thickness), ocean
%surface temperature, and vertical diffusivity near the surface
%(\reffig{fouryearadj}) after 4 years of integration the following
%mechanisms can be identified: near the ``observation'' (cross-section
%G), smaller vertical diffusivities lead to lower surface temperatures
%and hence to more ice that is available for export. Further away from
%cross-section G, the sensitivity to vertical diffusivity has the
%opposite sign, but temperature and ice volume sensitivities have the
%same sign as close to the observation.


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice_part2"
%%% End: 