% -*-LaTeX-*-
\chapter[Neutral vectors in a 3-layer QG model]{Neutral vectors and optimal forcing patterns in a 3-layer
quasigeostrophic model}
\label{c:3layer}

In chapter \ref{c:gm99}, we identified a mechanism of coupled atmosphere-ocean
interaction which leads to interannual variability.  In chapter \ref{c:elab},
found that the shape of the coupled mode is determined by nearly-resonant
``neutral vector'' patterns in the atmosphere.  Neutral vector patterns provide
the maximal atmospheric response to external forcing, thus providing an
amplification mechanism.  In this chapter, we take what we have learned about
the role of neutral vectors in coupled air-sea interaction, and attempt to
apply it to a more realistic model of atmospheric dynamics -- namely, the
3-layer quasigeostrophic model used by Marshall and Molteni (1993).

The model and neutral vector computation techniques are described in section
\ref{s:3l_model}.  We discuss the properties of the model's EOFs, neutral
vectors and optimal forcing patterns in section \ref{s:3l_solns}.  Some
mathematical insight into the connection between neutral vectors and EOFs is
given in section \ref{s:3l_nv_eof}.  We attempt to discover whether the optimal
forcing patterns remain optimal in the full nonlinear model in section
\ref{s:3l_nonlin}.

\section{Molteni's 3-layer quasi-geostrophic model}
\label{s:3l_model}

Molteni's model is discussed in some detail in Molteni's thesis (1994).  The
model is a 3-layer, global, spectral model at T21 resolution, with pressure as
a vertical coordinate.  The quasi-geostrophic potential vorticity equations are
discretized at the three pressure levels (200,500,800 mb), giving prognostic
equations for PV:
\begin{eqnarray*}
\deldel{t} q_1= -J(\psi_1,q_1)&- D_1(\psi_1,\psi_2) & + S_1 \\
\deldel{t} q_2= -J(\psi_2,q_2)&- D_2(\psi_1,\psi_2,\psi_3) & + S_2 \\
\deldel{t} q_3= -J(\psi_3,q_3)&- D_3(\psi_1,\psi_2)& + S_3
\label{eq:3l_molteqn}
\end{eqnarray*}
The $\psi_n$ and $q_n$ are the streamfunction and QGPV at each level, $D_n$
encapsulate various dissipative processes (see Molteni, 1994), and the $S_n$
are a constant PV source term. PV is defined as: 
\begin{eqnarray*}
q_1 = \nabla^2 \psi_1 &- R_1^{-2} (\psi_1 - \psi_2) &+ f \\
q_2 = \nabla^2 \psi_2 &+ R_1^{-2} (\psi_1 - \psi_2) - R_2^{-2}(\psi_2-\psi_3) &+ f \\
q_3 = \nabla^2 \psi_3 &+ R_2^{-2} (\psi_2 - \psi_3) &+ f(1+ h/H_0) \\
\end{eqnarray*}

where $R_1$ and $R_2$ are (spatially constant) Rossby radii of deformation, $h$
is the height of topography, and $H_0$ is a topographic scaling factor.

The model covers the entire globe, so it does not have artifical ``walls'' at
the equator or elsewhere which may spuriously reflect planetary waves.
However, since QG dynamics is not really appropriate near the equator (and the
constant values for the Rossby radii are certainly wrong there), the model's
behavior in the tropics should not be taken too seriously.

In the mid-latitudes, however, the model's dynamics are reasonable, and it can
be made to produce a very good mean flow field through careful specification of
the constant forcing fields $S_n$.  This is done by setting the $S_n$ equal to
the opposite of the average PV tendencies obtained by inserting observed
streamfunction fields into a version of the model equations from which the
$S_n$ are omitted.  This forces the model to have a stable climatology which is
near the observed fields used to generate the $S_n$.  The technique is similar
to the ``flux correction'' used to eliminate climate drift in coupled GCMs.
See Molteni (1994) for more details.  (***are more details in M\&M?)  The $S_n$
used for this study are generated using (find dataset source, months, 
years), so the model attempts to simulate the northern wintertime climatology.

\begin{figure}
\epsfysize=8in
\epsffile{figures/model_climate_summary.eps}
\label{f:3l_climate}
\caption[3-layer QG model mean climatology]{3-layer QG model mean climatology,
based on 5000 daily fields.  Top panel: mean streamfunction, 200 mb, contour
interval $20\EE{6}$ m$^2$/s.  Middle panel: mean streamfunction, 500 mb,
contour interval $10\EE{6}$ m$^2$/s.  Lower panel: Model streamfunction
standard deviation, 500 mb, contour interval $2\EE{6}$ m$^2$/s }
\end{figure}

Figure \ref{f:3l_climate} compares various model mean fields with observations;
these fields are computed from a 5000-day integration of the model.
Comparisons with observed wintertime mean streamfunction (not shown) show that
the model does a very good job at reproducing the basic state.  Model eddy
activity, as shown by streamfunction standard deviation, is less accurate, but
the model does have the right amplitude of eddy activity; we do see storm
tracks over the northern hemisphere oceans, at approximately the right
longitudes.  However, the Atlantic storm track in particular does not have a
sharp northern boundary, and Eddy activity is somewhat stronger than observed
over the pole.

(*note*: Should I include a figure showing observed fields to provide a
comparison, or should I assume people know the wintertime climatology of the
atm?)

We use this 5000-day run to compute empirical orthogonal functions of the model
output.  The EOFs shown in section \ref{s:3l_eofs} are computed over the entire
model domain, using daily data.  EOFs computed over weekly- or monthly-average
data are quite similar in shape.

\subsection{Computing neutral vectors} 

We compute neutral vectors using the same general technique as introduced in
section \ref{s:neutral}.  The only difficulty is to construct the $\ms{M}$
matrix.

This job is made much simpler by the existence of a linearized version of the
Molteni model code, written by (**whom?**).  We linearize about the mean state
of the 5000-day run shown in figure \ref{f:3l_climate}.  Using this linearized
code, we simply compute the tendencies of a complete set of orthonormal
spectral Green's function perturbations, and use those tendencies to build up
an explicit matrix $\ms{M}$ column-by-column.  We then use Matlab's eigensolver
to compute the smallest eigenvectors of $\ms{M}^{\dag}\ms{M}$.  We do not use
the Arnoldi technique, because the tendency matrix is dense (but smaller) when
expressed using a spectral basis.  However, the eigenvalues can still be found
fairly rapidly.

One factor remains to be determined: the definition of the inner products in
equation \ref{eq:MM_minimize}.  This boils down to the question: ``We want to
find the patterns whose tendency is smallest... but what do we mean by
smallest?''  For example, we could minimize the global integral of
streamfunction tendency, or the global integral of PV tendency, or an integral
with nonuniform weights in Euclidean or spectral space.  We will find that some
of our results are sensitive to this choice. Table \ref{t:ips} shows the
various inner products for which we have computed neutral vectors.  All of
these inner products were computed over a global domain, with equal weights for
the 3 layers.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|}\hline
Name & Minimization problem \\ \hline \hline
Psi  & $\lambda^2 =
\frac{\ip{\deldel{t}\mb{\Psi}}{\deldel{t}\mb{\Psi}}}
{\ip{\mb{\Psi}}{\mb{\Psi}}}$, \\
& $\ip{a}{b} = \int \int a b$ \\ \hline
Q & $\lambda^2 =
\frac{\ip{\deldel{t}\mb{q}}{\deldel{t}\mb{q}}}
{\ip{\mb{\Psi}}{\mb{\Psi}}}$, \\
& $\ip{a}{b} = \int \int a b$  \\ \hline
KE & $\lambda^2 =
\frac{\ip{\deldel{t}\mb{\Psi}}{\deldel{t}\mb{\Psi}}}
{\ip{\mb{\Psi}}{\mb{\Psi}}}$, \\
& $\ip{a}{b} = \int \int a \nabla^2 b$ \\ \hline
Mix & $\lambda^2 =
\frac{\ip{\deldel{t}\mb{\Psi}}{\deldel{t}\mb{\Psi}}_1}
{\ip{\mb{\Psi}}{\mb{\Psi}}_2}$, \\
& ${\ip{a}{b}}_1 = \int \int a b$; ${\ip{a}{b}}_2 = \int \int a \nabla^2 b$ \\ \hline
\end{tabular}
\end{center}
\caption{Definitions of inner products used in the study}
\label{t:ips}
\end{table}


\section{EOFs, Neutral vectors, and optimal forcing patterns}
\label{s:3l_solns}

\subsection{EOFs}
\label{s:3l_eofs}
The first few EOFs of the model are shown in figure \ref{f:3l_eofs}.  The EOFs
are different from those described by Molteni \etal (1998); this is probably
because Molteni \etal computed hemispheric EOFs of the \emph{eddy} fields
(i.e., zonal-mean components are removed), while we compute global EOFs of the
full streamfunction.  Also, Molteni's patterns are EOFs of observed fields,
while we compute EOFs of the model output -- and the model is, of course, not
perfect.

\begin{figure}
\label{f:3l_eofs}
\epsffile{figures/eofs_500mb_pol.eps}
\caption{Empirical orthogonal functions computed from daily samples of a
5000-day integration of the Molteni QG model.  }
\end{figure}

(**All the citations in the next few paragraphs can be found in that paper
which compares the NAO to the AO**).

Nonetheless, the model EOFs do resemble observations.  The first EOF has the
dipolar nature of the NAO in the Atlantic, but is much more zonally
extensive. It more closely resembles an ``annular mode'', or the ``Arctic
Oscillation'' (AO) (**cites**), which are commonly seen in full-hemisphere EOF
analyses of observations (**cites**) and models (**cites**).  

(**cite**) make a strong case that the NAO and the AO are really the same
phenomenon: an ``annular mode'' which is somewhat stronger in the Atlantic than
elsewhere.  We have noted that papers which use EOFs localized to the Atlantic
domain (Cayan, **date**) tend to find a more localized NAO pattern, while
hemispheric EOF calculations (**cite**) tend to find the annular AO pattern.
In appendix \ref{c:restrict}, we show why this occurs.

The first EOF explains only 12\% of the variance of the daily data, but this
fraction increases as we look at longer-timescale variability.  A
nearly-identical pattern (not shown) explains 18\% of the weekly-averaged, and
25\% of the monthly-averaged data.  The fraction of monthly variance
explained is is in good agreement with EOFs computed for monthly observations
(Cayan, **date**).  It is not immediately obvious why the first few EOFs should
increasingly dominate the variance as we look at longer timescales; we will
address this issue using ``neutral vector thinking'' in section
\ref{s:3l_nv_eof}.

The second EOF of the model is essentially identical to the planetary wave
pattern associated with the PNA (Pacific North-America) pattern (Wallace,
1995).

Like the first two modes, the third EOF also displays a dipolar pattern in both
the Pacific and Atlantic, but where EOF 2 has a ``checkerboard'' pattern of
highs and lows in the two basins
$\left(\begin{array}{cc}H&L\\L&H\end{array}\right)$, EOF 3 has a more zonal
structure $\left(\begin{array}{cc}L&L\\H&H\end{array}\right)$ which is shifted
south of EOF 1, and lies in quadrature with it.

Like EOF 1, the fractional variance explained by EOFs 2 and 3 also increases
for weekly and monthly-averaged data, but not as strongly.  For EOFs beyond
\#7, the fractional variance explained \emph{decreases} when averaging is
performed. The first dozen EOFs have the vast majority of their power in the
northern hemisphere.  Since southern-hemisphere variability is very weak, this
is not too surprising.

\subsection{Neutral Vectors \& Optimal Forcing Patterns}
\label{s:3l_nvs}

When we computed neutral vectors for the 3-layer model in the standard way
***what standard way?*** (not shown), we discovered that while the EOFs were
confined to the northern hemisphere, the neutral vectors were found in either
or both hemispheres.  The reason for this difference will be discussed in
section \ref{s:3l_nv_eof}.  We feel it is unlikely that a dynamical connection
exists between wave patterns in alternate hemispheres; it is more likely that
two separate neutral patterns in the two hemispheres share similar eigenvalues.
The SVD analysis cannot distinguish two modes with similar eigenvalues, and
will return two orthogonal linear superpositions of the two modes.  A ``rotated
neutral vector analysis'', along the lines of rotated EOFS (**cite**), might
help to separate the modes.

Instead of rotating the modes, we focus on the northern hemisphere by adding an
artificial damping term to the $\ms{M}$ matrix.  This term is proportional to
$\sin(\phi/2 - \pi/4)^6$ (where $\phi$ is latitude), and damps PV anomalies
with a timescale of 5 days at the south pole, 8 days at 45\degr S, 40 days
at the equator, and 1500 days at 45\degr N.  Thus, any otherwise-neutral mode
in the southern hemisphere will have a significant tendency due to this damping
effect, making it less neutral.  The general effect is to re-order the neutral
mode patterns, giving preference to northern-hemisphere modes.

In figures \ref{f:3l_kenv_1} through \ref{f:3l_kenv_3}, we show neutral vectors
and optimal forcing patterns computed using the ``kinetic energy'' norm
described by Molteni.  We describe the neutral vectors (right-hand column) in
detail in the following section; the optimal forcing patterns are discussed in
section \ref{s:3l_nvs_ofp}

\subsubsection{Neutral Vectors}
\label{s:3l_nvs_nvs}


\begin{figure}
\epsfysize=8in
\epsffile{figures/ke_nv1_sh_pol_bw.ps}
\caption{Right column: Neutral vector \#1 for model climatology, KE inner
product.  Left column: corresponding optimal forcing pattern.  Contour interval
is arbitrary, but consistent from level to level.}
\label{f:3l_kenv_1}
\end{figure}

\begin{figure}
\epsfysize=8in
\epsffile{figures/ke_nv2_sh_pol_bw.ps}
\caption{Right column: Neutral vector \#2 for model climatology, KE inner
product.  Left column: corresponding optimal forcing pattern.}
\label{f:3l_kenv_2}
\end{figure}

\begin{figure}
\epsfysize=8in
\epsffile{figures/ke_nv3_sh_pol_bw.ps}
\caption{Right column: Neutral vector \#3 for model climatology, KE inner
product.  Left column: corresponding optimal forcing pattern.}
\label{f:3l_kenv_3}
\end{figure}

(**Should I use a Miller rectangular projection for these plots (and EOFs))

The first neutral vector (figure \ref{f:3l_kenv_1}, right column) shows a
roughly zonally-symmetric pattern, with a negative center over the pole
surrounded by a positive annulus in high middle latitudes.  This annulus has
enhanced energy over western Siberia and the North Atlantic.  The whole pattern
is equivalent barotropic, and broadly resembles the first model EOF, although
it is shifted farther to the north.

The second neutral vector resembles the first EOF quite closely; we see broad
positive centers over the northern mid-latitude oceans, forming a
nearly-complete annulus about the globe, with a polar negative center.  Once
again, the mode is barotropic, with amplitude increasing with height.  The
projection of this pattern onto the NAO pattern is quite strong.

The third neutral vector displays a pattern nearly identical to the PNA pattern
(see *cite*).  We see a barotropic wavetrain extending from the subtropical
Pacific over the Gulf of Alaska, northern Canada, and into the subtropical
east Atlantic.  The subtropical Pacific maximum is shifted west across the
dateline compared to the observed PNA, but elsewhere, the resemblance between
neutral vector 3, EOF 2, and the PNA is very strong.

It should be quite obvious that these patterns strongly resembles the model
EOFs and observed patterns of low-frequency variability.  To emphasize this
point, and to demonstrate that these resemblances are not the result of chance
correlations in a low-dimensional system, we present figure \ref{f:3l_starmap}.
This figure shows the cross-correlations between model EOFs (figure
\ref{f:3l_eofs}) and the neutral vectors plotted in a KE norm (horizontal
axis).  The absolute value of the correlation is indicated by the size of the
dots.  Correlations less than 0.2 are not plotted.

\begin{figure}
\epsffile[112 225 499 616]{figures/KENV_EOF_xcor_starmap.eps}
\caption{Cross-correlations between model EOFs (vertical axis) and neutral
vectors computed using a KE norm (horizontal axis).  The absolute value of the
correlation is indicated by the size of the dots.  Correlations less than 0.2
are not plotted.}
\label{f:3l_starmap}
\end{figure}

The observations discussed above are borne out in this figure: neutral vectors
1 and 2 project onto EOF 1, and neutral vector 3 projects onto EOF 2.  We also
see that and neutral vector 1 projects onto EOF 3.  But more importantly, the
strongest correlations are among the first few EOFs and neutral vectors (the
large dots are clustered in the upper left corner).  If the correlations arose
by chance, we would expect this figure to show a random scattering of small
points throughout the domain.

Since the first few neutral vectors project strongly onto the EOFs, we should
expect that a substantial amount of the model's natural variability resides in
the subspace of the first few neutral vectors.  Figure \ref{f:eof_nv_varexpl}
demonstrates that this is the case.  The 3 most neutral vectors explain more of
the variance of a 5000-day model run than any other singular vector.  The first
three EOFs explain 20.8\% of the daily variance, while the first three neutral
vectors explain 14.2\% of the variance -- almost as much.\footnote{Since EOFs,
by definition, maximize explained variance, the variance explained by neutral
vectors must be smaller.}  When considering monthly-averaged data (not shown),
the results are even more striking: the first three EOFs (still computed from
daily data) explain 46\% of the monthly variance, compared to 37\% for the
first 3 neutral vectors.

\begin{figure}
\epsffile{figures/var_exp_eof_kenv.eps}
\caption{Comparison of fractional variance of daily model output explained by
EOFs (dots) with variance explained by neutral vectors (x's).}
\label{f:eof_nv_varexpl}
\end{figure}

Many papers (*cite*) have been written which project a model onto the first few
EOFs to generate a reduced-subspace model which encapsulates most of the
system's variability.  We find here that neutral vectors are equally good for
this purpose.  They also have the advantage over EOFs that they represent
dynamically-important modes of the system, rather than being empirically
selected.

\subsubsection{Optimal Forcing Patterns}
\label{s:3l_nvs_ofp}

The pattern which maximally excites neutral vector 1 shows two broad regions of
sensitivity of opposite signs: one focused on Kamchatka and extending over the
east Pacific and Siberia, and a second focused on the tropical Pacific and
covering Africa, Europe, the North Atlantic, and North America.  The pattern is
baroclinic, with an out-of-phase response in the lower two layers and nearly
zero sensitivity at the upper level.  This implies a sensitivity to low-level
heating -- which hints that this mode may be sensitive to forcing by SST
anomalies.

The second optimal forcing pattern shows a sensitivity to low-level PV forcing
over much of the North Pacific, particularly the Sea of Okhotsk, and to an
opposite sign of forcing over the North Atlantic and the Arctic.  Once again,
the mode is sensitive to low-level heating.

The optimal forcing pattern for neutral vector 3 also shows a global-scale
dipole, with sensitivity to low-level baroclinic forcing. The forcing centers
lie at the beginning and end of the PNA-like wavetrain.

If we considered only these optimal forcing patterns based on the kinetic
energy norm, we would conclude that the optimal forcing patterns are most
sensitive to broad-scale, low-level thermal forcing.  In some cases (like
figure \ref{f:3l_kenv_3}, in which forcing centers lie at the beginning and end
of the wavetrain), the spatial relationship of forcing to response makes sense,
while in others (like figure \ref{f:3l_kenv_1}, in which an east-west dipole
pattern gives rise to a zonally symmetric structure), it's not so obvious.

However, we find that while the neutral vectors are relatively insensitive to
the choice of inner product, the optimal forcing patterns look very different
for different inner products.  This is demonstrated in figure \ref{f:3l_psinv_1},
which shows the first neutral vector and optimal forcing pattern using the psi
norm (see table \ref{t:ips}).  

\begin{figure}
\epsfysize=8in
\epsffile{figures/psi_nv1_sh_pol_bw.ps}
\caption{Right column: Neutral vector \#1 for model climatology, Psi inner
product.  Left column: corresponding optimal forcing pattern.  Compare with
figure \ref{f:3l_kenv_1}}
\label{f:3l_psinv_1}
\end{figure}

The first neutral vector for the psi norm looks virtually identical to the
first neutral vector of the KE norm.  However, the optimal forcing pattern
looks radically different.  Rather than being characterized by a
planetary-sized dipole structure, we see many narrow, closely-spaced zones of
positive and negative sensitivity.  These tend to be strongest in the tropics,
where, as we have already mentioned, the model's dynamics are the least
believable.  The tendency for the optimal forcing to be low-level heating
remains, but is much less prominent.  With the exception of a strong
sensitivity to mid-layer forcing over eastern Siberia, the first optimal
forcing patterns of the KE and psi norms look completely different.

(*Should I move the following discussion  to a later section?*)

Why do the two norms display identical neutral vectors if their optimal forcing
patterns so different?  Consider the atmospheric response equation
$\ms{M}\mb{\Psi} = \mb{f}$.  A neutral vector has $\mb{f}$ small, so the
neutral vector algorithm selects the components of $\mb{\Psi}$ to ensure the
near-cancellation of the various terms in the $\ms{M}$ matrix.  This
cancellation is independent of the norm selected, so the neutral vector pattern
is not norm-sensitive.  However, relatively large changes in $\mb{f}$ can
result from relatively small changes in $\mb{\Psi}$, since the left side of the
equation is a difference of small terms.  If the norm penalizes one sort of
pattern more heavily, the amplitude of that pattern in $\mb{f}$ can be made
small with only small changes to $\mb{\Psi}$.

The psi norm applies equal weight to all wavenumbers, while the KE norm
penalizes high wavenumbers (for which $\nabla \psi$ is large) more heavily.  As
a result, the KE optimal forcing patterns are more broad-scale.  However, this
analysis cannot tell us whether one norm is ``better'' than another.

(*is this paragraph incoherent?*) While I was in the process of performing
these computations, the potential for norm-sensitivity was pointed out by Grant
Branstator (1999).  However, the norm-sensitivity we found was not at all
guaranteed.  A system in which less cancellation of terms within
$\ms{M}\mb{\Psi}$ was possible (for example, a heavily-damped system) would
have less perfect neutral vectors, and would also be less norm-sensitive.  A
trade-off exists between the perfection of neutral vectors and the sensitivity
of optimal forcing patterns.  Unfortunately, it turns out that this system does
not lie within a useful middle ground where both neutral vectors and optimal
forcing patterns are well-defined.  This is indicated by the value of
eigenvalue $\lambda_1$ corresponding to the most neutral pattern.  It
corresponds to a timescale of about two years -- incredibly long compared to
the dominant timescales of the terms in the PV equation (a week or two),
implying almost perfect cancellation of terms.  This maximizes the importance
of the neutral vectors in explaining model variability, but minimizes the
utility of the optimal forcing patterns.

This phenomenon can also be explained in terms of the condition number of the
$\ms{M}$ matrix (defined as the ratio of the largest and smallest signular
values).  When this number is large, the matrix is ``poorly conditioned'', and
the response $\psi$ is sensitive to small perturbations in $\mb{f}$ or
$\ms{M}$.  A singular matrix has an infinite condition number.  For our
$\ms{M}$, the condition number is of the order $10^4$, implying rather large
sensitivity to forcing, but we must stress that the matrix is \emph{not} close
to being singular to within machine precision.

Navarra (1993) followed an essentially identical neutral vector analysis on a
barotropic model, using the January 300-mb climatological flow.  His neutral
vector patterns look completely different from those found in this analysis,
and are much less similar to observed patterns like the NAO, PNA, and AO.
However, Navarro's work agrees that the neutral vectors strongly resemble the
model's EOFs, and that the condition number of the $\ms{M}$ matrix is
relatively large, leading to large sensitivity to forcing.  Interestingly,
Navarro's optimal forcing patterns show some similarity to those shown in
figure \ref{f:3l_psinv_1}, particularly in their fine-scale, zonally-oriented
bands of sensitivity in the tropics.

\section{Relationship between neutral vectors and EOFs}

The most important observation of the previous section is the close connection
between EOFs and neutral vector patterns.  In this section, we present a
mathematical explanation for this connection.\footnote{While writing, I was
informed that most of the results presented in this section have been published
by Navarro (1993).  Since many readers will be unfamiliar with that work, and
since my discussion of the importance of the $\mb{f}\mb{f}^{\dag}$ matrix goes
beyond Navarro's study, I present the derivation in full here.}

Suppose the atmosphere consists of a number of ``slow modes'' (stationary
planetary wave patterns like the NAO and PNA), which respond in a linear way to
stochastic forcing generated by ``fast modes'' (transient eddy forcing), such
that during any observation period $n$, the stochastic forcing $\mb{f}_n$
excites a planetary wave response $\mb{\Psi}_n$:
\[
\ms{M}\mb{\Psi}_n = \mb{f}_n \rightarrow \mb{\Psi}_n =  \ms{M}^{-1}\mb{f}_n
\]
or, defining matrices $\mb{\Psi}$ and $\ms{f}$ whose columns are the
$\mb{\Psi}_n$ and  $\mb{f}_n$:
\[
\mb{\Psi} =  \ms{M}^{-1}\ms{f}
\]
The EOFs are defined as the eigenvalues of the covariance matrix
$\mb{\Psi}^{\dag}\mb{\Psi}$:
\[
\mb{\Psi}\mb{\Psi}^{\dag} =  \ms{M}^{-1}\ms{f}\ms{f}^{\dag} \ms{M}^{-\dag}
\]
Consider for a moment the case where the eddy forcing is isotropic and
spatially uncorrelated, and the equations are normalized so that $\mb{f}$ has
unit amplitude: $\ms{f}\ms{f}^{\dag} = \ms{I}$.
Then
\begin{equation}
\mb{\Psi}\mb{\Psi}^{\dag} =  \ms{M}^{-1}\ms{M}^{-\dag}
\label{eq:3l_covmat}
\end{equation}
We may decompose the $\ms{M}$ matrix using its singular vectors thus:
\[
\ms{M} = \ms{U}\mb{\Lambda}\ms{V}^{\dag}
\]
Inserting this into (\ref{eq:3l_covmat}), and using the fact that $\ms{V}^{-1}
= \ms{V}^{\dag}$ and $\ms{U}^{-1}= \ms{U}^{\dag}$, we obtain:
\[
\mb{\Psi}\mb{\Psi}^{\dag} =
\ms{V}^{-\dag}\mb{\Lambda}^{-1}\ms{U}^{-1}\ms{U}^{-\dag}\mb{\Lambda}^{-1}\ms{V}^{-1}
= \ms{V}\mb{\Lambda}^{-2}\ms{V}^{\dag}
\]
Thus, the EOFs (which are the eigenvectors of the covariance matrix $
\mb{\Psi}\mb{\Psi}^{\dag}$) are identical to the singular vectors of the
tendency matrix.  Since $\mb{\Lambda}$ is raised to the -2 power, the dominant
EOFs correspond to the smallest singular vectors -- that is, to the neutral
vectors.

What happens when the transient eddy forcing covariance is \emph{not}
proportional to the identity matrix?  Then we have:
\begin{equation}
\mb{\Psi}\mb{\Psi}^{\dag} =
\left(\ms{M}^{\dag}\left(\ms{f}\ms{f}^{\dag}\right)^{-1} \ms{M}\right)^{-1}
\label{eq:3l_covmat_f}
\end{equation}
The matrix $\left(\ms{f}\ms{f}^{\dag}\right)^{-1}$ is a symmetric, positive
definite matrix (assuming the inverse exists).  It thus has the structure and
position of a weight matrix for an inner product between the $\ms{M}$ matrix
with itself.  Thus, in the presence of nonuniform stochastic forcing, we can
compare the EOFs to the singular vectors of $\ms{M}$ computed using this
unusual weight matrix $\left(\ms{f}\ms{f}^{\dag}\right)^{-1}$.  We demonstrate
this by performing an SVD decomposition of $\ms{M}$:
\[
\ms{M} = \ms{U}\mb{\Lambda}\ms{V}^{\dag}
\]
where now the orthormality of the $\ms{V}$ and $\ms{U}$ are defined using the
inner products
\[
\ms{V}^{\dag}\ms{V} = \ms{I} \qquad
\ms{U}^{\dag}\left(\ms{f}\ms{f}^{\dag}\right)^{-1}\ms{U} = \ms{I} 
\]
Then we find that
\[
\ms{M}^{\dag}\left(\ms{f}\ms{f}^{\dag}\right)^{-1}\ms{M} =
\ms{V}\mb{\Lambda}\ms{U}^{\dag}\left(\ms{f}\ms{f}^{\dag}\right)^{-1}\ms{U}\mb{\Lambda}\ms{V}^{\dag}
= \ms{V}\mb{\Lambda}^2\ms{V}^{\dag}
\]
And so, returning to (\ref{eq:3l_covmat_f}):
\[
\mb{\Psi}\mb{\Psi}^{\dag} =
\left(\ms{M}^{\dag}\left(\ms{f}\ms{f}^{\dag}\right)^{-1} \ms{M}\right)^{-1}
 = \ms{V}\mb{\Lambda}^{-2}\ms{V}
\]
Thus, the principal EOFs of this system are the neutral
vectors($\ms{V}$), where the SVD analysis used to compute the neutral vectors
uses the inverse of the stochastic forcing correlation matrix to normalize the
optimal forcing patterns ($\ms{U}$).  

What is the significance of this odd weight matrix
$\left(\ms{f}\ms{f}^{\dag}\right)^{-1}$?  Consider the simple case where each
vector element represents a location in space, and the stochastic forcing is
spatially uncorrelated but inhomogeneous.  Then the eddy forcing covariance
matrix will be diagonal, with larger elements on the diagonal where forcing is
strong.  Thus the weight $\left(\ms{f}\ms{f}^{\dag}\right)^{-1}$ is small where
the eddy forcing is large.  In computing neutral vectors, we want to minimize
the forcing needed to excite them.  A forcing pattern will be ``small'' with
respect to this norm when it has small amplitude where the weight is large, and
vice versa.  Thus, the SVD analysis selects neutral vectors whose forcing
patterns are localized at the site of large eddy forcing.

This allows us to explain the observation made in section \ref{s:3l_nvs} that
the EOFs all lie within the northern hemisphere, while the neutral vectors
reside in both hemispheres, unless we force them into the north using an
artifical hemispheric damping.  High-frequency eddy activity, and thus eddy
forcing, are far stronger in the northern hemisphere (where it's wintertime).
Thus, patterns sensitive to northern-hemisphere forcing will be driven more
strongly, and so northern-hemisphere modes will be more active, even though
they're no more ``neutral'' than southern-hemisphere modes.

To observe a closer connection between EOFs and neutral vectors, we should
compute neutral vectors using a norm weighted with the inverse eddy
forcing covariance, rather than any of the norms described in table
\ref{t:ips}.  It is not easy to build this forcing covariance matrix from
observations alone, but since eddy flux is proportional to eddy strength, using
the covariance of high-frequency eddy streamfunction might give good results. 

(*this is an obvious experiment to do.  Does it work?  Do I have time to do
it?)

\section{Response of the nonlinear model to ``optimal'' forcing}
\label{s:3l_nonlin}

We have found in section \ref{s:3l_solns} that the optimal forcing patterns are
difficult to define unambiguously, since they vary dramatically between
different definitions of the inner product.  Nevertheless, it is useful to find
out whether these patterns, which are optimal in forcing the linearized
stationary planetary wave model, are also optimal in forcing the full
non-linear, time-evolving model.  If this model were quasi-linear, we would
expect the results from section \ref{s:3l_solns} to hold equally well.
However, nonlinear effects may spoil this correspondence.  If the optimal
forcing patterns found using the linear model to not excite a strong response
in the nonlinear model, we certainly cannot expect them to tell us much about
the sensitivity of the true atmosphere to PV forcing.

We want to find the time-mean perturbation response to a constant ``optimal''
forcing perturbation.  We proceed by running three integrations of Molteni's
model.  First, we create a pair of control runs, where $S_n$ in
\label{eq:3l_molteqn} are unchanged from the specification described in
\ref{s:3l_model}.  These two runs are initialized with very slightly different
initial conditions.  The difference between the mean state of these two runs
will give us some idea of the uncertainty of the mean, with which we can
compare the experimental run.

In the experimental run, we perturb $S_n$ by a small amount in the direction of
the first ***