articles/ceaice/ceaice_adjoint.tex

\section{Adjoint sensitivities of the MITsim}
\label{sec:adjoint}

\subsection{The adjoint of MITsim}

The adjoint model of the MITgcm has become an invaluable
tool for sensitivity analysis as well as state estimation \citep[for a
recent summary, see][]{heim:08}. The code has been developed and
tailored to be readily used with automatic differentiation tools for
adjoint code generation. This route was also taken in developing and
adapting the sea-ice compontent MITsim, so that tangent linear and
adjoint components can be obtained and kept up to date without
excessive effort.

The adjoint model operator (ADM) is the transpose of the tangent
linear model operator (TLM) of the full (in general nonlinear) forward
model, in this case the MITsim. This operator computes the gradients
of scalar-valued model diagnostics (so-called cost function or
objective function) with respect to many model inputs (so-called
independent or control variables).  These inputs can be two- or
three-dimensional fields of initial conditions of the ocean or sea-ice
state, 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. At this problem dimension, perturbing
individual parameters to assess model sensitivities quickly becomes
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.

In anology to the TLM and ADM components of the MITgcm we rely on the
autmomatic differentiation (AD) tool ``Transformation of Algorithms in
Fortran'' (TAF) developed by Fastopt \citep{gier-kami:98} to generate
TLM and ADM code of the MITsim \citep[for details see][]{maro-etal:99,
  heim-etal:05}.  In short, the AD tool uses the nonlinear parent
model code to generate derivative code for the specified control space
and objective function. Advantages of this approach have been pointed
out, for example by \cite{gier-kami:98}.

Many issues of generating efficient exact adjoint sea-ice code are
similar to those for the ocean model's adjoint.  Linearizing the model
around the exact nonlinear model trajectory is a crucial aspect in the
presence of different regimes (e.g., is the thermodynamic growth term
for sea-ice evaluated near or far away from the freezing point of the
ocean surface?). Adapting the (parent) model code to support the AD
tool in providing exact and efficient adjoint code represents the main
work load initially. For legacy code, this task may become
substantial, but it is fairly straightforward when writing new code
with an AD tool in mind. Once this initial task is completed,
generating the adjoint code of a new model configuration takes 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
the 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 doing so is to complement the analysis of sea-ice
dynamics in the presence of narrow straits.  Lancaster Sound is one of
the main paths of sea-ice flowing through the Canadian Arctic
Archipelago (CAA).  Export sensitivities reflect dominant pathways
through the CAA as resolved by the model.  Sensitivity maps can shed a
very detailed light on various quantities affecting the sea-ice export
(and thus the underlying pathways).  Note that 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
\citep{mell:02, mich-etal:06,muen-etal:06}, which is not resolved in
our simulation.

The model domain is the same as the one described in \refsec{forward},
but with halved horizontal resolution.
The adjoint models run efficiently on 80 processors (as validated
by benchmarks on both an SGI Altix and an IBM SP5 at NASA/ARC).
Following a 4-year spinup (1985 to 1988), the model is integrated for four
years and nine months between January 1989 and September 1993.
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}.  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'' fresh water
export, that is the export of ice and snow converted to units of fresh
water,
%
\begin{equation}
J \, = \, (\rho_{i} h_{i}c + \rho_{s} h_{s}c)\,u
\end{equation} 
%
through Lancaster Sound at approximately 82\degW\ (cross-section G in
\reffig{arctic_topog}) averaged \ml{PH: Maybe integrated quantity is
more physical; ML: what did you actually compute? I did not scale
anything, yet. Please insert what is actually done.} over the final
12-month of the integration between October 1992 and September 1993.

The forward trajectory of the model integration resembles broadly that
of the model in \refsec{forward}. Many details are different, owning
to different resolution and integration period; for example, the solid
fresh water transport through Lancaster Sound is
%
\ml{PH: Martin, where did you get these numbers from?} 
\ml{[ML: I computed hu = -sum((SIheff+SIhsnow)*SIuice*area)/sum(area) at
$i=100,j=116:122$, and then took mean(hu) and std(hu). What are your numbers?]}
%
$116\pm101\text{\,km$^{3}$\,y$^{-1}$}$ for a free slip simulation with
the C-LSOR solver, but only $39\pm64\text{\,km$^{3}$\,y$^{-1}$}$ for a
no slip simulation. \ml{[Here we can say that the export through
  Lancaster Sound is highly uncertain, making is a perfect candidate
  for sensitivity, bla bla?]}

The adjoint model is the transpose of the tangent linear (or Jacobian) model
operator. It runs backwards in time, from September 1993 to
January 1989. During its integration it accumulates the Lagrange multipliers
of the model subject to the objective function (solid freshwater export),
which can be interpreted as sensitivities of the objective function
to each control variable and each element of the intermediate 
coupled model state variables.
Thus, all sensitivity elements of the coupled
ocean/sea-ice model state as well as the surface atmospheric state are
available for analysis of the transient sensitivity behavior.  Over the
open ocean, the adjoint of the bulk formula scheme computes
sensitivities to the time-varying atmospheric state.  Over ice-covered
areas, the sea-ice adjoint converts surface ocean sensitivities to
atmospheric sensitivities.

DISCUSS FORWARD STATE, INCLUDING SOME NUMBERS ON SEA-ICE EXPORT

\subsubsection{Adjoint sensitivities}

The most readily interpretable ice-export sensitivity is that to
effective ice thickness, $\partial{J} / \partial{(hc)}$.
\reffig{adjheff} shows transient $\partial{J} / \partial{(hc)}$ using
free-slip (left column) and no-slip (right column) boundary
conditions. Sensitivity snapshots are depicted for beginning of October 1992,
that is 12 months before September 1993 
(the beginning of the averaging period for the objective
function $J$, top),
and for Jannuary 1989, the beginning of the forward integration (bottom).
\begin{figure*}[t]
  \includegraphics*[width=\textwidth]{\fpath/adjheff}
  \caption{Sensitivity $\partial{J}/\partial{(hc)}$ in
    m$^2$\,s$^{-1}$/m for two different times (rows) and two different
    boundary conditions for sea ice drift. The color scale is chosen
    to illustrate the patterns of the sensitivities; the maximum and
    minimum values are given above the figures.
    \label{fig:adjheff}}
\end{figure*}

The sensitivity patterns for effective ice thickness are predominantly positive.
An increase in ice volume in most places ``upstream'' of
Lancaster sound increases the solid fresh water export at the exit section.
The transient nature of the sensitivity patterns 
(top panels vs. bottom panels) is also obvious:
the area upstream of the 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 the Barrow Strait
into the Viscount Melville Sound, and from there trough the M'Clure Strait
into the Arctic Ocean (the ``Northwest Passage''). \ml{[Is that really
  the Northwest Passage? I thought it would turn south in Barrow
  Strait, but I am easily convinced because it makes a nicer story.]}
Secondary paths are northward from the
Viscount Melville Sound through the Byam Martin Channel into
the Prince Gustav Adolf Sea and through the Penny Strait into the
MacLean Strait. \ml{[Patrick, all these names, if mentioned in the
  text need to be included somewhere in a figure (i.e. fig1). Can you
  either do this in fig1 (based on martins\_figs.m) or send me a map
  where these names are visible so I can do this unambiguously. I
  don't know where  Byam
  Martin Channel, Prince Gustav Adolf Sea, Penny Strait, MacLean
  Strait, Ballantyne St., Massey Sound are.]}

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

One peculiar feature in the October 1992 sensitivity maps (top panels)
are the negative sensivities to the East and to the West of the
Lancaster Sound.
These can be explained by indirect effects: less ice to the East means
less resistance to eastward drift and thus more export; similarly, less ice to
the West means that more ice can be moved eastwards from the Barrow Strait
into the Lancaster Sound leading to more ice export. 
\ml{PH: The first explanation (East) I buy, the second (West) I
  don't.} \ml{[ML: unfortunately, I don't have anything better to
  offer, do you? Keep in mind that these sensitivites are very small
  and only show up, because of the colorscale. In Fig6, they are
  hardly visible.]}

The temporal evolution of several ice export sensitivities (eqn. XX,
\ml{[which equation do you mean?]}) along a zonal axis through
Lancaster Sound, Barrow Strait, and Melville Sound (115\degW\ to
80\degW, averaged across the passages) are depicted as Hovmueller
diagrams in \reffig{lancasteradj}. These are, from top to bottom, the
sensitivities with respect to effective ice thickness ($hc$), ocean
surface temperature ($SST$) and precipitation ($p$) for free slip
(left column) and no slip (right column) ice drift boundary
conditions.
%
\begin{figure*}
  \includegraphics*[height=.8\textheight]{\fpath/lancaster_adj}
  \caption{Hovermoeller diagrams along the axis Viscount Melville
    Sound/Barrow Strait/Lancaster Sound. The diagrams show the
    sensitivities (derivatives) of the ``solid'' fresh water (i.e.,
    ice and snow) export $J$ through Lancaster sound
    (\reffig{arctic_topog}, cross-section G) with respect to effective
    ice thickness ($hc$), ocean surface temperature (SST) and
    precipitation ($p$) for two runs with free slip and no slip
    boundary conditions for the sea ice drift. Each plot is overlaid
    with the contours 1 and 3 of the normalized ice strengh
    $P/P^*=(hc)\,\exp[-C\,(1-c)]$ for orientation.
    \label{fig:lancasteradj}}
\end{figure*}
%
\begin{figure*}
  \includegraphics*[height=.8\textheight]{\fpath/lancaster_fwd}
  \caption{Hovermoeller diagrams along the axis Viscount Melville
    Sound/Barrow Strait/Lancaster Sound of effective ice thickness
    ($hc$), effective snow thickness ($h_{s}c$) and normalized ice
    strengh $P/P^*=(hc)\,\exp[-C\,(1-c)]$ for two runs with free slip
    and no slip boundary conditions for the sea ice drift. Each plot
    is overlaid with the contours 1 and 3 of the normalized ice
    strength for orientation.
    \label{fig:lancasterfwd}}
\end{figure*}
%

The Hovmoeller diagrams of ice thickness (top row) and sea surface temperature
(second row) sensitivities are coherent: 
more ice in the Lancaster Sound leads
to more export, and one way to get more ice is by colder surface
temperatures (less melting from below). In the free slip case the
sensitivities spread out in "pulses" following a seasonal cycle:
ice can propagate eastwards (forward in time and thus sensitivites can
propagate westwards (backwards in time) when the ice strength is low
in late summer to early autumn.  
In contrast, during winter, the sensitivities show little to now
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, Lancaster Sound channel system.
Consequently, the sensitivities do not propagate westwards (backwards in
time) and the export through Lancaster Sound is only affected by
local ice formation and melting for the entire integration period.

The sensitivities to precipitation exhibit an oscillatory behaviour:
they are negative (more precipitation leads to less export) 
before January (more precisely, late fall) 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
%
\ml{PH: Problem is, that's not true for the first two years (backward),
east of 95\degW, that is, in the Lancaster Sound.
For example, at 90\degW\ the sensitivities are negative throughout 1992,
and no clear correlation to ice strength is apparent there.}
except between 95\degW\ and 85\degW, which is an area of
increased snow cover in spring. \ml{[ML: and no, I cannot explain
  that. Can you?]}

%
Assuming that most precipation is snow in this area\footnote{
In the
current implementation the model differentiates between snow and rain
depending on the thermodynamic growth rate; when it is cold enough for
ice to grow, all precipitation is assumed to be snow.}
%
the sensitivities can be interpreted in terms of the model physics.
The accumulation of snow directly increases the exported volume.
Further, short wave radiation cannot penetrate the snow cover and has
a higer 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 (after
January).

On the other hand, snow reduces 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 and limits the ice growth from
below, so that less snow in the ice-growing season leads to more new
ice and thus more ice export.
\ml{PH: Should probably discuss the effect of snow vs. rain.
To me it seems that the "rain" effect doesn't really play a role
because the neg. sensitivities are too late in the fall,
probably mostly falling as snow.} \ml{[ML: correct, I looked at
NCEP/CORE air temperatures, and they are hardly above freezing in
Jul/Aug, but otherwise below freezing, that why I can assume that most
precip is snow. ]} \ml{[this is not very good but do you have anything
better?:]}
The negative sensitivities to precipitation between 95\degW\ and
85\degW\ in spring 1992 may be explained by a similar mechanism: in an
area of thick snow (almost 50\,cm), ice cannot melt and tends to block
the channel so that ice coming in from the West cannot pass thus
leading to less ice export in the next season.

\subsubsection{Forward sensitivities}

\ml{[Here we need for integrations to show that the adjoint
  sensitivites are not just academic. I suggest to perturb HEFF
  and THETA initial conditions, and PRECIP somewhere in the Melville
  Sound and then produce plots similar to reffig{lancasteradj}. For
  PRECIP it would be great to have two perturbation experiments, one
  where ADJprecip is posivite and one where ADJprecip is negative]}
  

%(*)
%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}

%Fig. XXX depcits 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}


%\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}}
%%}
%%
%%\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{arctic_topog},
%cross-section G) with respect to ice volume (effective 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"
%%% End: 