\documentstyle[epsf,psfig,12pt,titlepage]{article}
\topmargin=-.5in
\oddsidemargin=0in
\evensidemargin=0in
\textwidth=6.5in
\textheight=9.0in
\setcounter{totalnumber}{1}

\openup5\jot

\begin{document}

\bibliographystyle{siam}

\title{Dynamical Interpolation of Mesoscale Flows in the Topex/Poseidon Diamond
       Surrounding the U.S. JGOFS Bermuda Atlantic Time-series Site}

\author {D.J. McGillicuddy, Jr. and V.K. Kosnyrev \
           \thanks{Department of Applied Ocean Physics and Engineering,
                   Woods Hole Oceanographic Institution Woods Hole, MA 02543,
                   USA.  Tel: 508-289-2683 Fax: 508-457-2194
                   Email: dmcgillicuddy@whoi.edu (Corresponding Author).}
        }

\date {Draft Manuscript \\
       \today}

\maketitle

\begin{abstract}

\end{abstract}

\noindent
{\bf Introduction}

Ocean time series measurements at fixed locations are among the 
most valuable means for assessing temporal trends in physical,
biological and chemical properties of the water column.  Depth-profile
data at such sites provide a one-dimensional window into oceanic
variability on time scales ranging from the sampling interval to
the length of the record.  However, three-dimensional effects can 
pose some challenges for interpretation.  Without additional
information about the surrounding waters, it can be difficult to 
distinguish temporal changes from spatial heterogeneity moving 
past the site.  Moreover, if such three three-dimensional phenomena
play a role in maintaining the mean conditions, their 
signature in a one-dimensional time series may not be adequate to 
deduce the mechanisms by which this impact is expressed.  In such 
circumstances, it is of interest to obtain the spatial context in 
which the time series observations are embedded.

Motivation to provide such information around the U.S. JGOFS
Bermuda Atlantic Time-series Study (BATS) site arose from the
recognition that mesoscale phenomenology is relevant to both
aspects: mesoscale variability is aliased into the
time series record, and the underlying processes contribute
significantly to biogeochemical cycling in the Sargasso Sea.
Analysis of BATS hydrographic data with a one-dimensional
model revealed month-to-month variability that cannot be explained
by local air-sea fluxes; these anomalies are almost certainly
associated with mesoscale advection (Doney, 1996).  Comparison
of two nutrient profiles sampled one month apart in the summer of
1986 suggested an eddy-driven nutrient injection event that
could account for 20-30\% of the annual new production (Jenkins, 
1988).  Since that time, additional evidence has accumulated
which suggests that mesoscale eddies supply a significant portion
of the annual nutrient budget (McGillicuddy {\it et al.}, 1998).
Regional numerical simulations suggest eddy-induced upwelling
causes intermittent fluxes of nitrate into the euphotic zone
of magnitude sufficient to balance the nutrient demand implied
by geochemical estimates of new production (McGillicuddy and
Robinson, 1997).  Nitrate flux calculations based on satellite
altimetry and a statistical model linking sea level anomaly to
subsurface isopycnal displacements provide estimates of comparable
order (Siegel {\it et al.}, 1999).  Analysis of observations 
capable of resolving mesoscale fluctuations reveals that they
are replete with this phenomenology.  High resolution time
series from the Bermuda Testbed Mooring program (Dickey {\it et al.},
1999) documented a nutrient pulse and associated increase in
chlorophyll and particulate material associated with the passage
of an eddy (McNeil {\it et al.}, 1999).  Mesoscale biogeochemical
surveys carried out as part of the BATS validation activities
demonstrate that eddy-induced upward displacement of density 
surfaces can inject nutrients into the euphotic zone and result 
in the accumulation of biomass in the overlying waters (McGillicuddy 
{\it et al.}, 1999).  Finally, satellite data show covariation in
mesoscale ocean color and sea surface temperature patterns in the 
Sargasso Sea which are consistent with these ideas (McGillicuddy
{\it et al.}, accepted).

\nocite{doney96}
\nocite{jenkins88b}
\nocite{mcgetal98a}
\nocite{mr97}
\nocite{siegeletal99}
\nocite{dickeyetal98}
\nocite{mcneiletal99}
\nocite{mcgetal99a}
\nocite{mcgetal99b}

The technology for detecting oceanic mesoscale phenomena from
satellite-based altimetric measurements dates back to the late
1970s with the GEOS-3 (Tapley {\it et al.}, 1982) and Seasat
(Stanley, 1979) missions.  Since then, regular coverage over a
large portion of the world ocean has been provided by Geosat
[1986-1989] (Douglas and Cheney, 1990), the European Space 
Agency's ERS satellites [1991-present] (AVISO, 1997)
and Topex/Poseidon [1992-present] (Fu {\it et al.}, 1994).  As
with any mapping operation, competing interests of spatial and
temporal resolution determine the ground tracks of the satellite's
orbital configuration.  For example, the Topex/Poseidon 
(hereafter T/P) instrument completes its global coverage once 
every ten days; adjacent ground tracks are spaced approximately 
250km apart at the latitude of Bermuda.  In contrast, ERS-1 and 
ERS-2 tracks are much closer (about 70km around Bermuda) at the 
expense of a longer repeat cycle (35 days).  The T/P and ERS
missions occupy opposite ends of the space/time resolution
requirements needed for direct diagnosis of mesoscale processes
in this region.  Observations from the MODE and POLYMODE programs 
documented eddy sizes in the Sargasso Sea on the order of 100-150km, 
with synoptic dynamical time scales of 10-15 days (Richman {\it et al.}, 
1977; The MODE Group, 1978; Harrison and Heinmiller, 1983).  Thus,
the T/P repeat cycle can adequately resolve temporal variations in
the eddies near Bermuda, but the ground tracks are far enough apart
that an entire eddy can fit in between them.  Closer spacing of the
ERS tracks insures satisfactory spatial coverage, but the repeat
cycle is so long that important temporal dynamics can be missed.
Clearly some form of dynamical interpolation is needed to provide
a basis from which to diagnose mesoscale fluctuations in the
BATS time series.

\nocite{tapleyetal82}
\nocite{stanley79}
\nocite{douglascheney90}
\nocite{aviso97}
\nocite{fuetal94}
\nocite{richmanetal77}
\nocite{mode78}
\nocite{harrisonheinmiller83}

A variety of approaches have been used to incorporate satellite
altimetry into numerical models.  Excellent overviews of the
subject of ocean data assimilation are provided in Bennett (1992),
Wunsch (1996) and Malanotte-Rizzoli (1996); only a few illustrative
examples pertaining to the use of altimetry are mentioned here.
Some of the earliest applications with limited-area models utilized
simple methods of data insertion for model initialization
and boundary condition specification ({\it e.g.} DeMey and Robinson,
1987).  Although such approaches could be used to keep relatively
small regional models on track for short time periods ({\it e.g.}
Dombrowsky and DeMey, 1992), the need for long-term simulations
over large regions necessitated continuous assimilation of data
in the interior.  One of the most straightforward techniques for
doing so is to ``nudge'' the numerical model toward observations
via a forcing term added to the model equations.  Nudging with
altimeter data has been used to construct simulations of ocean
circulation in a variety of regions ({\it e.g.} Holland and
Malanotte-Rizzoli, 1989; Blayo {\it et al.},
1994, 1997; Capotondi {\it et al.}, 1995ab; Stammer, 1997).  
Optimal interpolation is a slightly
more sophisticated approach in which the error characteristics of
the observations are used to determine the weighting of model and
data in the blending procedure. [Example of OI with altimeter data?]
A different class of assimilation techniques involves use of the
model dynamics as a constraint in the estimation problem; this
includes the Kalman filter (and smoother) and the adjoint
method.   Although direct implementation of the
Kalman filter continues to be impractical for large-scale 
oceanographic problems, approximations to this algorithm have
proven to be effective in using altimeter data in both regional 
applications ({\it e.g.} Fukumori and Malanotte-Rizzoli, 1995)
and global general circulation models ({\it e.g.} Fukumori 
{\it et al.}, 1999).  The adjoint method has been used successfully
with altimeter data in several different contexts to invert for 
control variables (such as initial conditions, boundary conditions 
and model parameters) which minimize an objective function measuring
goodness of fit to the observations ({\it e.g.} Moore, 1991;
Schr\"oter {\it et al.}, 1993; Morrow and DeMey, 1995; Vogeler and 
Schr\"oter, 1999).

\nocite{bennett92}
\nocite{wunsch96}
\nocite{rizzoli96}
\nocite{demeyrobinson87}
\nocite{dombrowskydemey92}
\nocite{hollandrizzoli89}
\nocite{blayoetal94}
\nocite{blayoetal97}
\nocite{capotondietal95a}
\nocite{capotondietal95b}
\nocite{stammer97}
\nocite{fukumoririzzoli95}
\nocite{fukumorietal99}
\nocite{moore91}
\nocite{schroteretal93}
\nocite{morrowdemey95}
\nocite{vogelerschroter99}

Nearly all of the examples of altimetric data assimilation described
above deal with model domains that are regional and larger.  That is,
the operational area spans much more than a single set of adjacent
ground tracks.  The intent here is to provide mesoscale context for
time series observations at BATS-- a single point within the 
diamond-shaped domain bounded by two ascending and two descending
T/P lines.  Our approach is to configure a model domain so that its
boundaries are matched to the T/P diamond as closely as possible.
This facilitates specification of the required boundary conditions
more or less directly from altimetry.  For a domain of sufficiently
small dimensions (such as this one), evolution of the interior will
be almost entirely controlled by the boundary conditions.  In a sense,
the model is used as a tool to dynamically interpolate in between
T/P ground tracks.

[Overview]
This study consists of two main elements.  First, simulated data
are used to conduct an identical twin experiment to determine the 
expected error levels in the model interior given perfect boundary
condition information.  Second, the T/P diamond model is forced
with actual altimetric observations from 1992-1999.  Its performance
is evaluated with additional observations in the interior.

\vskip .2in
\noindent
{\bf Methods}

\vskip .2in
\noindent
{\it Model description}

A quasigeostrophic (QG) model is used for the assimilation experiments
described herein.  Relevance of the QG equations to open 
ocean mesoscale flows has been established on theoretical grounds 
(Bretherton and Karweit, 1975; Charney and Flierl, 1981), analysis of 
observations (McWilliams, 1976), and numerical simulation (Haidvogel, 1983).
A detailed description of the model is provided in
in Robinson and Walstad (1987) and references therein, so it will be 
described only briefly here.  The model has been applied in a number
of areas throughout the world ocean, with particularly extensive use 
in the Sargasso Sea.  Of particular interest here are ocean prediction 
experiments which demonstrate predictive skill in the region (Robinson
and Leslie, 1985; Carton, 1987; Walstad and Robinson, 1990).  In addition,
it has been used for detailed dynamical process studies of these flows
(Pinardi and Robinson, 1987; Spall, 1989).  The prognostic equation is

$$ \zeta_t + \alpha J (\psi , \zeta) + \beta \psi_x = F_{pqr} \eqno(1) $$

\noindent
where $\psi$ is the streamfunction and $\zeta$ is the dynamic vorticity
given by

$$ \zeta = \nabla^2_H \psi + \Gamma^2 (\sigma \psi_z)_z. \eqno(2) $$

\noindent
The Jacobian $J$ is defined as
$J(\psi , \zeta) = \psi_x \zeta_y - \psi_y \zeta_x$,
and $F_{pqr}$ represents
a Shapiro filter that is used to parameterize subgridscale dissipation.
The stratification is given by $\sigma(z) = {N_0^2 \over N^2}$ where $N$
is the Brunt-V\"ais\"al\"a frequency and $N_0$ its typical scale.
Computational parameters as well as the nondimensional numbers $\alpha$, 
$\beta$, and $\Gamma^2$ are described in Table 1.  Lateral boundary 
conditions are of the Charney, Fj\"ort\"oft, and von Neumann type: 
streamfunction is prescribed everywhere along the boundary, whereas
vorticity is specified only at inflow points. [surface and bottom
boundary conditions]. 

\nocite{brethertonkarweit75}
\nocite{charneyflierl81}
\nocite{mcwilliams76}
\nocite{haidvogel83}
\nocite{robinsonwalstad87}
\nocite{robinsonleslie85}
\nocite{carton87}
\nocite{walstadrobinson90}
\nocite{pinardirobinson87}
\nocite{spall89}

\vskip .2in
\noindent
{\it Simulated data and identical twin experiments}

In order to assess the error characteristics of the T/P diamond model
if it were fed perfect boundary condition information, identical twin
experiments were conducted with simulated data.  These ``data'' are
generated by forcing the model along the boundaries of a
large (1000 km$^2$) domain with time dependent statistical representations
of typical synoptic features in the region.  This boundary forcing
consists of Rossby wave fits to the MODE
data computed by McWilliams and Flierl (1976), including the surface
intensification of modal amplitudes documented by McWilliams (1974).
Starting from an arbitrary initial condition, the model is run out
to statistical steady state.  Far enough away
from the direct boundary forcing, the freely evolving interior takes on
a realistic character (Figure \ref{driver}).  This technique has been
used to produce environmental fields for forecasting studies in the region
by Robinson and Haidvogel (1980) and Miller and Robinson (1984), as well
as nutrient transport calculations by McGillicuddy and Robinson (1997).
This configuration was deliberately chosen over other alternatives (such
as spin-down in a doubly periodic domain) to facilitate long term simulations
from which robust statistics can be derived.

\nocite{mcwilliamsflierl76}
\nocite{mcwilliams74}
\nocite{robinsonhaidvogel80}
\nocite{millerrobinson84}
\nocite{mr97}

Simulated fields are extracted from a 250 km$^2$ interior subdomain,
which is of approximately the same dimensions as the T/P diamond model.
Results archived on a daily basis constitute a ``perfectly'' sampled 
representation of the ocean which obeys the model equations to within 
numerical accuracy.  This solution can thus be used as a basis to
evaluate the capability of the T/P diamond model to reconstruct evolution
of its interior based only on information provided at the boundaries.

\vskip .2in
\noindent
{\it Altimetric observations and the assimilation procedure}

Merged sea level anomaly (SLA) data from ERS-1, ERS-2 and T/P were
obtained via CD-ROM from the {\it Archiving, Validation and Interpretation 
of Satellite Oceanographic data} (AVISO) project (AVISO, 1997).  ERS and 
T/P data were merged with a global minimization algorithm using dual
crossover differences (Le Traon, {\it et al.}, 1995; Le Traon and
Ogor, 1998) which results in compatible orbit error characteristics.
Ground tracks in the vicinity of Bermuda are shown in Figure \ref{tracks}.
The BATS site lies near the southwestern edge of a T/P diamond, very
near and ERS crossover point.

\nocite{aviso97}
\nocite{letraonetal95}
\nocite{letraonogor98}

Field estimates of SLA were derived from the alongtrack observations
using an objective analysis (OA) algorithm based on that of Carter and
Robinson (1987).  The details of this procedure are provided in Siegel
{\it et al.} (1999) and therefore are presented here only briefly.
A correlation model of the form

$$ C(R,\Delta t) = {\rm exp} \biggl( {R \over c_1} \biggr)
                   {\rm cos} \biggl( 2 \pi {R \over c_2} \biggr)
                   {\rm exp} \biggl( -{\mid \Delta t \mid \over c_3} \biggr) $$

\noindent
was utilized, where the space/time distance
$R = \sqrt{(\Delta x - c_4 \Delta t)^2 + \Delta y^2}$
accounts for the observed westward phase propagation of features in
this region.  This function was fit to the SLA data to obtain
the empirical constants $c_1 = 175$ km, $c_2 = 575$ km, $c_3 = 26$ days and
$c_4 = -4.3$ km day$^{-1}$.  Daily SLA maps were produced on a grid of
approximately 15km spacing, which was subsequently interpolated onto a
10km mesh for incorporation into the model.

\nocite{carterrobinson87}
\nocite{siegeletal99}

Objectively analyzed altimetric observations provide a direct estimate
of the upper level streamfunction in the QG model via the relation

$$ \psi_1 = {g \over V_0 f_0 D} {\rm SLA} $$

\noindent
where $g$ is the acceleration of gravity and the remaining parameters
are defined in Table \ref{modpar}.  Near-surface streamfunction is
projected vertically onto the deeper levels using empirical orthogonal
functions (EOFs) derived from 0-4000m dynamic height anomalies
computed from BATS hydrography.

$$ \psi_k = \psi_1 { {\rm EOF}_k \over {\rm EOF}_1 } \qquad k=2,6 $$

\noindent
In this study only the first EOF is used, as it contains XX\% of
the variance (Figure \ref{eofs}).  Because the CTD data do not
extend all the way to the bottom, the barotropic contribution to
SLA variability is not accounted for in the EOF calculations.
For the purposes of this extension procedure, a constant partitioning
between baroclinic and barotropic modes is assumed, such that the
barotropic component is treated through an addition of an empirically
derived constant to the EOF profile.  This constant is evaluated
through sensitivity analysis (see below).

The full three-dimensional streamfunction field created
by vertical projection of objectively analyzed SLA provides all the
information necessary to compute the two terms of the dynamic vorticity
defined in Equation (2).  Spatial derivatives in these vorticity 
calculations tend to introduce significant noise into the fields,
which can have an adverse effect on the quality of the boundary
conditions.  To counterbalance this effect, the along-boundary
vorticity is smoothed with a two-dimensional (space-time) Gaussian
filter.  Sensitivity experiments indicated that the optimal
filtering parameters for this application are half-width scales
of 50km and 7.5 days.

The procedure described above thus provides a method for obtaining
all the initial and boundary condition information needed to drive
the dynamical model (Equations 2-3) from altimetric observations
alone.  The ``assimilation'' method is simply to initialize the
model with these analyses and provide boundary condition updates
on a daily basis (with linear interpolation in between).  Thus
there are no data assimilated in the interior other than those
used in the initial condition.  Beyond that, only boundary condition 
information is provided.  This has been called XX by Bennett.


\vskip .2in
\noindent
{\bf Results}

\vskip .2in
\noindent
{\it Assimilation of Simulated Data}

A set of three experiments were carried out to ascertain the degree to
which error characteristics of dynamical interpolation in a model
domain the size of the T/P diamond depend on the frequency at which
boundary condition information is supplied.  In each case, the model
is started at rest and then fed boundary conditions from the simulated
data set at the prescribed intervals (with linear interpolation in
between updates). 
Table \ref{simtable} summarizes the results of these computations for
boundary updating of 1, 4 and 10 day intervals.  Temporal evolution
of the streamfunction and vorticity errors is shown in Figure 
\ref{simdata}.  For the time scales associated with mesoscale flows
in this region, one-day boundary condition updating is essentially
continuous.  Ten-day updating reflects that which can be provided by
the T/P altimeter given its repeat cycle, and four-day updating
represents an intermediate case.  In all three experiments, there
is rapid adjustment of the streamfunction field, owing to the
elliptic nature of the equation.  Vorticity adjustment is slightly
more protracted due to the time required for information to
propagate into the domain.  Given that the model fields are fully 
equilibrated on a time scale of less than 100 days, the statistics 
reported in Table \ref{simtable} are computed for days 100-800.

The results of one-day and four-day updating experiments are nearly
identical, with normalized RMS streamfunction errors of less than 10\%.  
Examination of the fields shows that these are mostly composed
of slight errors in phase and amplitude.  Occasional events drive
the streamfunction error higher than average.  For example, around
day 500 a large anticyclone begins to propagate into the domain from
the east, intensifying as it moves westward.  This intensification
occurs mostly in the {\it interior} of the domain, and is largely
driven by interactions with the surrounding eddies which are mostly
{\it outside} the domain.  The geometry of this situation-- that a
highly dynamic event straddles the edge of the model domain-- degrades the
ability of the model to reconstruct the evolution of the interior
from boundary condition information alone.  In this particular instance,
this issue results in underprediction of the anticyclone's intensification,
causing the RMS streamfunction error to approach 20\% (Figure \ref{simpanel1}).
Note, however, that despite the increase in relative error,
the general characteristics of the simulated fields track that of
the reference ocean quite closely.  For animated visualizations
of the results, see {\tt http://science.whoi.edu/users/mcgillic/tpd/sim.html}.

Error characteristics for 10-day updating are noticeably degraded
in comparison with the prior cases.  Although the RMS streamfunction
error is acceptable on average (22\%), there are sporadic
episodes during which error levels climb beyond 50\% 
(Figure \ref{simdata}).  Just as in
the prior cases, these periods of anomalously high error are
associated with eddy-eddy interactions which straddle the model
boundary.  Whereas more frequent boundary updating is able to keep
the interior solution firmly on track, in certain circumstances
10-day updating allows for fairly significant departures from the
reference ocean.  Nevertheless, even for the incidents of highest
streamfunction error, the character of the
simulated fields is similar to that of the reference ocean
(Figure \ref{simpanel10}).

Scatterplots of predicted versus reference SLA for the three
experiments reveal yet another aspect of the error characteristics
(Figure \ref{simscat}).  Updating at 1-day and 4-day intervals
simply produce scatter about the 1:1 line, with slopes of the best
fit line which are not significantly different from 1.0 at the XX\% 
confidence interval.  However, 10-day updating introduces significant
smoothing of features, such that the amplitude of the predicted SLA 
in the interior is only 83\% of that in the reference ocean.

\vskip .2in
\noindent
{\it Assimilation of T/P and ERS Data}

Results in the preceding section were based on simulated data for
which the underlying physics (provided by the model) was perfectly
known and the boundary conditions were sampled without error.  Use of
real observations to simulate the natural system is significantly
more complicated, owing to imperfections in the model, noise in the
data, and errors in the approximations inherent in the assimilation
procedure.  Associated with such a simulation is a collection of
control parameters for which optimal values are not known
{\it a priori}.  

The approach here is to conduct a set of sensitivity
experiments in which these parameters are systematically varied
within the envelope of reasonable limits to ascertain an ``optimal''
set of values using objective metrics of simulation skill.  Clearly,
this method is not gauranteed to result in a truly optimal set unless
all possible combinations of values are tried (which is obviously
impractical).  Nevertheless, by executing a large number of runs (many
hundreds) it is possible to examine a significant portion of parameter
space and thereby assess the robustness of the results with respect
to these unknowns.  Rather than document the entire ensemble of runs
carried out in this undertaking, a base case will be presented
along with selected sensitivity experiments.  First, the central
run will be described which is the most faithful representation
of the natural system that we were able to achieve.  Then a series
of sensitivity experiments about that central case will be shown which
illustrate how the results depend on the particular choices of
control parameters. 

Results of the central simulation are available in animated form
on the web at {\tt http:XX}.
example snapshots at particular times are shown in Figure \ref{central}.
There are a variety of metrics which can be used to evaluate the
skill of these simulations.  The time series of normalized RMS
difference between simulated SLA and objective analysis of the
altimetric observations (averaged over the model domain, exclusive
of boundary points) provides a measure of how well the model tracks 
observed eddy structures in an integrated sense (Figure \ref{eval}).  Note
that this is not a perfectly independent metric; points close the edges of
the model domain will obviously be correlated with those used in the
boundary condition forcing.  The robustness of this metric relies
very heavily on ERS data to fill in the interior of the T/P diamond.
For time periods when ERS data are not available (as in 1994 and
the beginning of 1995), this measure is reduced to an assessment of
the degree to which dynamical interpolation of T/P data by the model 
matches purely statistical interpolation via objective analysis.
[By and large things are good.  Generally there are small phase 
errors, as on day XX (Figure \ref{central}, top row).  Sometimes
there are large errors (Figure \ref{central}, bottom row) which
happen when XX.  All in all the average error of 28% is not
too bad.]

A more independent test is provided by direct comparison of 
predicted SLA with individual ERS observations in the interior of
the model domain (Figure \ref{eval}c).  This requires sampling the model
output in space and time along each of the ERS ground tracks,
extracting the predicted SLA at the model gridpoint which is
closest to each of the ERS observations.  This comparison leads
to three quantitative metrics:  the RMS difference between the
predicted and observed values, their correlation coefficient,
and the slope of the best fit line (the y-intercept of the best
fit line is not included in this list because in most cases it
is not significantly different from zero).  

The BATS CTD data 
provides a completely independent measure of hindcast skill,
as these observations were not used at all in the boundary
condition forcing.  The time series of predicted SLA versus
that inferred from hydrographic profiles at BATS provides a measure
of how well the model is able to reconstruct mesoscale fluctuations
at this particular point.  Presentation of these two quantities in
ensemble format (Figure \ref{eval}d) provides an additional three metrics 
analogous to those introduced for the ERS comparisons shown in 
Figure \ref{eval}c.  Note that the CTD-derived SLA estimates are based on
observations which do not extend all the way to the bottom
and thus do not include barotropic contributions to the total
variance.  For the purposes of this comparison, model SLA estimates
are computed with the same reference level (4000m) used in the dynamic
height calculations with the CTD data.

Sensitivity Analysis

The central run evaluated in Figure \ref{eval} reflects the best
fit to the observations we were able to obtain through variation
of the model control parameters within reasonable limits.  These
control parameters include the deep amplitude of the EOF used in
the streamfunction projection scheme, the space and time scales
on which the vorticity boundary condition is filtered, and the
nondimensional parameters of the dynamical model itself 
(Table \ref{modpar}).  A series of experiments will 
now be presented to quantify the sensitivity of the results to
the particular choices of parameters used in the central run.

The various metrics of hindcast skill are well-behaved over
a wide range of choices for the EOF amplitude at depth 
(Figure \ref{sens_dea}).  For amplitudes less than the
central value of 0.3, RMS differences between predicted and
observed SLA values based on both altimetry and BATS hydrography
rise significantly.  Changes in the model performance are less
pronounced for amplitudes larger than the central value.
Thus, the central value appears to be at the center of a
local minimum in the RMS error measures.  The statistics of
the least squares fit between predicted and observed SLA 
values are relatively insensitive to the amplitude of the
EOF at depth.  The central value produces simulations with
correlation coefficients and slopes of the best fit lines
to the altimetric and hydrographic data which are at or
near their maximal values (although the differences between
runs are minor in these metrics).

Filtering the along-boundary vorticity turns out to be
a crucial step in the hindcasting scheme.  A time series
of the vorticity computed directly from the streamfunction
reveals the noisiness of the field introduced by the
spatial derivatives involved in the calculation (Figure \ref{vort}a).
Smoothing this time series in the space-time phase plane
with a Gaussian filter with half-widths of 50km and 7.5 days
reduces this noise considerably (Figure \ref{vort}b).  The 
necessity of this step is borne out by the serious degradation
of hindcast skill when the unfiltered vorticity is used
(Table) - xx happens. 
A suite of numerical experiments (not shown here) were used
to determine that these filtering parameters were optimal.

ND parameters

\vskip .2in
\noindent
{\bf Discussion}

1. simulated data

this is an error floor-- assumes perfect sampling and that the
ocean is perfectly QG.

\vskip .2in
\noindent
{\bf Conclusions}



\noindent
{\bf Acknowledgments}
We gratefully acknowledge the support of JPL, NASA and NSF.
Collaboration with Prof. David Siegel contributed a
great deal to this effort.  We are very thankful for Rod
Johnson's help with the BATS data; the exceptional efforts of the
BATS technicians are greatly appreciated.

Special thanks to Olga Kosnyreva for her help in data processing
and visualization of model output.

Satellite SLA products used in this study were supplied by the 
CLS Space Oceanography Division, Toulouse, France; 
ERS products were generated as part of the proposal ``Joint analysis of
ERS-1, ERS-2 and TOPEX/POSEIDON altimeter data for oceanic
circulation studies'' selected in response to the Announcement of Opportunity
for ERS-1/2 by the European Space Agency (Proposal code: A02.F105).


% We gratefully acknowledge the support of JPL contract 960517,
% NASA grant NAG5-6455, and NSF grant OCE-9725974.

\vskip .1in
\noindent
Correspondence and requests for materials to D.J. McGillicuddy, Jr.

\bibliography{ocean}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure Captions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\begin{figure}
\caption{Simulated non-dimensional streamfunction field at 100m
after 880 days of integration.  The entire model domain is 1000km on a
side, and the dimensions of the interior subdomain from which simulated
data are extracted are 250km by 250km.  Contour interval is 0.2$\deg$C.}
\label{driver}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Ground tracks for T/P (solid lines) and ERS (dotted lines)
altimetry.  Repeat cycles are 10 and 35 days, respectively.  The bold
line indicates the boundaries of the rectangular model domain fit to
the T/P diamond.}
\label{tracks}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{The first three empirical orthogonal functions (EOFs) of
dynamic height anomalies computed from BATS hydrography; also shown
are mean profiles of density and buoyancy frequency for the region.}
\label{eofs}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Time series of (a) streamfunction and (b) vorticity errors for
three assimilation experiments using simulated data.  Bold, dashed
and thin lines represent 1-day, 4-day and 10-day boundary condition
updating, respectively.  Error statistics represent the RMS difference
between simulated fields and the reference solution integrated over
the model domain.}
\label{simdata}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Top row: Simulated streamfunction and vorticity at 100m on 
day 539 from the experiment with daily boundary condition updates.
Middle row: the same fields from the reference solution.  Bottom row: 
difference between simulated and reference fields.}
\label{simpanel1}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Top row: Simulated streamfunction and vorticity at 100m on 
day 595 from the experiment with ten-day boundary condition updates.
Middle row: the same fields from the reference solution.  Bottom row: 
difference between simulated and reference fields.}
\label{simpanel10}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Predicted versus reference ocean SLA for (a) 1-day, (b) 4-day,
and (c) 10-day updating experiments.  RMS differences between predicted
and reference values are shown in the upper left corner of each plot.
Statistics of the least-squares fit is shown at the lower right:
correlation coefficient ($r$), y-intercept ($a$), and slope ($b$).}
\label{simscat}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Example snapshots of model output from the central run.}
\label{central}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Evaluation of hindcast skill for the central run:
(A) Time series of normalized RMS difference between 
predicted and objectively analyzed
SLA in the interior of the model domain.  Mean value for the entire
simulation is indicated in the upper left hand corner.  Note the
time period in 1994 and the beginning of 1995 for which ERS data are
not available.  
(B) Time series of predicted SLA versus that inferred from CTD observations
at the BATS site.  Note that this comparison does not treat barotropic
contributions to SLA due to the fact that the CTD data do not reach all
the way to the bottom (see text).
(C) Predicted SLA versus that observed from ERS in the interior of the
model domain.  The RMS error is indicated in the upper left corner,
and the statistics of a least squares fit between these two quantitites
are shown on the lower right ($r$ is the correlation coefficient, $a$ is
the y-intercept, and $b$ is the slope).
(D) Predicted SLA versus that inferred from BATS CTD observations.
This is an ensemble presentation of the time-domain comparison shown
in panel (B).  Statistics of the comparison are analogous to those
in panel (C).}
\label{eval}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Sensitivity of the hindcast skill to the EOF amplitude at
depth: (a) normalized RMS difference between predicted SLA and
objective analysis of T/P and ERS altimetric observations; (b) RMS
difference between predicted SLA and observations in the interior
of the model domain (both ERS and BATS hydrography); (c) correlation
coefficient from the least squares fit of the same; (d) slope of the best
fit line to the same.}
\label{sens_dea}
\vspace{.25in}
\end{figure}

\begin{figure}
\caption{Time series of along-boundary vorticity at 100m (a) unfiltered,
and (b) smoothed with a space-time Gaussian filter (see text).}
\label{vort}
\vspace{.25in}
\end{figure}


% \begin{figure}
% \caption{}
% \label{}
% \vspace{.25in}
% \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\begin{table}
\begin{tabular}{|ll|} \hline

\multicolumn{2}{|c|}{\it Computational Parameters} \\

horizontal resolution $\Delta x$         &    15.625 km         \\
vertical resolution                      &    6 levels          \\
time step $\Delta t$                     &    0.0 s             \\

 & \\

\multicolumn{2}{|c|}{\it Physical Scales} \\

time scale $t_0$                         &    0.0 days          \\
space scale $D$                          &    0.0 km            \\
vertical scale $H$                       &    0.0 m             \\
velocity scale $V_0$                     &    0.0 m s$^{-1}$    \\
buoyancy frequency scale $N_0$           &    0.0 s$^{-1}$      \\
Coriolis parameter $f_0$                 &    0.0 s$^{-1}$      \\

 & \\

\multicolumn{2}{|c|}{\it Non-dimensional Parameters} \\

$\alpha = V_0 t_0 D^{-1}$                &    1.5               \\
$\beta = \beta_0 D t_0$                  &    0.5               \\
$\Gamma^2 = f_0^2 D^2 N_0^{-2} H^{-2}$   &    0.3               \\ 


 & \\

\multicolumn{2}{|c|}{\it Additional Control Parameters} \\
\multicolumn{2}{|c|}{\it Used in Hindcasting}           \\

Deep EOF amplitude                       &    0.2               \\
Vorticity filtering                      &    50km, 7.5d        \\ \hline

\end{tabular}
\caption{Model parameters.}
\label{modpar}
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline
Run No. & Updating Freq. (days) & NRMS $\psi$ Error (\%) & 
                                  NRMS $\zeta$ Error (\%) \\ \hline
 1 &  1 &  6.8 & 28.0 \\ \hline
 2 &  4 &  7.8 & 29.4 \\ \hline
 3 & 10 & 21.9 & 40.1 \\ \hline

\end{tabular}
\caption{Assimilation experiments using simulated data in which the
boundary condition updating frequency is varied.  Streamfunction and
vorticity error statistics reported here are the temporal average
(days 100-800) of the RMS differences between simulated fields
and the reference solution integrated over the model domain, normalized
by the time-mean RMS of the reference fields.}
\label{simtable}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/diamd.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{tracks}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/topas1213.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{simdata}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/topas1_n_539.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{simpanel1}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/topas10_n_595.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{simpanel10}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/topas_cor1213.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{simscat}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

% \epsfxsize=16.0cm
% \epsfbox{figs/.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{central}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/dh3a15b5f.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{eval}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\epsfxsize=16.0cm
\epsfbox{figs/sens_dea.eps}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{sens_dea}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

% \epsfxsize=16.0cm
% \epsfbox{figs/}

\vfill
\centerline{McGillicuddy and Kosnyrev, Figure \ref{vort}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\end{document}
