\magnification=1040
%\magnification=1150
%\magnification=1200
\input /home/pierrel/texfile/macros/BOLD.tex
\input /home/pierrel/texfile/macros/eplain.tex
%\input /home/pierrel/texfile/macros/psfig.tex
\input /home/pierrel/texfile/macros/psfig.sty
%\input /home/pierrel/texfile/macros/newmacros.tex
\input /home/pierrel/texfile/macros/font.tex
\def\tab_src{/home/pierrel/scotia4/thesis_sav/work/chap5/fig/jaot/ren}
%\def\tab_src{/usr/pacific/pierrel/texfile/PhD/thesis/work/chap8/tab}
%\def\tab_src{/home/scotia4/pierrel/thesis_sav/work/chap5/fig/jaot/ren}
%\input newmacros.tex
%\input BOLD.tex
%
%\hsize=6.25truein
%\vsize=8.0truein
%\hoffset=0.2in
%\voffset=-0.1in
%
%\font\bb=msym10
%\def\bR{{\hbox{\bb R}}}
\def\bR{{\hbox{ R}}}
\def\o{{\rm o}}
\def\a{{\rm a}}
\def\b{{\rm b}}
\def\c{{\rm c}}
\def\t{{\rm t}}
\def\cmsecm1{cm$\,$s$^{-1}$}
\def\cm2secm1{cm$^2\,$s$^{-1}$}
\def\mm1secm1{m$^{-1}$$\,$s$^{-1}$}
\def\secm1{s$^{-1}$}
\def\tps{{\Bpsi}^{\t}}
\def\hps{\widehat{\Bpsi}}
\def\hpss{\widehat{\phi}}
\def\xa{\Bx^{\a}}
\def\xb{\Bx^{\b}}
\def\phib{\phi^{\b}}
\def\phii{{\phi_{\scriptstyle i}}}
\def\phiw{{\phi_{\scriptstyle w}}}
\def\varphii{{\varphi_{\scriptstyle i}}}
\def\varphij{{\varphi_{\scriptstyle j}}}
\def\varphiw{{\varphi_{\scriptstyle w}}}
\def\varphib{\varphi^{\b}}
\def\xbj{\Bx^{\b^{\scriptstyle j}}}
\def\xbone{\Bx^{\b^{\scriptstyle 1}}}
\def\xbq{\Bx^{\b^{\scriptstyle q}}}
\def\EbT{\BE^{\b^{\scriptstyle T}}}
\def\EaT{\BE^{\a^{\scriptstyle T}}}
\def\EcT{\BE^{c^{\scriptstyle T}}}
\def\EcsT{\BE^{c^{{\scriptstyle *}^{\scriptstyle T}}}}
\def\VbT{\BV^{\b^{\scriptstyle T}}}
\def\HpT{\BH^{p^{\scriptstyle T}}}
\def\Pb{\BP^\b}
\def\Pbp{\BP^{\b^{\scriptstyle p}}}
\def\Bap{\BB^{\a^{\scriptstyle p}}}
\def\xtidbj{\widetilde \Bx^{\b^{\scriptstyle j}}}
\def\tpsk{{\Bpsi}^\t_k}
\def\tpskl{{\Bpsi}^\t_{k-1}}
\def\tpskpl{{\Bpsi}^\t_{k+1}}
\def\hps{{\widehat{\Bpsi}}}
\def\hpsk{{\widehat{\Bpsi}_k}}
\def\hpskl{{\widehat{\Bpsi}_{k-1}}}
\def\hpskpl{{\widehat{\Bpsi}_{k+1}}}
\def\hdelk{{\widehat{\Bdelta}_k}}
\def\tdelk{{\Bdelta^\t_k}}
\def\tdelkl{{\Bdelta^\t_{k-1}}}
\def\hdelkl{{\widehat{\Bdelta}_{k-1}}}
\def\hetk{{\widehat{\Beta}_k}}
\def\tetk{{\Beta^\t_k}}
\def\tetkl{{\Beta^\t_{k-1}}}
\def\hetkl{{\widehat{\Beta}_{k-1}}}
\def\Up{\BUpsilon\,}
\def\trc{{\rm trc}}
\tolerance=10000 \hfuzz=5pt

%DEFINING THE CAPITOL LETTERS OF THE ALPHABET
\def\A{{\bf A}}
\def\B{{\bf B}}
\def\C{{\bf C}}
\def\D{{\bf D}}
\def\E{{\bf E}}
\def\F{{\bf F}}
\def\G{{\bf G}}
\def\H{{\bf H}}
\def\I{{\bf I}}
\def\J{{\bf J}}
\def\K{{\bf K}}
\def\L{{\bf L}}
\def\M{{\bf M}}
\def\N{{\bf N}}
\def\O{{\bf O}}
\def\P{{\bf P}}
\def\Q{{\bf Q}}
\def\R{{\bf R}}
\def\S{{\bf S}}
\def\T{{\bf T}}
\def\U{{\bf U}}
\def\V{{\bf V}}
\def\W{{\bf W}}
\def\X{{\bf X}}
\def\Y{{\bf Y}}
\def\Z{{\bf Z}}

\def\bPhi{{\BPhi}}
\def\bPsi{{\BPsi}}
\def\bpsi{{\Bpsi}}
\def\bDelta{{\BDelta}}
\def\bLambda{{\BLambda}}
\def\bUpsilon{{\BUpsilon}}
\def\bOmega{{\BOmega}}
%\def\bSigma{\hbox{\bgreek\char006}}
\def\bSigma{\BSigma}
\def\bPi{{\BPi}}
\def\bGamma{{\BGamma}}
\def\subphi{{\Phi}}
\def\bnabla{{\bf \nabla}}
\def\btau{{\Btau}}
\def\cen_rule{\nointerlineskip \centerline{\vbox{\hrule width6truein}}
\nointerlineskip}


%\baselineskip=18pt plus 1pt minus 1pt
\baselineskip=12pt plus 1pt minus 1pt
\parskip=5pt
\font\bigbf= cmbx10 scaled \magstep 1
\font\biggbf= cmbx10 scaled \magstep 2
\font\bigrm= cmr10 scaled \magstep 1
\font\smrm= cmr8
\font\bigsl= cmsl10 scaled \magstep 1
\font\caps=cmssbx10 at 12truept
\font\capsf=cmss10 at 9truept
\font\capsfb=cmssbx10 at 9truept
\font\capsl=cmsl10 at 9truept
\font\capsfn=cmss10 at 10truept
\font\scaps=cmssbx10 at 10truept
%\font\scaps=cmssbx10 at 10.6truept
\font\sscaps=cmssbx10 at 10truept

\vfill\eject
\pageno=1
\centerline {\bigbf ON THE MAPPING OF}
\vskip 0.3truecm
\centerline {\bigbf MULTIVARIATE GEOPHYSICAL FIELDS:}
\vskip 0.3truecm
\centerline {\bigbf STUDIES OF}
\vskip 0.3truecm
\centerline {\bigbf THE SENSITIVITY TO ERROR SUBSPACE PARAMETERS}
\vskip 0.5truecm
\centerline {\bigbf P.F.J.~Lermusiaux, C.J.~Lozano and D.G.M.~Anderson}
\centerline {\bigbf Harvard University}
\vskip 1truecm

\centerline {\bf Abstract}
\vskip 0.5truecm

{\tenpoint
\baselineskip=12pt plus 1pt minus 1pt

The influence of {\it a priori} parameters
on the error and variability subspace estimation and
corresponding three-dimensional, multivariate and multiscale
mapping methodology described in
(Lermusiaux {\it et al.}~(2000) is investigated.
The sensitivities of the subspace and of the
{\it a posteriori} fields and errors
to the size of the subspace, scales considered
and nonlinearities in the dynamical adjustments
are exemplified and studied.
The applications focus on the 
mesoscale to subbasin-scale physical fields 
in the northwestern Levantine during
February~10-March~15 and March~19-April~16, 1995.
Different cases are illustrated and compared among each other
as well as to a benchmark. For additional evaluation,
forecasts generated from the various
analyses are 
%qualitatively and quantitatively
compared to {\it in situ} and satellite data.
A back-of-the-envelope rule for estimating
the size of subspaces is also outlined.
Several properties of the multiscale
subspace estimation and analysis scheme are discussed.
The sensitivities to size confirm that the truncation to a
subspace is efficient and that the use of
quantitative criteria is important to determine adequate sizes.
The sensitivities to scales show that, for a given region,
smaller scales usually require larger subspaces, but that
dynamics can alter this property.
The sensitivities to the dynamical adjustment reveal that nonlinearities
can modify the dominant variability decomposition,
especially the dominant eigenvectors, and that
changes are largest for the features and regions with high shears.
Dynamical information is extracted from the estimated variability
variance fields, eigenvalue spectra, multivariate eigenvectors
and (cross)-covariance functions. Dominant dynamical balances
and the spatial distribution of hydrographic and velocity characteristic scales
are obtained for primary
regional features. In particular, the Ierapetra is found to be
close to gradient-wind balance and coastal-trapped waves are anticipated
to occur along the northern escarpment of the basin.

%\vskip 0.3truecm
%\noindent
%Keywords: 
%\item{}
}

\vskip 1truecm
\centerline{\bf Harvard Report in Physical/Interdisciplinary Ocean Science,
No.~60}
%No.~58}
%No.~57???}
%\centerline{Submitted to the {\it Journal of Atmospheric and
%Oceanic Technology}}
\centerline{December, 1998}

 
\vfill\eject
%\vskip 0.5 truecm
{\doublecolumns
\noindent
{\scaps 1. INTRODUCTION}
\vskip 0.3 truecm

In mapping geophysical fields and their uncertainties, it
is important to assess the sensitivity to {\it a priori}
parameters. The present work carries out such studies
for the {\it a priori} error and variability subspace estimation and
corresponding mapping methodology described in
a companion manuscript (Lermusiaux {\it et al.}, 2000, henceforth, LAL00).
With this three-dimensional, multivariate and multiscale approach (App.~A),
the dominant components of the {\it a priori}
error probability density functions of the fields or state to be mapped,
i.e.~the {\it a priori} error subspace,
are constructed in two parts. The observed portions of the
subspace are
specified based on measurements and multiscale statistical models.
The non-observed portions are then built by dynamical adjustments
to the observed portions.
Once this {\it a priori} error subspace is computed,
the data sets of interest can be mapped into gridded fields
and error properties by Bayesian estimation in the subspace.
Presently, the error statistics are assumed to be described
by covariance matrices and the mapping (analysis) is
reduced to a minimum error variance estimation problem. The 
%{\it a priori} error
subspace is then determined by the dominant
eigendecomposition of the {\it a priori} error covariance.

Two main objectives of this manuscript are
to investigate (i) the sensitivity of the state and error estimates
to the properties of the error subspace (e.g.~size, dynamics,
adjustment) and (ii) the sensitivity of
the subspace itself to its size and to the data and dynamics utilized in its
construction. The types of questions considered thus include:
How many, and why so many, eigenvectors should be employed to
represent most of a covariance matrix or to carry out a mapping?
What is the sensitivity of a variability or error subspace
to the parameters of the model employed
in the dynamical adjustments?
How, and by how much, should such
subspaces be dynamically adjusted? And finally,
what is the performance of forecasts generated from the various analyses?
The sensitivity studies are carried out by examining
the field and error estimates obtained from different {\it a priori} parameters.
For additional assessment, results
are also qualitatively and quantitatively compared to those
of the univariate horizontal (2D) scheme of
the Harvard Ocean Prediction System
(HOPS, e.g.~Lozano {\it et al.}, 1996; Robinson, 1996; Lermusiaux, 1999).
Since this benchmark (Sects.~3--4) is a multiscale analysis,
the present three-dimensional, multivariate and multiscale
mapping of all physical fields and their error covariances 
is referred to in short as the ``multivariate 3D'' scheme.
Ultimately, forecasts generated from the various analyses are
compared to {\it in situ} and satellite data,
using root-mean-square differences and
correlation coefficients for measures of skill. The third
objective is to extract some dynamical
knowledge from the variability
subspaces and fields computed.
%The questions considered
%include: how can this be done? and what kind
%of dynamical properties can be obtained?
This is carried out by studying the
spatial distribution of the estimated
hydrographic and velocity variances,
shapes and amplitudes of dominant multivariate eigenvectors,
and structures of (cross)-covariance functions.

The applications illustrated are in the Levantine sea
and focus on the mesoscale to subbasin-scale physical errors and fields.
The measurements were collected during the LIW95 campaign
(Roether {\it et al.}, 1996;
Malanotte-Rizzoli {\it et al.}, 1996),
carried out from January to April 1995 in the
Eastern Mediterranean by the
Physical Oceanography of the Eastern Mediterranean (POEM) group.
To introduce this region and its contemporary
upper-thermocline features, the results of a multivariate 3D estimation of
the potential density at 105~m and
velocity at 5~m are plotted on Fig.~1.
Note that several surface circulation features,
because of  their approximate geostrophic balance,
are easily distinguished along the steepest slopes of isopycnals.
To date, most investigations in the Levantine
that combine data and dynamics have resolved
the subbasin-scales and are usually based on
this close-to-geostrophy hypothesis
(e.g.~Hecht {\it et al.}, 1988; Milliff and Robinson, 1992;
\"Ozsoy {\it et al.}, 1993; Robinson and Malanotte-Rizzoli, 1993;
Malanotte-Rizzoli {\it et al.}, 1999).
%; and the references cited therein).

%The present work is concerned with several areas of research.
There are various techniques for the analysis of geophysical fields,
some of which also estimate uncertainties.
A relatively comprehensive overview of such techniques
and of issues involved in determining accurate {\it a priori}
error statistics
%for oceanic and atmospheric systems
is presented in LAL00 and references cited therein.
Determining (error) covariance or structure functions, and estimating
their parameters and covariance matrix form, is also
related to this work.
Recent studies include: in atmospheric science,
Rabier {\it et al.}~(1998), Riishojgaard (1998),
Dee and da Silva (1999), Gaspari and Cohn (1999), Gneiting (1999),
Franke and Barker (2000);
in physical oceanography,
Cressie and Huang (1999), Miller and Cornuelle (1999),
Menemenlis and Chechelnitsky (2000), Molinari and Festa (2000);
and, in biological oceanography, Abbott and Letelier (1998).
Today, error models, in particular the efficient representation of error
covariance matrices, are essential for error predictions
(Ehrendorfer, 1997; Lermusiaux, 1997, henceforth, LER97)
and data assimilation (Robinson {\it et al.}, 1998;
Chin {\it et al.}, 1999; Derber and Bouttier, 1999;
Lermusiaux, 1999a,b; Lermusiaux and Robinson, 1999;
Verron {\it et al.}, 1999;
Mitchell and Houtekamer, 2000; Tippett {\it et al}, 2000).

%The remainder of the paper is organized as follows.
In what follows,
Sect.~2 outlines a simple rule for estimating the size of error subspaces.
In Sect.~3, variations of the {\it a priori} error subspace
with its size and with the scales in the data,
and the changes in the {\it a posteriori} fields
as a function of these variations are illustrated for 
multivariate 3D analyses of temperature and salinity.
In Sect.~4, such sensitivity studies are
carried out, but for the complete physical fields:
effects of subspace size, scales considered
and nonlinearities in the dynamical adjustments
on the subspace itself and on the
{\it a posteriori} fields and errors
are investigated in detail.
In doing so, the subspaces are also utilized to 
help derive some dynamical 
properties (e.g.~characteristic scales, dominant balance, 
dynamical identity)
of geophysical features and phenomena estimated.
%In Sect.~4, the sensitivity of the {\it a priori} error subspace
%estimate to its size and to the nonlinear terms of the momentum
%equations used in the adjustments is first exemplified.
%The variations of the {\it a posteriori} fields
%$\xa$ with the properties
%of the error subspace (e.g.~size, linearly/nonlinearly
%adjusted) are then shown and studied. 
%The univariate 2D mapping
%scheme of HOPS is employed for comparisons and evaluations. 
In Sect.~5, the performance of forecasts initialized based
on the present field and error mapping scheme is evaluated.
The summary and conclusions are in Sect.~6.
In the Appendices, algorithms and notation are summarized,
relevant scalings and non-dimensional numbers
listed, and main parameters of the dynamical model overviewed.


\vskip 0.5 truecm
\noindent
{\scaps 2.~THE SIZE OF ERROR SUBSPACES: BACK-OF-THE-ENVELOPE
CALCULATION}
\vskip 0.3 truecm

A simple rule giving an approximate order of
magnitude for the size of multivariate error
subspaces is outlined. It is based on experience with
more rigorous and quantitative calculations (LER97, 
Lermusiaux, 1999a,b). It is exemplified for the case where
errors are the variability that occurs in a region
over specific time and space scales.

Consider first one state variable, e.g.~temperature or zonal velocity.
For the vertical variability,
2 to 3 empirical orthogonal functions (EOFs) or dynamical modes are often found
to explain most of the significant variance
(e.g.~De Mey and Robinson, 1987;
Haney {\it et al}, 1995; De Mey, 1997; Wunsch, 1997;
Pedder and Gomis, 1998;
von Storch and Frankignoul, 1998).
For the horizontal variability, most studies to date focus on
surface fields. In that case, the number of dominant
patterns or modes usually
varies between 15 to about 50: for example, consider a
surface field in the Gulf Stream or Atlantic region
(e.g.~Thacker and Lewandowicz, 1996, 1997;
Everson {\it et al.}, 1997; Kaplan {\it et al.}, 1997).
For the 3D variability, an approximate subspace size is simply the product of
vertical and horizontal requirements. With
the above values, this product yields a size of 30 to 150.

In the case of several state variables, a simple hypothesis
is to assume that the size of the subspace 
necessary to describe most of a normalized form of the 
multivariate variability
is proportional to the number of variables. 
For example, considering the primitive equation (PE) model\footnote{$^1$}
{\hsize=6.5truein
\eightpoint \hskip -23pt   
The PE model employed (App.~C) governs
five state variables, henceforth called the ``PE fields'': at any
time, it represents four three-dimensional fields,
the temperature $T$, salinity $S$
and internal horizontal velocities $\widehat u$ and $\widehat v$
(zonal $x$ and meridional $y$ directions),
and one two-dimensional field, the barotropic transport stream
function $\psi$.}
of HOPS, the approximate 
size of a mesoscale error
subspace is then 135 to 650 for the above values, which
is of the order of $10^2$ to $10^3$.
In summary, the order of magnitude of
the size $p$ of an error subspace
is often well approximated by,
$$
{\cal O}(p) \simeq \sum_{v} n_h^v \times n_z^v~,
\eqno(1)
$$
where $n_h^v$ ($n_z^v$) is an estimate of the horizontal
(vertical) requirements for the state variable~$v$.
Similarly, a range of likely $p$ is obtained
from the products of the minimum
and maximum horizontal and vertical requirements,
i.e.~$\sum_{v} \min (n_h^v \times n_z^v)\leq p \leq
\sum_{v} \max (n_h^v \times n_z^v)$.
Such back-of-the-envelope rules don't require
extensive 3D computations: values can be quickly obtained
for the specific region and processes
considered. They provide quick estimates of the
number of state samples, breed vectors or singular vectors
(e.g.~Molteni {\it et al.}~1996; Toth and Kalnay, 1997;
Barkmeijer {\it et al.}~1998)
that are needed to explain most of the variance
in an analysis or ensemble forecast, and
are thus very useful today.

The above sizes, even though large when compared to classic 1D
requirements, are small when compared to the
number of discrete state variables. To exemplify this property,
consider a covariance matrix $\BB$ for a horizontal,
univariate tracer field.
In its non-dimensional form (from here on denoted by $^*$),
the covariance function employed is a ``Mexican hat''
(e.g.~Carter and Robinson, 1987; Louis {\it et al.}, 1997),
i.e.~the negative of the second derivative of a two-dimensional
Gaussian function.
Its scalable parameters are the zonal $l_a^x$ and meridional $l_a^y$
zero-crossings, zonal $l_b^x$ and meridional $l_b^y$
$e$-folding decay-scales and
decorrelation time $\tau$. 
The elements of $\BB^*$ are then of the form
$C^*({\Br_1,\Br_2})=(1-a^2)\;e^{-b^2/2}$,
where $a^2 = (\Br_1-\Br_2)^T \BL_a^{-2} (\Br_1-\Br_2)$
and $b^2 = (\Br_1-\Br_2)^T \BL_b^{-2} (\Br_1-\Br_2)+
\left({\Delta t_i/ \tau}\right)^2$ are scalars,
and $\Br_1$ and $\Br_2$ denote horizontal position vectors,
with $\Br = (x,y)$.
The 2-by-2 anisotropic scaling matrices contain
the local zero-crossings,
$\BL_a={\bf diag} (l_a^x,l_a^y)\,$, and local length-scales,
$\BL_b={\bf diag} (l_b^x,l_b^y)\,$.
The scalar $\Delta t_i$ is the
interval between the data time $t_i$
and estimation time $t_0$ (App.~A).
The discretized domain utilized covers the
northwestern Levantine ($500$~km in the $x$ and
$520$~km in the $y$ direction) with a uniform resolution
of $10$~km (App.~C). For one 2D tracer field,
the size of $\BB^*$ is thus $[2073 \times 2073]$.

The normalized cumulative
eigenvalue spectrum of $\BB^*$ corresponding to subbasin-scales
($l_a^x = l_a^y = 200$~km, $l_b^x = l_b^y = 100$~km and
$\tau=\infty$) and mesoscales
($l_a^x = l_a^y = 60$~km,
$l_b^x = l_b^y = 30$~km and $\tau=\infty$) are given
on Fig.~2 (for convenience, scales are chosen
homogeneous).
% in space).
Figure 3 shows a few horizontal,
non-dimensional eigenvectors for the subbasin-scale case
(mesoscale eigenvectors are in LER97 and LAL00).
For the subbasin-scales,
15 eigenvectors represent $85\%$ of the variance,
36 vectors 99\%, and 50 vectors 99.8\%. 
For the mesoscales, larger
numbers are required: 15 eigenvectors represent
20\% of the variance, 50 vectors 56\%, and 100 vectors 83\%.
For at least 90\% of the variance (see Fig.~2),
19 vectors are necessary for the subbasin-scales
and 128 for the mesoscales, which is
only 0.7\% and 4.7\% of the number of
discrete state variables, 2703, respectively.
Note that in the statistical curves of Fig.~2,
all realizations of the 2D subbasin-scale/mesoscale fields 
are accounted for. However, in the real ocean,
during a specified time scale, the dynamical evolution is a constraint.
%not all realizations occur due to dynamical constraints.
Gyres, meanders or eddies are not everywhere
at the same time. Hence, for a given percentage of variance,
purely statistical sizes (Fig.~2) are usually
upper-bounds on the dynamical size needed at a given time. 

\vskip 0.5 truecm
\noindent
{\scaps 3. MULTIVARIATE THREE-DIMENSIONAL ANALYSIS OF TRACER FIELDS:
SUBSPACE SIZE}
\vskip 0.3 truecm

Hydrographic data (dynamical tracers $T$ and $S$) are
frequent observations in oceanography. 
The sensitivity of the {\it a posteriori} state $\xa_\trc$
to the size of the {\it a priori} error subspace is thus first
exemplified for the case of a global mapping of
temperature and salinity fields.
The {\it a priori} error subspace is constructed assuming
that the {\it a priori} error covariance of
these tracers $\BB_{\trc}$ is ``observed'' in the sense of LAL00,
i.e.~historical (synoptic) data are available to specify
its dominant eigendecomposition $\BB_{\trc}^p$. 
Variations of $\BB_{\trc}^p$ with
the scales considered in the data are illustrated.
Some differences between the
multivariate 3D scheme and Kalman
update (see App.~A.4) are also exemplified.

The data to be gridded are 146 $T$ and $S$ profiles. They
originated from 49 CTD (Conductivity-Temperature-Depth) and
97 salted AXBT (Aircraft-Deployed Expendable Bathythermographs)
observations gathered in the Levantine basin during February 10--18 
and on March 15, 1995, respectively\footnote{$^2$} {\hsize=6.5truein
\eightpoint
\hskip -23pt The data manipulations, which include
a few horizontal and vertical extensions,
and the estimation of the salinities of the AXBTs, are
described in LER97.} (Fig.~4a).
Their 5~m mesoscale coverage is estimated on Fig.~4b.
The numerical domain (App.~C) contains
$2073$ horizontal grid points and $20$ vertical
levels (leading to $108120$ tracer state variables); its
bottom topography is plotted on Fig.~4c.
The 146 profiles lead to 5840 data-points (dimension of $\By^\o$,
eq.~A1b), i.e.~the number of $T,S$ scalar data
residuals on vertical levels.
In all examples to follow, the
first {\it a priori} estimate $\xb_\trc$
($w=0$ in Apps.~A.2--3) is the horizontal average of $\By^\o$.
The scales to be mapped range from subbasin-scales to mesoscales.

The univariate 2D mapping scheme of HOPS, employed as a benchmark,
assumes {\it a priori} that two independent scales are present
and that errors are homogeneous horizontally.
It uses a successive correction,
minimum error variance technique (e.g.~Daley, 1991; Lorenc, 1992):
the largest scales are gridded first and the resulting fields
constitute the background in the mapping of the smaller scales.
Successively for each scale, global horizontal analyses
of the scale-filtered $T$ and $S$ data are thus carried out 
at various levels. 
 
For fair comparisons with this benchmark, 
the present multivariate 3D scheme also follows
a two-scale successive correction approach ($w=1,2$ in Apps.~A.2--3)
and also assumes horizontally homogeneous {\it a priori} errors,
scale by scale. Two multivariate 3D mappings
(App.~A.3) are thus carried out successively.
For each scale, $\BB_{\rm trc}^p$
is constructed based on a statistical model fit to data
in the horizontal and on data EOFs in the vertical (App.~A, Table~A.a).
In this Section, the data used to do so is $\By^\o$.

\vskip 0.5 truecm
\noindent
{\bf 3.1 Subbasin-scale}
\vskip 0.3 truecm 

{\capsfn Univariate 2D benchmark}.~Figure 5 shows four
subbasin-scale ($w=1$) univariate 2D analyses for $T$ and $S$.
The {\it a priori} horizontal covariances were specified
using a ``Mexican hat'' (Sect.~2),
of zero-crossings $l_a^x = l_a^y = 200$~km, decay-scales
$l_b^x = l_b^y = 100$~km and decorrelation time-scale
$\tau=50$~days, all estimated based on
the subbasin-scale data residuals.
The data error covariance at
data-points was assumed diagonal (at each level, $\BR=r\BI$).
The corresponding horizontally uniform, non-dimensional (0 to 1)
variance was calibrated to $r^*=0.03$, mainly
based on experience.

{\capsfn Multivariate 3D scheme}.~The 
subbasin-scale ($w=1$) multivariate 3D results
are summarized by Figs.~6--7 (construction of
$\BB_{\rm trc}$) and Figs.~8--10 (mappings).
The dominant eigendecomposition of $\BC^z_{\trc}$
(eq.~A10) is illustrated (Fig.~6) by its normalized cumulative
variance and four dominant $T,S$ eigenvectors, i.e.~four EOFs
of the subbasin-scale data residuals (eqs.~A7--A9).
The first multivariate vector ($64\%$ of the variance)
is surface-intensified,
vanishing around 500--600~m, with
a slightly larger amplitude in $T$ than in $S$.
For $T$, the maximum is around $150$~m.
For $S$, it is more pronounced and deeper, near 225~m.
The dominant subbasin-scale $S$~variability thus
occurs at intermediate water depths, below that of $T$.
On average, the vector best explains a vertical projection
of the $T,S$ subbasin-scale
variability in the northwestern Levantine: 
e.g.~a displacement, strengthening or weakening of the
Rhodes Gyre (e.g.~Milliff and Robinson, 1992).
The second vector (25\% of the variance) exhibits $T,S$ values
of opposite signs (adding effects in density),
a larger amplitude in $T$ than in $S$ and a
sustained decrease to zero around 225~m,
depth of the $S$ maximum on the first vector. It relates to
%atmospheric forcings and 
upper-layer variability during February--March 1995
(Malanotte-Rizzoli {\it et al.}~1996; LER97).
The third vector explains 5\% of the 
variance, with $T,S$ in phase from the surface down
to 700--800~m depth (compensating effects in density). Its
shape is such that it can for example
account for a tilting or stretching/squeezing of the 
cyclonic gyres and anticyclonic cone-shaped vortices (Brenner, 1993).
%and for the corresponding convections and mixing events.
The fourth vector only explains 2\% of the variance. 
%It may relate
%to (de)-couplings between the variations of tracers
%within 0~to~225~m and within 225~to~1000~m depths.

The horizontal covariance $\BC^{{\Br}^{\scriptstyle *}}_{\trc}$ (eq.~A11)
is computed based on the subbasin-scale Mexican hat 
used in the univariate 2D analyses.
From (eqs.~A13--15), the subbasin-scale vertical (eq.~A10)
and horizontal (eq.~A12) decompositions are then combined.
The resulting normalized
variance is explained by a few hundred 3D vectors,
as shown on Fig.~7.
Using (eq.~A16) in (App.~A.3), the subbasin-scale
mapping ($w=1$) is finally carried out.

Figure 8 plots parts of the state $\xa_\trc$ for the case where
$p=500$ in (eq.~A16), i.e~$\BB_{\trc}$
in (eq.~A13) is limited to its dominant 500 eigenvectors
which explain 99.6\% of the {\it a priori} error variance (Fig.~7). 
The size of $\xa_\trc$ being 108120,
this analysis is $216$ times cheaper than a full covariance formulation
(App.~A.4) and about $12$ times
cheaper than a Representer method (Bennett and Chua, 1994),
the number of Representers being 5840. For example,
the elapsed-time was $15$ minutes on a Sun Ultra and
the full formulation would have thus lasted about 2.25 days
(see LAL00 for more on computational costs). 
Comparing results (Figs.~8 and 5), the $T,S$ maps
are alike. For the subbasin-scale,
the assumptions of zero vertical covariance and of
uncorrelated $T$ and $S$ error fields does not appear
restrictive. This is mainly because there is 
here more than enough data coverage for the subbasin-scale.
The vertical and $T$--$S$ correlations are as expected stronger
with the 3D scheme, but the differences between fields
of Figs.~5 and 8 are within data error bounds.
%Without additional data,
%the best of the two estimates can not be determined.

Figure 9 shows the subbasin-scale $\xa_\trc$ for the case where
$p$ is only 250 in (eq.~A16).
The relative Root Mean Square (RMS) 
differences between the $T,S$ maps of Fig.~8
(500 vectors) and Fig.~9 (250 vectors) are less than $1\%$.
This is an interesting result. Once most (e.g.~here $\ge 98\%$,
see Fig.~7) of the error covariance has converged,
the remaining eigenvectors can be ignored if their eigenvalues 
are negligible in comparison to data error eigenvalues (see App.~A.4).
Considering costs, the mapping with $p=250$ would be
$10^3$ times cheaper than the full covariance mapping.
Finally, using $10$ vectors as on Fig.~10 is not sufficient. 
Indeed, for 99\% of the variance, the quick product of 
the horizontal and vertical requirements (Sect.~2) indicates here
a total of $36\times 8=288$ vectors
(see eq.~1 and use Fig.~2a and Fig.~6),
i.e.~an order of $10^2$ multivariate 3D vectors.
This back-of-the-envelope calculation agrees with the 
numerically computed spectra of Fig.~7.

\vskip 0.5 truecm
\noindent 
{\bf 3.2 Mesoscale} 
\vskip 0.3 truecm
 
{\capsfn Univariate 2D benchmark}.~Figure 11 shows four $T$ and $S$
mesoscale ($w=2$) univariate 2D analyses. 
%at the same levels as on Fig.~5.
Zero-crossings were estimated based on data
at $l_a^x = l_a^y = 60$~km,
decay-scales at $l_b^x = l_b^y = 30$~km and 
decorrelation time-scale at $\tau=7$ days.
The non-dimensional data error variance $r^*$ 
for $\BR=r\BI$ was calibrated to 0.03, as in Sect.~3.1.

{\capsfn Multivariate 3D scheme}.~In the mesoscale
3D mapping, $\xb_\trc$ is the subbasin-scale state of Fig.~8. 
The construction of the
dominant mesoscale covariance
$\BB_{\rm trc}$ is illustrated by Fig.~12.
In the vertical (Fig.~12a), the 
variance of the mesoscale data residuals (eq.~A10)
decays less rapidly with the eigenvalue number than
in the subbasin-scale case (Fig.~6).
In the horizontal, $\BC^{{\Br}^{\scriptstyle *}}_{\trc}$
(eq.~A11) is computed based on the mesoscale Mexican hat 
used in the univariate 2D analyses.
From (eqs.~A13--15), the mesoscale vertical (eq.~A10)
and horizontal (eq.~A12) decompositions are then combined.
As shown on Fig.~12b, the resulting normalized
mesoscale variance is mostly explained by a few 3D vectors,
but compared to the subbasin-scale (Fig.~7),
slightly more vectors are needed.
Using (eq.~A16) in (App.~A.3), the mesoscale
mapping ($w=2$) is finally carried out.

Figures 13--14 show $\xa_\trc$ obtained 
using the dominant $500$ and $750$ mesoscale vectors, respectively.
Differences with the fields of Fig.~11
are small, at several locations within data error bounds.
However, the 3D scheme which uses
the significant EOFs of the mesoscale residuals (eqs.~A7--10)
yields larger vertical correlations 
in $\xa_\trc$ than the 2D scheme does.
For example, in the upper 0~to~50$\,$m mixed-layer,
the measurement noise and
unresolved vertical scales are filtered more on
Figs.~13--14. Small mid-depth features
of significant vertical extension, e.g.~see $T,S$ around 225~m,
also appear better represented.
From our experience, for mesoscale mappings with
data resolutions similar to those of Fig.~4a,
a disadvantage of the univariate 2D 
analysis is in fact its sensitivity to
environmental noise.
Finally, considering mesoscale error subspace sizes,
the relative RMS differences between maps
of Figs.~13 and 14 are less than 3\%.
Even though 750 vectors explain 5\% 
more {\it a priori} mesoscale
variance than 500 vectors do (Fig.~12b),
most of the additional 250 vectors have eigenvalues
smaller than the local data error variance. They 
thus mainly yield mesoscale corrections within data error bounds.

For the 3D scheme, two main properties
(see App.~A.4) have been exemplified.
First, once most of the {\it a priori} error covariance has converged
(e.g.~99\%), the remaining eigenvectors can be truncated
in most cases. Second,
the {\it a priori} error vectors of eigenvalues negligible
with respect to data errors 
have negligible impacts on the {\it a posteriori} state.
Even if the {\it a priori} error covariance has
not converged, these vectors can usually be omitted.

\vskip 0.5 truecm
\noindent
{\scaps 4. MULTIVARIATE THREE-DIMENSIONAL ANALYSES OF PE FIELDS}
\vskip 0.3 truecm

Primitive equation fields and their errors are now considered.
The examples shown are still for
%the subbasin-scales to mesoscales in
the northwestern Levantine basin, but the focus is on the period
from March 19 to April 16, 1995
%after March 15, 1995 
(Sect.~3 considered tracer fields and
only used data collected before and on March 15).
Four mapping cases are discussed (Table~1):
one univariate 2D benchmark and three variations
of the multivariate 3D scheme.
\bigskip

\eject}

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf
 
\centerline{Table 1. The four analyses selected and their main differences}
\vskip -0.7truecm
$$\vbox{\halign{
#\hfil&\quad#\hfil&\quad#\hfil&\quad#\hfil&\quad#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
Case& Scheme Employed & Size of {\it a priori} &
Velocity Covariances & Velocity Dynamics Utilized\cr
&& Error Subspace \cr
\noalign{\vskip2pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip3pt}
\cr
~ 1: & Univariate 2D (benchmark) & full covariances
& Not built: direct integration & Geostrophy\cr
& & (all 2D vectors) & of thermal-wind equations\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
~ 2: & Multivariate 3D & $p=322$ & Built via an ensemble of & Linear momentum PE model\cr
& & & adjustment PE integrations\cr
~ 3: & Multivariate 3D & $p=322$ & \qquad ~~ As above &
Nonlinear momentum PE model\cr
~ 4: & Multivariate 3D & $p=400$ & \qquad ~~ As above &
Nonlinear momentum PE model\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
\vskip-0.5truecm
{\doublecolumns
The univariate 2D benchmark (Case 1) first carries out
horizontal analyses of the $T$ and $S$ data.
This is done in two successive scale corrections, as in Sect.~3.
%the gridded subbasin-scale tracer fields, result of
%the first stage, are the background in the second (mesoscale) stage.
The total velocity is then computed
by assuming geostrophic balance with the
gridded tracers and integrating the thermal-wind equations
up and down from a level of reference.

For fair comparisons, the present multivariate 3D scheme (Cases~2--4)
also uses a two-scale successive correction approach
($w=1,2$ in App.~A.2) to
construct the {\it a priori} error subspace.
For each scale, the tracer portion $\BB_{\trc}$ of the
{\it a priori} error covariance $\BB$ is specified as in Sect.~3.
The non-observed total velocity portions are then built in 
statistical and dynamical accord with $\BB_{\trc}\,$,
through an ensemble of adjustment momentum integrations 
(App.~A, Table~A.b).
These integrations lead to an ensemble of PE adjusted fields.
Variability samples, i.e.~the differences between these fields
and the {\it a priori} state $\xb$, are then computed,
normalized and organized by singular value decomposition (SVD). 
Once a convergence criterion determines that
the number of samples is large enough 
to explain most of the variance,
an estimate $\BB^p$ of the dominant components of $\BB$
is obtained. In building this covariance,
for the subbasin-scales, $\xb$ is
the horizontal average of the tracer data, with null velocity;
for the mesoscales, $\xb$ is
the result of the subbasin-scale mapping.

The multivariate 3D cases of Table~1 were selected 
to illustrate the main tendencies in an extensive set of computations.
They exemplify the sensitivities of the error subspace and state
$\xa$ to: (i) the nonlinear terms of the
momentum equations used in the adjustment of the non-observed
velocity variability to the observed tracer variability (Cases~2 and 3),
and (ii) the size of the subspace, i.e.~the number $p$
of {\it a priori} error eigenvectors utilized
in the mapping (Cases~3 and 4).
%; the rationale for the values 322 and 400 is in Sect.~4.2).
Case 4 is the mapping briefly exemplified in LAL00.

In the remainder of this Section,
the PE field and error covariance estimates of the four cases
are compared and studied.
For conciseness, the subbasin-scale mappings are briefly
overviewed in Sect.~4.1;
the detailed presentation is restricted to 
the mesoscale stage $w=2$.
The estimates of the {\it a priori} error covariance for
the mesoscale fields,
$\BB^p= \BE \BPi \E^T$, are compared in Sect.~4.2;
the {\it a posteriori} fields, in Sect.~4.3;
and the {\it a posteriori} mesoscale error covariances,
$\Bap=\BE^\a\BPi^\a\EaT$, in Sect.~4.4.
The dynamical possibilities raised by the study of 
error and variability subspaces
are emphasized in Sects.~4.2 and 4.4.

\vskip 0.5 truecm
\noindent
{\bf 4.1 Subbasin-scale}
\vskip 0.3 truecm

For each case (Table~1),
the data $\By^\o$ employed for the subbasin-scale stage
are 512 profiles of the MODB (Med5) winter
climatology for the Levantine (Brankart, personal communication)
and the T/S profiles
obtained from the {\it in situ} CTD and AXBT data
gathered during February 10--18 and March 15, 1995 (Fig.~4a).
For coherence with the recent evolution
of the deep Eastern Mediterranean waters (Roether {\it et al.},~1996), the 
%climatological
MODB data were cut at 1250~m depth: only {\it in situ} data are
utilized below 1250~m. In all cases,
$\xb$ is in hydrostatic equilibrium: its tracer 
fields are the horizontal averages of 
the data $\By^\o$ and its velocity is null.

For the univariate 2D scheme (Case 1), 
the tracer fields are estimated from $\By^\o$ as
in Sect.~3.1, using the same {\it a priori}
horizontal error covariances of $T$ and $S$.
The total velocity fields are assumed in
geostrophic balance with these gridded tracers,
from a level of no motion at 600~m (since the 2D scheme is
univariate and geostrophy is linear, velocities are only
computed after the mesoscale tracer analyses, Sect.~4.3).
For the multivariate 3D Cases~2--4, 
the $\BB_{\trc}^p$'s are also
as in Sect.~3.1, except for the variations of the size $p$.
Three estimates of the
subbasin-scale $\BB^p$ are then built as indicated in
Table~1, using the PE model of App.~C.
With these $\BB^p$'s,
three mappings of the data $\By^\o$ for this first stage
are finally carried out, following (App.~A.3).

\vskip 0.5 truecm
\noindent
{\bf 4.2 Mesoscale: data and {\it a priori} error covariance}
\vskip 0.3 truecm

%During the quality control, 4 profiles were taken off for this mapping. 
%Hence, the 236 CTD's used here + the 49 used before + the 4 missing
%lead to the total of 289: 36+49+4=289
%see /home/scotia4/pierrel/liw/data/hydro/README)

For each case (Table~1), the data to be gridded
for the mesoscale correction are 
236 CTD profiles gathered
from March~19 to April~16 (Fig.~15).
They lead to 9440 data-points (dimension of $\By^\o$).
%the number of $T,S$ scalar observations on vertical levels.
The field estimation is centered on March~27.

The four $\xb$'s (Cases~1 to 4) are the four subbasin-scale fields 
of Sect.~4.1; they contain no
mesoscale by construction. Their {\it a priori} error is thus
mesoscale variability: i.e.~for Cases~2 to 4, $\BGamma=\BI$ in 
(eq.~A23). The observations used to estimate the
corresponding covariances consist of
the 289 CTD profiles gathered from Feb.~10 to April 16;
they include the data $\By^\o$ and are denoted 
by $\By^\o_{\rm ctd}$. These four covariance estimations are
now described.

{\capsfn Univariate 2D benchmark}.~The {\it a priori}
horizontal error covariances of $T$ and $S$
are specified using a Mexican hat, as in Sect.~3.2.
Its zero-crossings and decay-scales 
are calibrated on basin average to the mesoscale residuals
of the 289 CTD profiles, which yields $l_a^x = l_a^y = 60$~km and
$l_b^x = l_b^y = 30$~km, respectively.
The non-dimensional data error variance $r^*$
for $\BR=r\BI$ is adjusted to 0.03.
Velocity errors are not computed;
velocities are directly estimated from geostrophic balance
with a level of no motion at 600~m.

{\capsfn Multivariate 3D scheme}.~In the construction of
the mesoscale $\BB^p= \BE \BPi \E^T$,
%the mesoscale tracer error covariance $\BB_{\trc}$ 
%(size $[108120 \times 108120]$) is first assumed
%observed in the sense of LAL00 and 
the dominant decomposition of its tracer submatrix,
$\BB_{\trc}^p = \BE_{\rm trc} \BPi_{\rm trc} \BE_{\rm trc}^T$, 
is computed first (App.~A, Table~A.a). This is done by combining
the vertical EOFs of the mesoscale tracer residuals
($\By^\o_{{\rm ctd}_{\scriptstyle i}} - {\cal H}_i(\xb)$, $i=1,\ldots, 289$,
eq.~A7) with the eigendecomposition
of an analytical horizontal tracer covariance (of elements specified using
the mesoscale Mexican hat of the benchmark, see above).
The resulting cumulative error variance
%as a function of the number of the 3D tracer eigenvector
is a curve similar to that of Fig.~12b
(it differs slightly because of the different $\xb$ and
data $\By^\o_{\rm ctd}$ employed).
For the sizes $p$ utilized (Table~1),
the dominant 322 tracer eigenvectors explain
67.2\% of the total tracer variance, the dominant 400, 72\%,
and dominant 500, 76.5\%.
The complete $\BB^p$'s are then obtained from $\BB_{\trc}^p$, by
adjustment momentum integrations (App.~A, Table~A.b).
%In the present cases, the $q$ perturbed tracer fields were created 
%using $\Bu^j={\sqrt q} \Be^j$.
%
The $j=1, \ldots, r$ dynamical adjustments lasted for 2 model-days,
with $r$ increasing based on the convergence criterion 
$\rho \ge \alpha$ (eq.~A22), with the limit $\alpha$ set to 0.93.
This coefficient $\rho$ assesses the added value of new parallel batches of
adjustment runs; some of its values are listed in Table~2.a.

%Initial location of Table 2
%\noindent
For Case~2, $\rho \ge \alpha$ was achieved after 
322 runs, for $p=322$ and $\tilde p = 288$.
For Case~4, 400 runs were necessary, based on
$p=400$ and $\tilde p = 360$. This larger dimension is due to the
nonlinear momentum terms which create flow variability 
not estimated by the linear adjustments of Case 2.
Since Case 3 is simply the nonlinear extension of Case 2,
convergence is not an objective.
%and it is not exactly achieved ($\rho_{288-322}=0.928$).
Overall, the rate of increase of $\rho$ diminishes
as $\tilde p \le p$ grows: this is partly
because the singular values of $\BB$ decay rapidly.
%and because $\tilde p$ becomes $\sim  p$.

Table 2.b estimates the overall
similarity among the three $\BB^p$'s, using (eq.~A22).
Comparing Cases~2 and 3, averaged on volume and state variables,
nonlinear momentum terms only change
about 2.7\% of the global, non-dimensional PE mesoscale variability.
In the present examples, since $\Bu^j={\sqrt q} \Be^j$ in (eq.~A17),
the tracer covariance in a $\BB^p$ is exactly
$\BB_{\trc}^p$. Therefore, the
covariances involving total velocities in 
the $\BB^p$'s of Cases~2 and~3
differ by about $(1-\rho) n/n_{\rm vel}=5.4\%$,
where $n$ is the state vector size (218,943) and 
$n_{\rm vel}$ the number of velocity variables (110,823).
This $5.4\%$ difference is an indicator of
the volume-averaged importance of the nonlinear momentum
terms in estimating the total velocity covariance.
Its local variations will be exemplified in Sect.~4.2.a.
Comparing Cases~3 and~4, they are about
13.2\% apart. The tracer
eigenvectors of high number (small eigenvalues) here still
affect the error subspace estimate\footnote{$^3$}
{\hsize=6.5truein
\eightpoint \hskip -23pt 
This 13.2\% rise is partly due to the 4.8\% increase
in normalized 3D tracer variance explained
(from 67.2\% at $p=322$ to 72\% at $p=400$), see above.}.
%This tracer effect is here accentuated 
%by the numerical PE model employed which can sometimes
%overestimate of the covariances involving total velocities.}.
Finally, comparing Cases~2 and~4, their differences combine the
nonlinear and size effects (2.7\% and 13.2\%).
%\bigskip

\eject}

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf

\centerline{Table 2. Convergence of and similarity between error subspaces}
\vskip -0.5truecm
$$\vbox{\halign{
#\hfil&\quad#\hfil&\quad#\hfil&\quad#\hfil&\quad#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip3pt}
\cr
a. Convergence coefficient \cr
\hfill ~~ for Case 2: & $\rho_{254-288}=0.919$ & $\rho_{288-322}=0.930$\cr
\hfill ~~ for Case 3: & same values as for Case 4.\cr
\hfill ~~ for Case 4: &
$\rho_{288-322}=0.928$ & $\rho_{322-360}=0.929$ & $\rho_{360-400}= 0.935$\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
b. Similarity coefficient \cr
\hfill ~~ between Cases 2--3: & $\rho=0.973$\cr
\hfill ~~ between Cases 2--4: & $\rho=0.845$\cr
\hfill ~~ between Cases 3--4: & $\rho=0.868$\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
\vskip-0.5truecm
%\noindent

{\doublecolumns
In what follows, the three $\B^p= \BE \BPi \E^T$'s
(Cases~2--4) are further intercompared by discussing
and studying their error variances (Sect.~4.2.a),
covariance eigenvectors (Sect.~4.2.b)
and covariance functions (Sect.~4.2.c).
These estimates give an important advantage
over the univariate 2D scheme (Case 1)
which only provides 2D error variance fields.

\vskip 0.3 truecm
\noindent
{\scaps 4.2.a A priori mesoscale error variance}

{\capsfn Univariate 2D benchmark}.~The 
{\it a priori} tracer error variances are set
to the horizontally-averaged
variance of the residuals of
$\By^\o_{\rm ctd}$: they are depth-dependent.
The corresponding velocity error variances
% are not computed, but they 
are set to these of velocities in geostrophic balance
with the tracers.

{\capsfn Multivariate 3D scheme}.~The error variances form
the diagonal of the $\BB^p$'s. For fair comparisons
with the benchmark, tracer variances are
% by construction
horizontally uniform in the data domain (Fig.~15).
For Cases 2 and 3 ($p=322$),
the $T$ ($S$) error standard deviation at 5~m is
$0.40^\circ$C (0.080~PSU), while at 500~m, it is
$0.17^\circ$C (0.039~PSU). 
For Case~4 (see LAL00), standard deviations
are logically slightly larger:
to first order\footnote{$^4$}
{\hsize=6.5truein
\eightpoint \hskip -23pt
Assuming $x_2 > x_1 > 0$,
to link ratios of standard deviations
to ratios of variances, and inversely, one can use: 
${x_2^2 - x_1^2 \over 2 x_2^2} \le
{x_2 - x_1 \over x_2} \le
{x_2^2 - x_1^2 \over 2 x_1^2}$.
%i) if ${x_2 - x_1\over x_2} = \beta$, then
%${x_2^2 - x_1^2 \over x_2^2} \le 2 \beta$;
%and, ii) if $x_2^2 - x_1^2= \alpha$,
%then ${\alpha \over 2 x_2^2} \le
%{x_2 - x_1 \over x_2} \le
%{\alpha \over 2 x_1^2}$.
}, using the variance numbers given above,
by 3.3\% to 3.6\% ($100 \times (72-67.2)/2$, divided
by 72 and 67.2).
The total velocity variances are however non-uniform in 3D.
By construction, large {\it a priori} errors (Table 3)
correspond to the locations and features (Fig.~1)
estimated to have large mesoscale flow variability.

%Initial location of Table 3
%\noindent
For Case 2 (linear momentum),
at the surface, the largest standard deviations
are as expected above the Rhodes basin (flanked by steep topography)
within the deep, subbasin-scale Rhodes Gyre.
Other features of above average
flow variability are, in decreasing
order, the Ierapetra, Mersa Matruh Gyre (one of
its lobes is in the domain) and
central branch of the Mid-Mediterranean Jet.
At 500~m depth, the largest $\widehat u, \widehat v$ deviations
are near the inflow/outflow
of the Kasos and Karpathos Straits,
and again above the Rhodes basin.

For Case 3 (nonlinear momentum, $p=322$), at the surface,
it is not the Rhodes basin, but the high velocity, smaller size
Ierapetra which is estimated to dominate.
With respect to Case 2, the nonlinear correction in
standard deviation for the Ierapetra is 50\% for the surface 
internal velocity and 10\% for the barotropic transport (Table~3).
For the larger Mersa Matruh Gyre, the correction is smaller,
especially in $\psi$ (30\% more for
$\widehat u, \widehat v$, but only 5\% more for $\psi$).
For the Mid-Mediterranean Jet and
subbasin-scale Rhodes Gyre, the nonlinear terms barely increase the 
barotropic and surface flow variabilities.
At 500~m depth, as for Case 2, the two Straits and Rhodes
basin dominate the internal deviations,
with a nonlinear correction of 3\% and 7.7\%, 
respectively. However, Case 3 also
estimates that the Ierapetra and an anticyclone of the
Mersa Matruh-Shikmona Gyre (MM-S-G) complex
have above average internal deviations at 500~m.
For these deep mesoscale vortices,
nonlinear terms are thus important at depth.
Overall, the average increase in standard deviation is 3\% for $\psi$,
8\% for the surface $\widehat u, \widehat v$, and 4.6\%
for the 500~m $\widehat u, \widehat v$ (percentages relative to Case 2).
These upper-ocean computations agree with the 
volume-averaged velocity variance increase of 5.4\% 
(estimated above from Table~2.b).
For Case 4 (nonlinear momentum, $p=400$), results
are similar to these of Case 3, except that
standard deviations are a bit larger (LAL00),
to first order and overall, by about 6.6\% to 7.6\%
($13.2/2$ to $13.2/2/0.868$, see footnote~4).
%\bigskip

\eject}

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf

\centerline{Table 3. {\it A priori} mesoscale error
standard deviations (St.Dv.) for velocity:
estimates of Cases 2 and 3/4}
\vskip -0.5truecm
$$\vbox{\halign{
#\hfil&\quad#\hfil&\quad#\hfil&~#\hfil&\quad#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
Case& Depth & Features (and horizontal averages) &
$\psi$ error St.Dv.&Internal velocity error St.Dv.\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip3pt}
\cr
~ 2. & Surface & Above the Rhodes basin: & 0.85~Sv & \quad 15.5~\cmsecm1\cr
& (5~m) & Ierapetra: & 0.8~Sv & \quad 14.5~\cmsecm1\cr
& & Mersa Matruh Gyre: & 0.8~Sv & \quad 13~\cmsecm1\cr
& & Mid-Mediterranean Jet: & 0.68~Sv & \quad 13.5~\cmsecm1\cr
& & \hfill Surface average: & 0.65~Sv & \quad 12.5~\cmsecm1\cr
\noalign{\vskip3pt}
& 500~m & Kasos and Karpathos Straits: & & \quad 1.7~\cmsecm1\cr
& & Above the Rhodes basin: & & \quad 1.3~\cmsecm1\cr
& & \hfill 500~m average: & & \quad 0.65~\cmsecm1\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip2pt}
~ 3./4.  & Surface & Ierapetra: & 0.87/0.9~Sv & \quad 22/23~\cmsecm1\cr
& (5~m) & Above the Rhodes basin: & 0.87/0.9~Sv & \quad 16/16~\cmsecm1\cr
& & Mersa Matruh Gyre: & 0.84/0.86~Sv & \quad 17/18~\cmsecm1\cr
& & Mid-Mediterranean Jet: & 0.69/0.72~Sv & \quad 14/15~\cmsecm1\cr
& & \hfill Surface average: & 0.67/0.7~Sv & \quad 13.5/14~\cmsecm1\cr
\noalign{\vskip3pt}
& 500~m & Kasos and Karpathos Straits: & & \quad 1.75/1.8~\cmsecm1\cr
& & Above the Rhodes basin: & & \quad 1.4/1.5~\cmsecm1\cr
& & Ierapetra: &  & \quad 0.95/1~\cmsecm1\cr
& & Anticyclone in MM-S-G complex:~ & & \quad 1/1.05~\cmsecm1\cr
& & \hfill 500~m average: & & \quad 0.68/0.7~\cmsecm1\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
\vskip-0.5truecm

{\doublecolumns
Two general remarks can be made: 
(i) for accurate flow variability variances, not only
the background density field, but also the regional topography,
nonlinear and diffusion terms (App.~C) should be considered;
(ii) even though subspaces are used,
the present global estimates of mesoscale variances
are in accord with and complement
previous local estimates (e.g.~\"Ozsoy {\it et al.}, 1993;
Robinson and Golnaraghi, 1993).

\vskip 0.3 truecm
\noindent
{\scaps 4.2.b A priori mesoscale error eigenvectors and spectrum}

{\capsfn Case 2}.~Figure 16 illustrates
the first and second non-dimensional singular vectors,
columns of $\BE^*$ in (eqs.~A21--23).
As expected from the variance results (Sect.~4.2.a),
for linear momentum dynamics,
the dominant mesoscale variability during early Spring
correspond to large-mesoscale variations of the Rhodes Gyre.
In the past few decades, the Rhodes Gyre has indeed been observed
to be the main subbasin-scale feature of the region 
(e.g.~Milliff and Robinson, 1992; Nittis and Lascaratos, 1998).
The two vectors shown explain 1.08\% and 1.05\%
of the variance explained by the 322 vectors, respectively.

The first vector (Panels 1a--b) covers the principal core of the Gyre,
above the Rhodes basin. Horizontally (Panel 1a), 
the main $T$ and $S$ structures are almost axisymmetric,
with horizontal decay-scales of about 40~km (radius of 60~km).
Their shapes are close to large-mesoscale Mexican hats.
%, which is expected by construction,
%even though only 322 PE vectors are employed.
On the open-ocean side,
the main $\psi$ structure
%(Panel 1a)
is also close to a
Mexican hat, indicating a predominant geostrophic equilibrium
(from the isocurve $\psi=0$,
the ratio of minimum to maximum $\Delta \psi$ is about 4,
as for the tracers).
Close to the coastline, the $\Delta \psi$'s vanish,
mainly because of decreasing depth and increasing coastal friction.
Vertically, the $T$ structure is surface intensified (Panel 1b),
with an uniform extremum in the upper 20 to 200~m depth
and a significant extension down to about 900--1000~m above
the Rhodes basin. The $S$ structure has a
subsurface extremum within 250~m to 320~m and
one significant zero-crossing near 100~m.
The signs of $T$ and $S$ structures reverse
at about 1500~m and 1300~m, respectively. Below,
the largest $T$ and $S$ amplitudes are each
near 2300~m, but these local extrema are
much smaller than the global upper-thermocline extrema
(about 100 and 20 times smaller, respectively).
Overall, where $T$ and $S$ are both relatively large,
they are in phase, compensating each other in density (at 5~m, they are
in opposition, but $T$ strongly dominates).
The surface $\widehat u$ (Panel 1a) and
cross-sections in $\widehat u$ and $\widehat v$ (Panel 1b)
confirm the predominant geostrophic equilibrium.
In addition to the two main lobes, a weak
topographic wave pattern is visible in the $\widehat u$ map
and sub-surface tracer cross-sections. Its average amplitudes are
about 20\% of those of the two main lobes.
Vertically, the $\widehat u$ and $\widehat v$ structures
%(Panel 1b)
reverse sign near 700~m. The extrema
below this reversal are at about 1500~m, but their
amplitudes are 10 times smaller than the surface ones.
Using $\psi$ (Panel 1a), a minimum in total velocity
is expected near 900~m.

The second vector (Panels 2a--b) mainly
relates to displacements of the core of the Gyre, along
the north-south axis of the Rhodes basin.
Except for the main double lobes in $T$, $S$ and $\psi$,
and corresponding triple lobes in $\widehat u$ and $\widehat v$,
most properties are as those of the first vector:
e.g.~$\psi$, $T$ and $S$
horizontally in phase where they are all important;
depths of zero-crossings for
$T$, $S$, $\widehat u$ and $\widehat v$ overall unchanged; 
predominant balance for $\psi$, $T$, $S$, $\widehat u$ and $\widehat v$ 
remains geostrophic equilibrium, etc. Nonetheless, such a
displacement of the Gyre's core is logically suggested to
interact more with topographic wave patterns (e.g.~Smith, 1983)
than in the first vector.
Along the northern escarpment, the dominant wavelength 
of the pattern (Panels 2a--b)
is about 160~km (scale of 25~km). Possible wave types and 
time scales are further discussed below.

The other dominant vectors (LER97) 
also relate to the main features (Fig.~1 and Table 3).
This is because to explain a specific variation of,
or interaction among, features,
a group of vectors is usually required.
Considering for example the vectors 3 to 16,
their largest amplitudes assign 
the vectors: (5--9,11--16) to
the Ierapetra; (3--7,9,12--16) to
the Rhodes Gyre; (3--6,13--14) to
wave patterns along the northern
escarpment and Mid-Mediterranean Ridge;
(12--14) to the Mersa Matruh Gyre;
(8--16) to the Mid-Mediterranean Jet;
(3--4,6,9--16) to the West Cyprus Gyre;
(4,6,9,12--16) to the Asia Minor Current; and
(9--10,12--16) to the Mersa Matruh-Shikmona Gyre complex.
%Mersa Matruh Gyre (
%lack of data: 5,7
Since each vector is constructed to be orthogonal to all the
more dominant ones and 3D dynamical
couplings can be complex, a detailed subspace analysis
requires additional computations (rotations, projections, wavelets, etc)
beyond the present scope.
Considering finally the rate of decay of the error spectrum, 
the error variance
explained by the 322 vectors is distributed such that
the dominant 50, 100 and 200 vectors explain 41.4\%,
66.7\% and 89.4\% of the total variance, respectively.

{\capsfn Case 3}.~Figure 17 is as Fig.~16, but for nonlinear 
momentum adjustments (Table~1).
In accord with the variance results (Table~3),
the feature of dominant mesoscale variability is then the Ierapetra.
The two vectors shown explain 1.59\% and 1.24\%
of the variance explained by the 322 vectors,
respectively. These percentages are larger than in Case 2 because
nonlinear terms steepen the start and flatten the upper-middle
of the eigenvalue spectrum.
The first vector corresponds to the vortex itself,
the second to eastward-westward displacements (for
the sign plotted, it is a westward motion).

In comparison with Case 2, the
horizontal $T$ and $S$ structures of the first vector (Panel 1a)
are, when significant, closer to Mexican hats 
(almost no wave patterns) and of
lesser width (radius near 50~km, decay-scale of 30~km).
The lobes far from the vortex
in the surface $S$ are mainly due to the truncation
and have a maximum amplitude that is
insignificant based on $\rho=0.928$ (Table 2.a).
They are excited by radiating nonlinear effects
(absent in Case 2)
and can be seen where amplitudes are relatively small.
Vertically (Panel 1b), $T$ and $S$ are more stretched
than in Case 2 (Fig.~16), as one would expect for a ring-like vortex.
The $T$ is still surface intensified, but its uniform extremum is
deeper, within 50 to 250~m. This agrees with previous results:
e.g.~the shallow surface mixed-layer of the Rhodes Gyre
(e.g.~Malanotte-Rizzoli {\it et al.}, 1996) and the structure
of Levantine warm core eddies (e.g.~Brenner, 1993). The dominant
feature in $T$ extends down to about 600~m,
main depth of the present Ierapetra,
instead of about 1000~m for the Rhodes Gyre (Fig.~16).
The main structure in $S$ has
a zero-crossing near 80~m and
a subsurface extremum within 250~m to 340~m.
The sign in $T$/$S$ reverses
at about 1300~m/1100~m, which is
shallower than for the Rhodes Gyre.
The extrema below are also shallower,
at about 1700~m for both $T$ and $S$, but they
remain smaller (50 and 20 times, respectively)
than the upper-thermocline extremum. Overall,
$T$ and $S$ are still in phase when they are both
important. The $\psi$ component (Panel 1a) is less
close to geostrophic equilibrium than on Fig.~16,
as shown for example by the asymmetries of its pattern.
The $\widehat u$ and $\widehat v$ components (Panels 1a and 1b) 
confirm this fact.
The nonlinear, inertial and viscous terms of the momentum
equations as well as topography modify the classic
geostrophic antisymmetries of $\widehat u$ and $\widehat v$
with respect to the $x$ and $y$ axes. The horizontal lobes
(Panel 1a) present sharp gradients and 
are tilted, with their heads wider than their tails.
The maximum velocity is reached near a radius of 30~km.
Vertically (Panel 1b), $\widehat u$ and $\widehat v$
reverse sign near 400~m, indicating a minimum in total
velocity near 600~m based on $\psi$ (Panel 1a). The 
$\widehat u$ and $\widehat v$ extrema
below this reversal are near 1600~m and deeper, but 
their amplitudes are about 10 percent of the surface ones.
The variations of the deep counter-flow and corresponding 
density feature are thus weak. However,
they have a large vertical extent (here 1500~m) and
can influence the stability of the Ierapetra,
as found by forecast sensitivity studies in LER97 
and argued in general by Dewar and Killworth (1995) and 
Killworth, {\it et al.}~(1997).

Considering the second vector on Fig.~17, 
displacements of the Ierapetra
can excite radiative momentum patterns  (Panels 2a--b).
Within the fixed mesoscale tracer variability, these
adjustment patterns logically correlate most to density variations
in quasi-geostrophic balance, hence
the $T$ and $S$ radiation patterns on
Panels 2a--2b. The horizontal $\widehat u$ 
is close to a ``45 degree shamrock'' which indicates a
thermal-wind balance (Daley, 1991). However, 
at all depths, $\widehat u$ and $\widehat v$ have clear 
asymmetries and amplitude variations which
result from local peculiarities of
the topography, density, nonlinearity and bottom friction
(see dynamics below). In the vertical, most 
properties are as those of the first vector.

Since the two dominant vectors are modified,
the others also differ, even though most of
subspace itself remains unchanged (see Table~2.b).
Considering again the vectors 3 to 16,
their largest amplitudes assign
the vectors: (3--6,8,14--16) to the
Ierapetra; (6--16) to the Rhodes Gyre;
(4,5,8) to wave patterns along the northern escarpment
and Mid-Mediterranean Ridge;
(3--4,7--10,12) to the Mersa Matruh Gyre;
(8--12,14--16) to the Mid-Mediterranean Jet;
(10--14,16) to the West Cyprus Gyre;
(10,12,16) to the Asia Minor Current;
and (8--16) to the
Mersa Matruh-Shikmona Gyre complex.
%lack of data: 4,5,11
Overall, the dominant 50, 100 and 200 vectors explain 41.1\%,
66.2\% and 89.4\% of the error variance
explained by the 322 vectors, respectively.
These percentages are close to these of Case 2:
the result is that here nonlinear terms
impact the dominant vectors of the subspace
but don't change much of the overall spectrum.

{\capsfn Case 4}.~The error subspaces of Cases 3 and 4 were
estimated to be similar at 86.8\% (Table~2.b). 
Up to an arbitrary sign, their dominant vectors~1 to 15
are almost identical (the vectors for Case~4
are only slightly less noisy and are thus not shown). In particular,
the two dominant vectors (LAL00) confirm that the Ierapetra
dominated the nonlinear mesoscale variability during
the Spring of 1995.
The vectors of Cases 3 and 4 only start to differ somewhat
by vector 16, mainly in the small
mesoscales. The additional vectors 323 to 400
also include 3D patterns of relatively large horizontal 
scales and usually ``wiggly'' vertical structures, or vice-versa.
This agrees with the characteristics of the
tracer vectors 323 to 400 used
in the adjustment integrations and indicates
that the ordering of scales in the Levantine sea is complex.
The decay of the error spectrum is similar to that of Case 3:
the 50, 100 and 200 dominant vectors explain
38.4\%, 62.7\% and  85.2\% of the error
variance explained by the 400 vectors, respectively.

%The other vectors relate to the main features,
%for example, the 
%Ierapetra (4--6,8,14--16), (4 is north-south motion of Ierapetra)
%Mersa Matruh Gyre (3--4,7--9,12),
%R wave: 4,7,8
%lack of data: 4,5,11
%Rhodes Gyre (6,9--16),
%Mid-Mediterranean Jet (8--12,14--15)
%West Cyprus Gyre (10--14)
%Asia Minor Current (10,12,16)
%Mersa Matruh-Shikmona Gyre complex (8--15)

{\capsfn Dynamical properties}.~From the dominant
mesoscale vectors and fields obtained, the types and scales
of possible linear wave patterns along the northern escarpment
and Mid-Mediterranean Ridge (e.g.~Case 2), and 
the dynamical balance of the Ierapetra (e.g.~Cases 3--4)
can be briefly investigated. To do so,
some scaling estimates are necessary.
Note that the values obtained from Cases~3 and~4
are similar, indicating a robustness of the results.

Using the definitions of App.~B and Fig.~16, 
locally for a wave packet along the escarpment,
$\rho_0=1028$~kg/m$^3$, $L=25$~km,
$D$ is about 1000~m, $H$ is on Fig.~4c,
$f_0=8.45 \; 10^{-5}\,$\secm1
and $\beta_0=1.86 \; 10^{-11}\,$\mm1secm1 (for $\theta_0=35.5$).
Computing a reference $\rho_s(z)$ by horizontal average
of the density field, this leads 
an average $N_s$ within $D$ of
about $2.5 \; 10^{-3}\,$\secm1, a local
$R_D = N_s D / f_0 \simeq 30$~km
and $R_{\rm ext} \simeq 1200$~km.
For $U$, based on the second vector (Fig.~16, Panels 2a--b),
the sum of the surface internal and external components
yields: $U\simeq ( 3.9 \times 1.5 \; 10^{-2} + 
5.6 \; 10^{7} \times 0.5 \; 10^{-2} /H/L ) \sqrt{1893}
\simeq 3\,$\cmsecm1 
($1.5 \; 10^{-2}$ and $0.5\; 10^{-2}$ are the
non-dimensional amplitudes,
$3.9$ and $5.6 \; 10^{7}$ the normalization factors,
and 1893 the eigenvalue). With contributions from
other relevant vectors, the characteristic velocity $U$ for
a wave packet is about $5\,$\cmsecm1.
Hence, the local values of the non-dimensional numbers are (see App.~B):
$Ro = 0.015$, $Bu_{\rm ext} = 2 \; 10^3$, $Bu = 1.4$,
$\beta = 0.23$, $\alpha \simeq 25$ and $s \simeq 0.6$
(excepted just east of the island
of Rhodes where $\alpha$ and $s$ reach 50 and 1, respectively).
The small $Ro$ and large $Bu_{\rm ext}$ indicate that
the rigid-lid linear PE regime is an acceptable approximation.
From $\beta$, the planetary vorticity gradient 
is relatively small and from $\alpha \gg 1$,
the topographic vorticity gradient is large: 
$\beta_0$ can be neglected compared to $\vert \nabla_h H \vert {f_0 \over D}$.
From $Bu \sim 1$, the length-scales of the
horizontal flow and internal stretching 
(stratification's strength) match. Since $s$ is close to one,
quasi-geostrophy may not apply. Finally, for the model topography (Fig.~4c),
$s$ and $\alpha$ are 
relatively uniform along the escarpment waveguide (about 60~km wide).
Based on these numbers and interpretations, the
topographic wave packets are of the 
coastal-trapped type (Huthnance, 1992; 
Pedlosky, 1987; Brink, 1991).
Note that here, wave patterns arise
during rigid-lid momentum adjustments (fixed
tracer disturbances): for the barotropic
mode along the escarpment,
to first-order in $s$, the frequency is then
of the order ${\vert \nabla_h H \vert f_0/D} \over  k$,
where $k = 1/L $ is the wave number.
This leads a period of about 1.45~days, i.e.~of the order   
of 1 to 2 days.
%Using Pedlosky (1987, section 6.15), the non-dimensional
%potential-vorticity equation has the general solution 
%$\psi_{\rm QG}= (-1)^j {\rm sin} [(j+{1\over 2} ) \pi z] 
%{\rm e}^{ i ( k x + l y - \sigma t)} $,
%with $\sigma_j = - {\beta k \over k^2+l^2 + (j+{1\over 2})^2
%\pi^2/Bu}$

In the ocean, wave packets would be generated
by complete dynamical adjustments, and also
by winds, atmospheric pressure and coastal flows.
A complex wave spectrum around
the above barotropic values can be expected. First,
the natural topography is more jagged and at several locations
steeper than the model topography, which we estimate could
shift the period by about 50\%.
Second, several wave types are possible along the
escarpment and/or Mid-Mediterranean Ridge: in addition to
a suite of coastal-trapped waves (e.g.~Leblond and Mysak, 1971)
like Kelvin waves (internal, external or double)
and continental shelf waves (barotropic or baroclinic),
inertial and Rossby waves can be generated.
For complete statements on scales, shapes and forcings,
additional observational, analytical and numerical
studies are necessary (e.g.~Huthnance, 1978;
Holland and Webster, 1994).
The dominance of such wave patterns
and associated mesoscale variations
in the vectors related to the Rhodes Gyre (see Cases 2--4)
suggests that they could be in a large part
responsible for the observed multiple centers of the Gyre
(Milliff and Robinson, 1992).
In particular,
%the mesoscale vectors indicate that
the main core of the Gyre is not stationary,
but undergoes translations and precessions
of synoptic time-scales (seen in our simulations).

%Note that in the region, the orientation of lobes (30 degrees 
%southeast--northwest for $\widehat u$ and
%60 degrees, southwest--northeast for $\widehat v$)
%is in part determined by the duration of the
%adjustments and thus the inertial terms,
%but not by the intrinsic properties of the
%{\it a priori} PE fields. If the duration is
%only one day, the locus of the extrema is still the same,
%along the northern coastline, but the orientation is ....;
%the importance of the corresponding inertial oscillations
%is addressed in ...

Using Fig.~17 and local scale analysis,
an approximate dynamical balance for
the Ierapetra can be obtained.
One has $\rho_0=1028$~kg/m$^3$, $L=20$~km,
$D=400$~m, $H \simeq 2500$~m (Fig.~4c),
$f_0=8.15 \; 10^{-5}\,$\secm1
and $\beta_0=1.89 \; 10^{-11}\,$\mm1secm1 (for $\theta_0=34.1$).
The vertical averages within $D$
of $N_s$ and of the horizontal $\Delta \rho$ are
$2.5 \; 10^{-3}\,$\secm1 and 0.15~kg/m$^3$, respectively.
Hence, $R_D =12$~km
and $R_{\rm ext}=770$~km.
From Table~3, the mesocale component of the total surface $U$ is
25~\cmsecm1. Adding the
subbasin-scale component (Sect.~4.1) yields
$U\simeq 45\,$\cmsecm1.
%$U\simeq (4.0 \times 8.0 \; 10^{-2} + 
%5.9 \; 10^{7} \times 2.0 \; 10^{-2} /H/L ) \sqrt{2864} = 18$~\cmsecm1 
Hence, for the Ierapetra, 
the relevant non-dimensional numbers are (see App.~B):
$Ro=0.28$, $Bu_{\rm ext}=1500$, $Bu=0.38$,
$\beta=0.017$, $\alpha=1.8$ and $s=0.5$.
The size of $Ro$ attests that nonlinearities
are locally important while $\beta$ and $Bu_{\rm ext}$ 
confirm that the planetary vorticity and stretching of the free-surface 
do not matter. In fact, the above numbers
suggest that the dominant horizontal 
momentum balance within the vertical extent
$D$ is the gradient-wind balance:
solving for $U$ in,
${U^2 \over L} - f_0 U + 
{\Delta \rho \, g D \over \rho_0\, L} \simeq 0$,
gives a root $U \simeq 50$~\cmsecm1 
(accounting for the uncertainty in the characteristic scales,
this root lies between 40 and 70~\cmsecm1).
In the vertical, since $Ro$ and $Bu$ are of the same order,
$W$ is large enough for the vertical convergence and horizontal divergence
to balance ($W \simeq U D/L$) and the horizontal changes in density
to match the vertical ones ($\Delta \rho \simeq \Delta \rho_s$).
The values of $\alpha$ and $s$ indicate
that the local topographic valley (Fig.~4c) is a
factor in maintaining the location of the vortex.
Finally, the viscous terms are only 
important in the bottom layers and at the
sub-mesoscales: the horizontal momentum Ekman number
$Ek_h=A_h/f_0 L^2$ is 0.92 for scales under 20~km but
decays to $5\; 10^{-4}$ at 50~km (see App.~C).

\vskip 0.3 truecm
\noindent
{\scaps 4.2.c A priori mesoscale error covariance functions}

The mesoscale covariance
function estimates of Cases 2--4 are now compared.
%(Figs.~18--19). 
The covariance
between the surface (5~m) temperature at (33.80N,27.85E) and the
other state variables is considered as an example.
The plots correspond to elements of a row of $\BB^p$
(size ${218943  \times p}$).
The position (33.80N,27.85E) is about 80~km east-southeast of
the Ierapetra.
%, on the outer edge of the
%anticyclone of the Mersa Matruh-Shikmona Gyre complex
%and in a branch of the Mid-Mediterranean Jet,
%which is locally flowing northeastward.

{\capsfn Case 2}.~Figure 18 shows that
the horizontal $T$--$T$ and $T$--$S$ Mexican hat covariances
are well represented.
Vertically (Panel b), the characteristics of the tracer covariances
are: $T$--$T$ and $T$--$S$ uniform within 
a 30~m surface mixed-layer; subsurface extremum in
$T$--$S$ near 300~m;
$T$--$T$ and $T$--$S$ positive down to about 1200~m where they
change sign; etc.
The $T$--$\psi$ covariance (Panel a) indicates
that, locally, quasi-geostrophic balance is a good
approximation.
However, the horizontal $T$--$\widehat u$ and $T$--$\widehat v$ covariances
are not exactly classic geostrophic double-lobe structures,
antisymmetric with respect
to the $x$ and $y$ axis, respectively.
At the surface (Panel a), the main reason is that inertial terms 
matter at (33.80N,27.85E).
In fact, if only the inertial and geostrophic
terms were important, after 2 days
of adjustment integration, with $2\pi/f\simeq0.89\,$days locally,
the $T$--$\widehat u$ pattern
at 5~m would be tilted clockwise by almost 45 degrees.
With the horizontal, vertical and bottom diffusions,
this tilt at 5~m is here reduced to about 30 degrees
(Panel a). Note that even though
inertial oscillations occur (e.g.~in numerical simulations,
the center of the Ierapetra at the surface rotates in small circles at
frequencies close to $f$), their effects on covariances
are not always desired\footnote{$^5$}
{\hsize=6.5truein
\eightpoint \hskip -23pt Inertial
oscillations arising in adjustments
can have undesired effects on covariances
% (here 30 degrees tilt).
They can be made insignificant by allowing the duration of the
adjustment integrations to be random within an inertial period.
A location where this would be useful is
inside the Rhodes Gyre where nonlinear terms are weak:
the tilt in the 5~m $T$--$\widehat u$ and
$T$--$\widehat v$ covariances there mainly depends on the duration
of adjustment, due to inertial rotations.
However, at the Ierapetra, where nonlinear terms
dominate, the duration of adjustment did not influence
the tilt in the velocity covariances (the tilt is desired).}.
In the vertical (Panel b), $T$--$\widehat u$
and $T$--$\widehat v$ 
have a zero-crossing within 500 to 600~m.
From the $T$--$\psi$ pattern,
the total velocity covariances near (33.80N,27.85E)
have a vertical minimum near 700~m.
%This minimum is deeper than that of the Ierapetra (Sect.~4.2.b).
Below the zero-crossings, $T$--$\widehat u$ and
$T$--$\widehat v$ are opposite to their
surface structures, but have substantially smaller amplitudes,
with maxima near 10\% of the surface values.


{\capsfn Cases 3 and 4}.~By construction,
the tracer auto-covariance functions for Cases~3 and~2
(Fig.~18) are equal.
Those of Case 4 are very similar (see LAL00). The main difference is a
relative volume-averaged increase in variance
of $100 \times (72-67.2)/67.2=7.1\%$ (Sect.~4.2),
which is distributed in the 3D $T$--$T$ and $T$--$S$ fields
according to the shape of the tracer vectors 323 to 400.
Figure 19 thus only illustrates
cross-covariances with velocity fields.
Comparing Cases 2 (Fig.~18) and 3 (Fig.~19a), overall,
the nonlinear adjustment of velocities
modifies the structures of the
tracer-velocity cross-covariances, but
amplitudes remain similar.
Focusing on differences, the horizontal $T$--$\widehat u$ 
structure has a smaller tilt, of about 15-to-20 degrees clockwise:
the nonlinear terms dissipate some inertial oscillations.
Cross-covariances now extend to the Ierapetra
and their gradients are tighter due to nonlinearities.
The $T$--$\psi$ structure has become asymmetric,
with cyclonic $T$--$\psi$ correlations at
% the location of
the Ierapetra.
In the vertical, structures in the upper-layers
extend deeper than in the linear adjustment case.
The other vertical changes are mainly due to the
lesser 15-to-20 degree rotation in the horizontal.
Considering the velocity cross-covariances
of Cases 3 (Fig.~19a) and~4 (Fig.~19b),
shape and amplitude
differences are small. This agrees with the
corresponding global similarity coefficient (Table 2.b).

\vskip 0.5 truecm
\noindent
{\bf 4.3 Mesoscale: {\it a posteriori} PE fields}
\vskip 0.3 truecm

At the mesoscale stage ($w=2$), for all four cases (Table 1),
the mapping procedures are as overviewed in Sect.~4.1.
The differences lie in the
{\it a priori} states (the subbasin-scale estimates),
data set and mesoscale covariance parameters (Sect.~4.2).
To exemplify the sensitivity to subspace parameters,
the focus is on Cases 2 and 3 for the effects
of nonlinear momentum terms, and on Cases 3 and 4
for the effects of the size of the subspace.
%The PE fields of Cases 1 and 4 are compared in more
%details in LAL00.

Figures 20--22 illustrate the results, variable by
variable. It is useful to first summarize some
physical reasons for the differences between the fields
before (Figs.~11 and 13--14) and after March 15
(Fig.~1 and Figs.~20--22).
For the region of common data (Figs.~4 and 15),
variations are mainly due to:
the flattening of the temperature and salinity front
along the northern coastline
which provides part of the support for
the Asia Minor Current; 
the contraction of the cyclonic Rhodes Gyre
towards a rounder shape
which is especially seen at 105 and 225~m depths; 
the weakening of the cyclonic West Cyprus Gyre;
the stretching with some surface weakening of the 
Ierapetra eddy; the mesoscale meanderings of the
sampled portions of the Mersa Matruh Gyre;
the deep anticyclonic convections in the 
Mersa Matruh-Shikmona Gyre complex;
and the formation of a surface pycnocline
and thus surface mixed-layer. By late
March, the depth of this overall warming 
surface mixed-layer varies from 5 to
10~m at the top of the two
subbasin-scale cyclonic gyres,
%(RG, WCG)
to about 50 to 100~m around the anticyclonic gyres and rings.
%(IE, MMG, eddies and rings of the MM-S-G complex).
For more on the physical processes
during the spring of 1995, we refer to LER97.

{\capsfn Tracer fields}.~Figure 20 illustrates the $T,S$ components of
$\xa$ obtained by multivariate 3D analysis (App.~A.3)
in the Case 2. By construction, the tracer
fields for Case 3 (not shown) are identical to
these of Case 2 (up to machine precision).
Comparing Cases~2 (Fig.~20) and~4
(Fig.~10 in LAL00), their two tracer fields
have small differences,
in accord with the 4.8\% increase
in the 3D tracer variance explained (from 67.2\% to 72\%).
Differences are largest mainly in the surface layers and
in the small mesoscales, which is logically the main
variability spanned by the additional
vectors 323 to 400 used in Case~4 (Sect.~4.2.b).
The results of the
comparison with the benchmark Case~1 (see Figs.~9--10 in LAL00)
are analogous to these of Sect.~3. The large mesoscales
to the subbasin-scales have small differences, within data error bounds.
The sub-mesoscales to small mesoscales present larger
differences, especially in the surface layers. This is because
many of these small scales are environmental noise
for the available data resolution (Fig.~15a)
and because the 400 eigenvectors only
explain 72\% of the total tracer variance (Sect.~4.2).

{\capsfn Transport fields}.~Figure 21 compares
the barotropic transport stream functions.
Panel~(a) is the result of Case~1
which computes $\psi$ by vertically averaging the
total velocities in thermal-wind balance with
the horizontal $T$ and $S$ analyses (see Fig.~9 in LAL00),
assuming a level of no motion at 600~m.
Panel~(b) is the $\psi$ component
of $\xa$ for Case 2: the transport is then globally estimated
from the $T$ and $S$ data (Fig.~15) using App.~A.3.
Panel~(c) is $\psi$ for Case 3, Panel (d) $\psi$ for Case 4.
The advantages of Cases 2--4 (Fig.~21b--d) are that
their $\psi$ estimates, and corresponding uncertainties (Sect.~4.4),
are in accord with the dominant SVD of the PE adjusted perturbations.
As described in LAL00, most of the regional dynamical features
can be distinguished in these $\psi$ fields.
However, in Case 1, some features likely do not
have an acceptable $\psi$ estimate.
Comparing Panels~(b) and~(c),
the effects of nonlinear terms
in the momentum adjustments agree with the local
mesoscale standard deviations of Table 3.
Differences are largest for the Ierapetra and 
Mersa Matruh Gyre.
Comparing Panels~(c) and~(d), increasing the
size of the subspace leads to small 2 to 10\%
variations in $\psi$, in accord with the global
covariance similarity of 86.8\% (Table~2.b). The largest changes
are at the Ierapetra, above the Rhodes basin and along the 
Asia Minor Current and Mid-Mediterranean Jet.

%An advantage of the present methodology is that both $\psi$
%and its uncertainty are estimated,
%in accord with the dominant covariances of the
%{\it a priori} missing PE variability.
%For example, by nonlinear and stochastic effects
%in the adjustment integrations, parts
%of this PE variability can be uncorrelated to the observations
%(here, Fig.~15) and thus logically remain pure
%uncertainty {\it a posteriori}.

{\capsfn Internal velocity fields}.~Figure 22
compares the surface internal velocities.
Panel (a) is the result of Case~1,
which computes $\widehat u, \widehat v$ in
thermal-wind balance with the horizontal
$T$ and $S$ analyses.
Panel (b) are the surface $\widehat u, \widehat v$
components of $\xa$ for Case 2:
$\widehat u, \widehat v$ are then globally estimated
from the $T$ and $S$ data (Fig.~15) using App.~A.3.
Panel (c) is the same, but for Case 3;
Panel (d) is for Case 4.
The differences between Cases~1 and~2 are small
and mainly due to differences in their
respective tracer fields. Introducing the nonlinear
terms in the momentum equations (Cases 2 to 3)
mainly affects the Ierapetra and Mersa Matruh Gyre.
The other features in the surface circulation
remain relatively unchanged. Increasing the size
of the subspace (Cases 3 to 4) leads to variations
in the small mesoscales, in accord with the
variability spanned by the additional vectors considered.

Overall, the regional increases in velocity amplitudes
from Cases 2 to 3, and from Cases 3 to 4,
logically agree with the {\it a priori} standard deviations
of Table 3 and global values obtained
from Table 2.b. The comparison of Cases 1 and 4 reveals that
their $\widehat u, \widehat v$ fields differ
mainly because of the different tracer fields (see above)
and because of ageostrophic motions in several regions
and features (see LAL00 for details).
The variations between Panel~a and Panels b--d show
that nonlinear terms (e.g.~centrifugal force)
are necessary to explain velocities
at the Ierapetra: as concluded in Sect.~4.2.b,
the gradient-wind balance is locally
a good approximation of the
momentum equations but geostrophy is not.


\vskip 0.5 truecm
\noindent
{\bf 4.4 Mesoscale: {\it a posteriori} error covariance}
\vskip 0.3 truecm

%\vskip 0.3 truecm
\noindent
{\scaps 4.4.a A posteriori mesoscale error variance}

{\capsfn Univariate 2D benchmark}.~The
{\it a posteriori} tracer error variances result
from the 2D analyses, e.g.~the 5~m errors on Fig.~15b.
These error fields give an estimate of the horizontal
mesoscale coverage of the data $\By^\o$ (236 CTDs).
The velocity error variances
% are not computed operationally,
%but they 
are set to these of velocities in geostrophic balance.

{\capsfn Multivariate 3D scheme}.~The {\it a posteriori}
error standard deviations
%(square-roots of the diagonal elements of $\BE^\a \BPi^\a \EaT$)
for Cases~2 and~3 are shown on Figs.~23 and 24
(for Case~4, see Fig.~13 of LAL00).
An important result is that error
properties are overall as one would
expect dynamically even though
only 322 or 400 vectors are employed.
The surface tracer errors of Cases~2--4
(Figs.~23a--24a) are
similar to these of Case~1 (Fig.~15b).
The differences are due to the use in Cases~2--4
of {\it a priori} $T$ and $S$ covariances
with non-zero multivariate and vertical correlations.
By construction, the tracer errors 
for Cases 2 and 3 are identical
(up to machine precision). 

In general, the structures of error
standard deviations of total velocities (Figs.~23--24)
differ from these of tracers which mimic the
patterns of hydrographic ship-tracks.
Differences in shape are much larger for
$\widehat u, \widehat v$ errors (under geostrophic influence)
than for $\psi$ errors (under hydrostatic influence).
The relationships among the velocity error patterns,
regional dynamical features, and
hydrographic horizontal and vertical resolutions
(un-sampled regions, distances between ship-tracks,
spacings between profiles along ship-tracks,
non-uniform profile lengths, etc) 
are discussed in detail in Sect.~5.b.ii of LAL00.
These discussions apply here. The
velocity errors for Cases 2--4
differ from these of Case~1 (not shown) because, in Cases 2--4,
{\it a priori} error covariances involving velocities
are non-homogeneous and anisotropic in 3D
and because the dominant momentum balance is not everywhere
geostrophy (diffusion and nonlinearities matter locally).
In particular, comparing Cases 2 and 3 (Figs.~23--24),
the nonlinear momentum terms impact
the {\it a posteriori} $\psi$ errors and
$\widehat u, \widehat v$ errors at various depths and locations.
Differences are largest
near the Ierapetra, Mersa Matruh Gyre and
Mersa Matruh-Shikmona Gyre complex which are the
features found to be the most nonlinear {\it a priori} (Sect.~4.2). 
Error patterns in Case~3 are broader and
more uniform (less depend on data location) than in the linear 
Case~2 because of nonlinear mixing.
The transport errors (Figs.~23b--24b) have similar shapes
but amplitudes in Case~3 are on average larger by about 4 to 6\%
(in accord with the {\it a priori}
variance increase of 5.4\%, Sect.~4.2.a), with peak
increases at about 10 to 20\% for the Ierapetra.
The effects of nonlinearities are larger on
internal velocity errors. At 5~m (Figs.~23a--24a),
the $\widehat u, \widehat v$ errors
are about 30 to 50\% greater at the
Ierapetra, and about 10\% greater in the Kasos Strait,
around the Rhodes Gyre and in regions of low data
resolution east of 28E.
At 410~m and 1500~m, $\widehat u, \widehat v$ errors
remain of broader extent and larger (e.g.~see the Ierapetra,
anticyclone in the Mersa Matruh-Shikmona Gyre complex and
around the Rhodes Gyre),
but the increase (around 5--10\%) is usually smaller
than at the surface.
Comparing Cases~3 (Fig.~24) and~4 (Fig.~13 in LAL00),
the error amplitudes of Case~3
are slightly smaller, by an amount that logically agrees with the
global similarity coefficient (Table 2.b)
and smaller regional {\it a priori} error estimates (Table~3).
Patterns for the two cases are also very close,
excepted at the small mesoscales because of the additional
space spanned by the vectors 323-to-400.

\vskip 0.3 truecm
\noindent
{\scaps 4.4.b A posteriori mesoscale error eigenvectors and spectrum}

The first and second non-dimensional
{\it a posteriori} error singular vectors for Cases~2 and~3
are illustrated on Figs.~25 and 26, respectively
(for Case 4, see Fig.~14 of LAL00).
In the linear Case 2, the first vector accounts for
5.7\% of the normalized {\it a posteriori}
error variance; the second for 5.1\%.
For Case 3, these numbers are 5.2\% and 4.8\%
respectively: they are smaller because additional
kinetic energy due to the nonlinear terms
remains un-measured in the
{\it a posteriori} eigenvalue spectrum.
For Case 4, values are 4.25\% and 4.05\%,
which are smaller than for Case 3
due to the larger subspace size.
Overall, the eigenvalue spectra differ
for the same reasons, but by much smaller amplitudes.
In Case~2, the dominant 50
vectors explain 75.8\% of the
variance explained by all 322 vectors; the dominant 100
explain 87.2\%. In Case~3, these numbers remain
nearly unchanged, 75.5\% and 87.2\%.
In Case~4, they are logically slightly smaller:
71.5\% for 50 vectors
and 84.5\% for 100 vectors out of all 400 vectors.
For all Cases, in comparison to {\it a priori}
values (Sect.~4.2.b), the analysis has here
steepened the top and flattened the middle and end of the
error spectrum.

Comparing {\it a posteriori} vectors one by one,
the three Cases have almost the same first 
vector (Panels 1 on Figs.~25--26).
However, due to the nonlinear adjustments, other vectors
are reordered in Cases~3--4: e.g.~up to
the effects of the orthogonality constraint and
arbitrary sign, the second vector of Case 2 is as the third
of Case~3 (not shown) and  inversely
%the third of Case~2 as the second of Case~3 
(Fig.~26). Finally, the dominant vectors of Cases~3
and 4 (Fig.~14 of LAL00)
are almost identical, indicating a
convergence of the subspace.

For each of Cases 2--4,
the data $\By^\o$ have redistributed the variance
in the subspace. Globally,
the dominant {\it a posteriori} vectors have little
in common with the {\it a priori} ones.
In Case 2 (Fig.~25), the two dominant vectors
are associated with the low data resolution along 30.75E (Fig.~15a)
and lack of data in the northeast
corner (north of 35.5N, east of 29.5E).
%The zonal wavelength of the second vector in that corner are half
%of that of the first vector due to the orthogonality constraint.
In Cases~3 (Fig.~26) and 4 (Fig.~14 of LAL00),
this holds for the first vector, but not
for the second. Because of the nonlinear adjustments,
the second vector corresponds
to the low data resolution patches centered on 29E,
which is across the Mersa Matruh-Shikmona Gyre complex
where nonlinearities matter (Sect.~4.2).
%
%The regional properties of the data coverage and dynamics strongly
%influence the {\it a posteriori} vectors. A detailed
%analysis of these influences is given
%for Case 4 in LAL00;
%a similar analysis applies for the Cases plotted here.
%
Locally, the properties of the dynamics and data coverage also
affect these vectors. A study
of regional dynamical influences is given for Case 4 in LAL00;
similar statements apply to the Cases shown here.
Focusing on the effects of data coverage, 
considering for example the first vector
of Cases 2--4 (e.g.~Figs.~25--26, Panels 1a--b),
for each variable, the amplitude of
the pattern south of 35.5N (northern 
ship-transect) is at least twice as small
as it is north of 35.5N (no data in this corner). 
The zero-crossings (maxima) of the
$T$ and $S$ patterns of this first vector
are also aligned with the 
(lack of) meridional ship-transects (Fig.~15):
from 34N to 35.5N included, zero-crossings
are every 0.5 degrees of latitude, while north of 35.5N,
the last lobe covers 0.8 degrees of latitude, 
across the no-data corner. 
%These properties have nothing to do with physics.
The profile depths also influence error patterns.
For example, east of 30.25E, all profiles are about 1000~m deep,
except for two at (31E,35.5N) and (31E,34.5N)
which are 2000~m deep. Logically, these two deep profiles have
reduced the error on $\widehat u$ at 35N
but not on $\widehat v$ (see Panels~1b from 1000~m to 2000~m).
Analogous facts hold for other vectors.
Other lessons that can be learned from the analysis of
{\it a posteriori} error covariances, including
guidance on the design of sampling strategies,
are discussed in LAL00.

\vskip 0.3 truecm
\noindent
{\scaps 4.4.c A posteriori mesoscale error covariance functions}
 
The comparison of Figs.~18 and 19 (Sect.~4.2.c) with the
{\it a posteriori} error covariance functions plotted on 
Figs.~27 and 28 clearly shows that
error covariance functions are also modified by
data. All amplitudes are reduced in accord with the data
coverage and measurement error model (App.~A and LER97). 
At the position (33.80N,27.85E),
both the horizontal and vertical correlation scales
have been reduced to small mesoscale values.

{\capsfn Case 2}.~Keeping data coverage (Fig.~15a) in mind,
the horizontal structures (Fig.~27a) clearly show the
influences of the low data resolution south of 34N
within 28--29E and of the three ship-tracks along
26E, 26.5E and 27E. In the vertical, the
$T$--$T$ and $T$--$S$ structures have stronger negative
lobes and are more surface intensified
because {\it a priori} $T$ and $S$ data errors
are modeled to decay with depths (LER97 and Fig.~31
in Sect.~5.2). The $T$--$\widehat u$ and
$T$--$\widehat v$  structures have however
remained similar and surface intensified,
except that the
negative lobe of the $T$--$\widehat u$ structure
is now to the east, in accord with thermal-wind effects
and the lower data resolution in that direction.

{\capsfn Cases 3 and 4}.~Only Case 4 is shown (Fig.~28).
The covariances for Case~3 are very similar,
except that amplitudes are smaller
by an amount in accord with the global
similarity coefficient of Table 2.b
(these amplitude changes can also be seen by comparing
Figs.~27 and~28).
In the horizontal (Fig.~28a), the structures of Case~4
are similar to these of Case~2 (Fig.~27a),
except that gradients are tighter
due to nonlinear and larger subspace size effects.
In the cross-covariances with velocity fields,
differences are largest where features have been shown
to be most nonlinear (e.g.~Table~3).
In the vertical, gradients are also tighter
for all covariances.
The $T$--$\widehat u$ and
$T$--$\widehat v$ structures have remained
similar, but the $T$--$T$ and $T$--$S$ structures
are more complex, in accord with the
additional small mesoscales explained by the
vectors 323--400.

\vfill\eject
%\vskip 0.5 truecm
\noindent
{\scaps 5. QUANTITATIVE EVALUATION OF FIELDS AND IMPACT
ON FORECAST PERFORMANCE}
\vskip 0.3 truecm

To complement the sensitivity studies
(Sects.~3--4), the two mapping schemes (Cases~1--4) and
error subspaces of Cases~2--4 are now
evaluated by measuring the quality of forecasts initialized
from the four different analyses.
Starting from each gridded state (e.g.~Figs.~20--22),
the numerical PE model (App.~C) is first
integrated forward in time.
The forecasts are then compared
to {\it in situ} hydrographic data
and satellite sea surface temperature (SST) data.
Due to the interactions between tracers and velocities, 
the effects of the initial flow fields 
and velocity error subspaces are also assessed.

Even though the {\it in situ} data are local and
intermittent, during March 26--April 8,
samplings are relatively intensive (see
Fig.~31a below). The corresponding CTD profiles are
used for the forecasts' evaluations.
Presently, daily satellite images are often cloudy and can contain
unphysical, non-uniform variations, sometimes
of more than 1$^\circ{\rm C}$. Filtering and 
time-averaging can smooth out this
noise, but can also smooth out some oceanic features
or patch together sparse daily snapshots (which can
occur when a storm passes over the domain in a few days).
The high-resolution images provided
by GISIS (1995) are here interpolated onto the
numerical model grid, masking cloudy areas.
During the week of March~27--April~3, on three days,
different halves of the domain are clear and 
contain limited noise (e.g.~Fig.~29a). This leads
to a relatively good coverage
for the week-composite (Fig.~30a), even though one
should beware of non-synoptic patching. These interpolated day and week
images are used for the forecasts' evaluations.

An issue is that the interval during which forecasts
can be compared is limited by the predictive
capability of the prediction system.
Beyond this limit, nonlinear integrations of data and model errors
could reverse conclusions
(i.e.~indicate that the worse mapping is best). For the area,
period, regime of interest and data available,
the present limit is near 10 days (LER97).
Forecasts are thus not compared beyond that.
Remaining issues are to define the error and choose its
measure. Since comparing forecasts to analyses
of data can be misleading,
data and forecasts are here directly compared at data points.
These data-forecast residuals
%at data points
can be a good estimate of the forecast error, in so far as
they are larger than the estimate of the errors in the data.
They are measured
by root-mean-squares (RMS), correlation coefficients
and pattern correlation coefficients (PCC).

\vskip 0.5 truecm
\noindent
{\bf 5.1 Satellite SST data}
\vskip 0.3 truecm

Considering first daily images (e.g.~Fig.~29),
it is challenging to determine qualitatively which of the
forecasts is closest to the satellite SST, even though
forecasts differ. Quantitatively, their global scores
are indeed similar (see Table 4 below).
Focusing on differences, overall, gradients 
are more accurate in the multivariate 3D scheme (Cases~2--4)
than in the benchmark (Case 1). The branch
of the Mid-Mediterranean Jet near 33N (on Fig.~29a, seen starting at 29E)
is also better represented in Cases~2--4 (Figs.~29c,d,e).
The lobes of the Mersa Matruh Gyre are visible in Cases~1 and~4 
(Figs.~29b,e), but not so much in the Cases with
the smaller subspace size (Figs.~29c,d).
For the Ierapetra and its surrounding cold vortices,
Cases~1 and~4 have similar large mesoscales.
At the small mesoscales, Cases 1 and 4 are better than 
Cases 2 and 3; overall, the Ierapetra gradients of Case 4 are
the most accurate. Comparable comments hold for the
other daily images (not shown).

Considering the week composite (Fig.~30), similar statements can
again be made, including for the Rhodes Gyre which is now
visible on the SST. Of course, time-scales longer than a week
are smoothed out in the forecast averages (Figs.~30b--e).
The sub-mesoscales and localized
steep gradients (due to patching) 
that remain on the week SST image are noise,
which should not be considered in this week-average comparison.

In Table 4, the skill of forecasts generated from the different analyses
is quantitatively evaluated. In all computations,
only cloud-free points are considered. For the three relatively
clear days (Table~4a), the correlation coefficient between the
satellite SST day composite and the estimated temperature
field at 5~m for that day is given. To obtain these
coefficients, the horizontal mean of each field is removed first.
The week composite
%\bigskip

\eject}

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf

\centerline{Table 4. Day and week measures of skill for temperature:}
\centerline{satellite SST vs. estimated T fields at 5~m}
\vskip -0.5truecm
$$\vbox{\halign{
#\hfil&~#\hfil&~~#\hfil&~~#\hfil&~~#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
a.~{\capsl Day Cor.~Coef.} & Case & March 27 & March 28 & April 3\cr
~ Day composite vs.  & & Nowcast & 1-day forecast & 7-day forecast\cr
\noalign{\vskip2pt}
& ~ 1 & ~~ 0.57 & ~~ 0.76 & ~~ 0.82\cr
& ~ 2 & ~~ 0.67 & ~~ 0.77 & ~~ 0.81\cr
& ~ 3 & ~~ 0.67 & ~~ 0.76 & ~~ 0.81\cr
& ~ 4 & ~~ 0.66 & ~~ 0.76 & ~~ 0.81\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
b.~{\capsl Week measures} & Case &  Cor. Coef. & ~~ PCC & RMS Diff.\cr
\noalign{\vskip2pt}
~ (March 27--April 3 & ~ 1 & ~~ 0.76 & ~~ 0.18 & ~ 0.3912\cr
~ week composite vs. & ~ 2 & ~~ 0.75 & ~~ 0.20 & ~ 0.3894\cr
~ average of eight & ~ 3 & ~~ 0.75 & ~~ 0.21 &  ~ 0.3870\cr
~ estimated T fields)& ~ 4 & ~~ 0.76 & ~~ 0.20 & ~ 0.3895\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
\vskip -0.5truecm
{\doublecolumns
\noindent
(Table 4b) is then 
compared to the week average of the day-to-day
temperature forecasts at 5~m.
To take advantage of the larger coverage (Fig.~30a),
more measures are evaluated.
In addition to the correlation coefficient, a PCC and 
the RMS of the residuals between the clear satellite data
and averaged forecasts are also calculated.
The PCC is defined by, 
$
{\rm PCC} \doteq {
(\T^S-\T^{b})^T\, ({\widehat \T}-\T^{b})
\over ||\T^S-\T^{b}||_2\; ||{\widehat \T}-\T^{b}||_2}
$,
where the vector $\T^S$ is the satellite week composite,
${\widehat \T}$ the estimated week average and
$\T^{b}$ the background set to the
satellite month composite for March 1995
(only common clear points are considered).

%Initial location of Table 4
%\noindent
For the day snapshots (Table~4a), the main difference
is for the nowcast on March 27. The estimates of Cases~2--4
are about 15\% better than that of the univariate 2D Case~1.
On other days, skill measures are very similar.
A surprising result is that there is very little
difference among Cases~2--4. This is 
mainly due to the uncertainty in the satellite data,
even though it also reflects the subspace convergence.
Another effect of data uncertainties is the
increasing correlation with time.
For the week composite (Table~4b), correlation
coefficients are very close for all cases.
The PCC numbers indicate that Case~1
is the least effective, while Case~3 is best
%(even though differences among Cases~2--4 are insignificant,
(this shows that adding more vectors
can at times decrease performance).
These PCC numbers are relatively low. This is because they focus
on the small mesoscale (the satellite month composite is removed)
for which the skill is smaller than for the large mesoscale
and subbasin-scale due to the initial data resolution,
and also because skin effects are in the week and month composites
but not in the average ${\widehat \T}$ fields.
The RMS numbers indicate a RMS error slightly 
larger for Case 1 than for Cases 2--4,
but differences are not very significant
(the estimated error in the {\it in situ}
data at the surface is around $0.3^\circ$C, see Fig.~31 below).
For the day images,
the RMS and pattern correlation coefficient (not given)
lead to the same conclusions.
Overall, the satellite data
indicate that the skills of
Cases 2--4 are similar and only 
a few times slightly larger than the skill of Case 1.

\vskip 0.5 truecm
\noindent
{\bf 5.2 {\it In situ} data}
\vskip 0.3 truecm

Considering the {\it in situ} $T$
and $S$ data (Fig.~31a), time is divided in periods
of three days. For the data collected in each
of these periods, data-forecast residuals are computed
at the intersections of the CTDs with the model grid, using
the forecast for the center of the period (details
on the measurement models, see App.~A, are in LER97 and LAL00).
These residuals are measured by RMS, averaging level by level.
On Fig.~31b, the resulting RMS-error (RMSE) estimates
are plotted for the temperature and salinity of Cases~1 and~4,
on 12 horizontal model levels (from 0 to 500~m depth).
Overall, forecast errors increase with forecast durations
and decrease with depths in proportion to the amplitude
of the variability.
For March 26--28 (Rhodes Gyre area) and
March 29--April~1 (near 30E),
most RMSEs are below the estimated error
standard deviations of the data. Both schemes
likely issue good forecasts. At certain depths,
Case~1 is better than Case~4, and vice versa.
Focusing where RMSEs are significant, Case~4 is better
than Case~1 at the surface, while the opposite holds near 500~m. 
For April~3--5 (data along the Cretan Arc), RMSEs are much larger,
in accord with the forecast duration but also with the lesser quality of the
numerical model in this region. For temperature, Case~1 is better
than Case~4 at all depths. For salinity, Case~4 is substantially better than
Case~1 in the surface layers, but is worse below.
For April~6--8 (period with most data),
in the surface layers, Case~4 is better for
temperature, but not for salinity. At mid-depths,
the situation reverses. Near 500~m, Case 1 is best overall
but errors become smaller than the estimated data error.

These results first confirm that the two schemes have
similar skill. They also restate
that the univariate 2D scheme can be too
sensitive to small scale noise (e.g.~Sect.~3.2), hence the
usually lesser skill in the surface layers. They show
that the present multivariate 3D scheme can sometimes be less
effective at depth. To improve this, either the
subspace can be enlarged or the statistical
specification of the initial covariance decomposition
(Apps.~A.1--2) refined.
For example, in coastal regions, the lack of
better performance at depth is
likely due to the tracer covariances which were simplified
(e.g.~horizontally homogeneous, vertical
structures dominated by open-ocean data instead of coastal ones)
for fair comparisons with the benchmark (see Sect.~4.2).

%East of the RG is 
%an open-ocean location where the ageostrophic components
%are not negligible since the horizontal density gradients
%in this area are weaker at all depths than at other locations
%around the RG. By conservation of mass, the RG is closed,
%in part by cyclostrophic/advection balance.

%\vfill\eject
\vskip 0.5 truecm
\noindent
{\scaps 6. SUMMARY AND CONCLUSIONS}
\vskip 0.3 truecm

This foregoing study investigated the influence of {\it a priori} parameters
on the error and variability subspace estimation and
corresponding mapping methodology described in
(Lermusiaux {\it et al.}, 2000, henceforth, LAL00).
The sensitivities of the subspace and of the
{\it a posteriori} gridded fields and errors
to the size of the subspace, scales considered
and nonlinearities in the dynamical adjustments
were exemplified and studied for 
mesoscale to subbasin-scale physical fields 
in the northwestern Levantine during
February~10-March~15 and March~19-April~16, 1995.
The {\it a posteriori} temperature, salinity and velocity
fields and errors obtained
were evaluated by comparison among each other
and by using the univariate 2D analysis scheme of HOPS
as a benchmark. Forecasts generated from the various
analyses were
qualitatively and quantitatively compared to
{\it in situ} and satellite data.
A simple back-of-the-envelope rule for the
required size of error subspaces was also outlined
and utilized.

Several three-dimensional, multivariate and multiscale
properties of the subspace estimation algorithm and analysis
scheme were illustrated and discussed. The global mapping
was shown to: estimate simultaneously all physical fields and errors
from the tracer data; usefully decompose multivariate error covariances;
lead to velocity fields
that are close to being PE adjusted and
{\it a posteriori} velocity errors that are PE consistent;
improve vertical correlations and the
filtering of the environmental noise;
and, yield more accurate field gradients.
Overall, in comparison to the geostrophic benchmark, the largest effects
were found where topography, nonlinear and diffusion terms mattered.
In particular, {\it a posteriori} errors were
found to be complex, reflecting
the properties of both the available data sets and
regional dynamics (not always close to geostrophic).
Based on the forecasts from the various analyses,
the multivariate 3D scheme performed as well
as the benchmark. Similar conclusions have been
drawn in meteorology (e.g.~Andersson {\it et al.}, 1998).
Presently, this is likely in part because of the limited data
and because the 3D scheme
was simplified for fair comparisons.
%in general, all of its capabilities,
%e.g.~non-homogeneous horizontal tracer covariances and multiscale
%interactions, should be utilized.

From the studies of the sensitivity to scales, it
was confirmed that for a fixed domain size,
smaller scales usually require larger subspaces. This is
valid especially when subspaces are specified based only on 
common statistical models.
When a dynamical model is used or when the subspace itself is dynamical,
the ordering of scales becomes more complex
and the size requirements more variable.
From the studies of the sensitivity to the size of the subspace,
it was exemplified that
truncating the {\it a priori} error covariance
to a converged subspace is efficient and that
{\it a priori} error eigenvectors of eigenvalues
negligible with respect to data errors can be neglected (App.~A.4).
In the Levantine example,
increasing the subspace size usually improved
the analysis at depths and the small mesoscales.
Overall, the variations of {\it a posteriori} fields
with the parameters of the error subspace
were as indicated by a global similarity coefficient (Table~2):
if {\it a priori} error covariance matrices differed by about
5\% (in the standard deviation sense), so did the fields and 
{\it a posteriori} errors.
From the studies of the sensitivity to the linear-nonlinear
dynamics of the adjustments, it was found however
that even when subspaces are globally similar,
some of their dominant vectors can be quite different.
The nonlinear momentum correction modified
the dominant vectors from large-mesoscale variations
of the Rhodes Gyre to mesoscale variations of the Ierapetra eddy.
It also modified the shape of error spectra,
the structures of covariance functions
and some properties of the {\it a posteriori} fields and errors.
Relative linear-nonlinear differences were computed
and as anticipated found to be largest for the features and regions
with high shears, around deep eddies or gyres
(e.g.~Ierapetra, Mersa Matruh Gyre), along fronts
(e.g.~Mid-Mediterranean Jet) and in Straits.
%at the location of where cyclostrophic balance dominates.

%The computational costs of the two schemes are discussed in LAL00.
%
%The filtering of the environmental noise
%is facilitated by the use of the dominant,
%three-dimensional and multivariate, error eigenvectors.
%The increased noise filtering ability of 3D schemes
%was also recently shown in univariate analyses,
%mainly using diagonal weights, by Lu and Browning (1998).
%The result is extended here to
%multivariate analyses using non-diagonal covariances.
%
%Difference lin nonlinear is: advection terms set to zero.
%Hence, difference between geostrophy (Case 1) and linear ESSE (Case 2)
%for internal velocities is: momentum diffusion and momentum boundary layers
%(bottom, coast).

Several dynamical properties
were obtained by studying the {\it a priori} estimates of the
tracer and velocity variances, spectra of eigenvalues,
shapes and amplitudes of dominant multivariate eigenvectors,
and structures of (cross)-covariance functions.
The spatial distribution
of the characteristic hydrographic and velocity
scales and dominant dynamical balances were quantitatively
obtained for most features of the northwestern Levantine.
In particular, the Ierapetra eddy was found to be close
to gradient-wind balance in the Spring of 1995.
Coastal-trapped waves were also inferred to be likely
along the northern escarpment of the basin and
to account for a significant portion of
the global variations of momentum and density fields
in the region. Some 
interactions among features were also discussed.
% 
%analysis of vertical and 3D EOFS.
% 
%To explain a specific variability,
%a group of vectors was required.

Many physical and computational properties of the present
methodology are discussed in LAL00, including
general directions for improvements, several algorithm
variations and problem-dependent specifics.
Focusing on research directions
associated with the present sensitivity analyses,
a first issue relates to data availability and 
accurate multiscale statistical models.
Mainly for fair comparisons with the univariate 2D scheme
but also for the lack of sufficient data, the multivariate 3D scheme
was simplified: a successive correction approach was used
instead of a general one for multiple interacting scales
and {\it a priori} tracer covariance functions
were assumed homogeneous and isotropic in the horizontal
instead of being non-homogeneous and anisotropic in 3D as
velocity covariances were.
Both of these assumptions can be removed
in principle (LAL00).
In regions with more data than here (to allow accurate fits),
the performance of the 3D scheme may then be improved.
Another direction concerns algorithms for efficient
dynamical adjustment of the error subspace estimate.
For example, specific models or processes
(e.g.~free surface, coupled biochemical irreversibilities, 
stochastic internal waves) may require specific
adjustment procedures.
Additional convergence criteria (LER97) should also be investigated.
For example, our experience in other regions
suggest that taking into account each variable
separately as well as all variables as a whole (eqs.~A22)
is useful. In these criteria, ideal normalizations
for specific purposes should also be further researched.

%
%3. Related problems, open issues: ---------
%
%Here only a first-guess. But can be used in a subsequent ensemble
%to introduce anisotropy or multiscales in the
%simplified statistical model fit to data.
%
%Talk about sea surface data extension
%and the relations between C and P.
%And that P constructed here could be used as C or a combination
%of C and P in OI approaches.
%
%Add vertical extension (Fischer)
%Can be used as a measurement model.
%
%inverse methods general circulations
%Add some Wunsch 1996 comments (level of no motion).

%Having about the equivalent of about two people working
%full time on model improvements,
%you obtains, on average, 10 to 15 percent increase in forecast
%quality. This is what we get here.
%
%Biological model can be used and validated
%by creating the non-observed phytoplankton and zooplankton 
%variability from observed hydrographic and nitrate variabilities.

It is important to note that most
results obtained here apply to the analysis step
of data assimilation schemes combining
data and forecasts based on their 
respective error subspaces.
Hence, considering time evolutions,
the research opportunities include:
subspaces for multiresolution interpolation and nested grids,
efficient error initializations for real-time error forecasting
and adaptive sampling (e.g.~Lermusiaux, 1999b; Arango {\it et al.}, 2001),
model improvements based on accurate error estimates (LER97),
dynamical studies based on the
quantitative organization of multivariate variability subspaces
(e.g.~Lermusiaux, 2001; Lermusiaux and Robinson, 2001),
and complex interdisciplinary applications
(e.g.~Robinson and Lermusiaux, 2001).

%The efficient computation of accurate
%initial conditions is important.
%This is commonly accepted for the initial state, but
%it is perhaps more controversial for the initial errors.
%In fact, for operational meteorological constraints,
%initial error conditions have sometimes been
%usefully tailored towards the best error forecast.
%Such attempts include the use of singular vectors
%for maximum error growth (e.g.~Molteni {\it et al.}~1996;
%Barkmeijer {\it et al.}~1998) and the breeding method
%(Toth and Kalnay, 1997) in
%ensemble forecasts. Following LER97, the
%present aim is towards the fields and errors that are
%the closest (in the sense of the error measure chosen)
%to the true ones, at all times, hence initially.
%This is especially important in oceanography,
%where an initial error can have
%as much value as an error forecast.

\vskip 0.5truecm
\noindent
{\bf Acknowledgements}

PFJL is grateful to Professor A.R.~Robinson,
Dr.~Pat J.~Haley and Mr.~Wayne G.~Leslie
for their interest and encouragements.
Mr.~R.~Meisner (DLR) is thanked for
answering questions on the GISIS filtering
of the satellite data and Ms.~M.~Armstrong
for preparing some figures.
We thank the numerous individuals who participated in the
collection of the LIW95 data and the international
POEM collaboration sponsored by the IOC of UNESCO.
The work of PFJL and CJL was supported in part by
the Office of Naval Research under grants
N00014--95--1--0371, N00014--97--1--0239 and
N00014--97--1--1018 to Harvard University.
%6.1 N00014--95--1--0371,     
%6.2 N00014--97--1--0239 and
%LOOPS: N00014--97--1--1018
%Uncertainty: N00014--00--1--0771

\vfill\eject
%\vskip 1.truecm
\centerline{APPENDIX A}
\centerline {\bf A priori error subspace estimation}
\centerline {\bf and mapping methodology}
%Notation and Evaluation of the {\it a posteriori} ESSE state}
\vskip 0.3truecm

This Appendix summarizes algorithms introduced
in LAL00 for the estimation (initialization) of 3D
multivariate and multiscale geophysical fields and their dominant errors.
An important step is the construction of the dominant
%JAOT: The focus here is on the dominant
{\it a priori} error covariance.
Details and properties of these algorithms are in LAL00.
Wherever possible, the notation of Ide {\it et al.}~(1997) is adopted.
The translation into other standards (e.g.~Jazwinski, 1970; Gelb, 1974;
Daley, 1991; Bennett, 1992; Wunsch, 1996;
Robinson {\it et al.}, 1998) is easily carried out.

The framework is that of continuous-discrete estimation.
The true (superscript $^t$) state vector $\Bx^{\t} \in\bR^n$ is
assumed subject to
the stochastic dynamical and measurement models,
respectively,
$$
\eqalignno{
d \Bx^{\t}=& \, {\cal M}(\Bx^{\t},t)\, dt+d \Beta^{\t}(t)~,&
{\rm (A1a)}\cr
\By^\o=& \,{\cal H} \big(\Bx^{\t}(t_0)\big)+\Bepsilon~,&
{\rm (A1b)}}
$$
where ${\cal M}$ is the dynamics operator,
$\Beta^{\t}$ a random process of zero mean and covariance matrix $\BQ$,
$\By^\o\in\bR^m$ the data vector,
${\cal H}$ the measurement operator and
$\Bepsilon$ a random process of zero mean and covariance matrix $\BR$.
Specifics of (A.1.a--b) are in LER97 and LAL00.
The time $t_0$ for the estimation is fixed.
An unbiased estimate of $\Bx^{\t}(t_0)$ is denoted by $\Bx \in\bR^n$.
The state error covariance is defined
by $\BP\doteq {\cal E}\big \{  \big(\Bx-\Bx^{\t}(t_0)\big)
\, \big(\Bx-\Bx^{\t}(t_0)\big)^T  \big \}$ $\in\bR^{n\times n}$,
where ${\cal E}\left \{ \cdot\right \}$
is the statistical mean operator.
For the present minimum error variance criterion (Sect.~1),
an error subspace is the portion of $\bR^n$ of dimension $p$
which is spanned by the eigenvectors
corresponding to the dominant $p$ eigenvalues of a (normalized)
error covariance. With this rational reduction,
the optimal estimate becomes:
$$
\eqalignno{
\bigg\{
\Bx^\a \biggm |
\min_{\Bx}\,\, {\rm tr}\, [\BP^{\a^{\scriptstyle p}}]\, ,
\; \, & {\rm knowing}\; [\By^\o, \BR]\; ~~~~ & \cr
& {\rm and}\, \left[{\xb, \BP^{\b^{\scriptstyle p}} }\right] \bigg\} ,~
& ~~ {\rm (A1c)} }
$$
in which $\Bx^{\t}(t_0)$  is subject to (A1a--b).
The superscripts~($^\b$) and ($^\a$)
refer to estimates before mapping ({\it a priori})
and after mapping ({\it a posteriori}).
The superscript~($^p$) indicates a rank-$p$ approximation:
e.g.~the {\it a priori}/{\it a posteriori} principal error covariance
$\BP^{\b^{\scriptstyle p}}$/$\BP^{\a^{\scriptstyle p}}$
is the rank-$p$ eigendecomposition
of $\BP^\b$/$\BP^\a$, respectively.
What follows summarizes the: approach and algorithms
for specifying the observed portions of
error covariances; present algorithm; mapping
scheme; and, relations with the Kalman update.
Quantities marked with asterisks are normalized.

\vskip 0.3truecm
%\vfill\eject
{\bf A.1 Observed portions of the {\it a priori} error
covariance: general specification}

Following LAL00,
the value of a field $\phi$ at location $(\Br,z)$,
where $\Br=(x,y)$ and $z$ are
horizontal and vertical positions in a
suitable coordinate system, is denoted by $\phi(\Br,z)$.
The {\it a priori} error covariance function
for the fields $\phi(\Br_1,z_1)$ and $\varphi(\Br_2, z_2)$ is
then,
$$
\eqalignno{
C_{\phi \varphi}(\Br_1,\Br_2,z_1,z_2) & \doteq 
{\cal E}\bigg \{ \big(\phib(\Br_1,z_1)
-\phi^{\t} (\Br_1,z_1) \big) ~~~  &
\cr
& \big(\varphib(\Br_2,z_2) -\varphi^{\t} (\Br_2,z_2) \big)
\bigg \} .
& ~ ({\rm A2}) }
$$
An often efficient representation of (A2) is
obtained by expanding the {\it a priori} errors into,
$$
\eqalignno{
\phib(\Br,z)- \phi^{\t}(\Br,z) =&
\sum_{i=0}^{\infty} \phi_i(\Br)\,\, Z_\phii^*(z) ,
~~~ & \rm{(A3a)} \cr
\varphib(\Br,z)- \varphi^{\t}(\Br,z) =&
\sum_{i=0}^{\infty} \varphi_i(\Br)\,\, Z_\varphii^*(z) ,
~~~ & \rm{(A3b)}}
$$
where the $Z_\phii^*(z)$'s and
$Z_\varphii^*(z)$'s are vertical functions,
normalized in some suitable sense.
Substituting (A3a--b) into (A2) yields,
$$
\eqalignno{
C_{\phi \varphi}(\Br_1,\Br_2,z_1,z_2) = \sum_{i,j=0}^{\infty,\infty} &
A_{\phii,\varphij}\; C^*_{\phii,\varphij}(\Br_1,\Br_2) ~~~~ & \cr
& Z_\phii^*(z_1)\;  Z_\varphij^*(z_2) ,
& ({\rm A4a})}
$$
where the factor
$A_{\phii,\varphij} = A_{\varphij,\phii}$ is the horizontal average
of ${\cal E} \{\phi_i(\Br) \varphi_j(\Br)\}$ and
$C^*_{\phii,\varphij}(\Br_1,\Br_2)
\doteq {\cal E} \{\phi_i(\Br_1) \varphi_j(\Br_2)\}^*$
the non-dimensional horizontal cross-covariance function
associated with $\phii$ and $\varphij$.
%in general non-homogeneous and anisotropic.
The representation (A4a) is efficient if,
$$
\eqalignno{
C_{\phi \varphi}(\Br_1,\Br_2,z_1,z_2) \simeq \sum_{i,j=0}^{I,J} &
A_{\phii,\varphij}\; C^*_{\phii,\varphij}(\Br_1,\Br_2) ~~~~ & \cr
& Z_\phii^*(z_1)\;  Z_\varphij^*(z_2),
&({\rm A4b})}
$$
with small indices $I$ and $J$
(i.e.~$A_{\phii,\varphij}$'s decay rapidly
with increasing $i$ and $j$).
\vfill\eject}
%\bigskip

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf
\centerline{Table A.
{\it A priori} error subspace construction for one scale ($w$ is omitted)}
\vskip -0.7truecm
$$\vbox{\halign{
#\hfil&\quad#\hfil&~#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip5pt}
{\capsfb a. Observed Portions}\cr
~Tracer Vertical EOFs\cr
\hfill Compute scale-restricted {\it a priori} tracer residuals:&
$
\Bd_i= \By^\o_i - {\cal H}_i(\xb)~, \quad i=1, \ldots ,s$~.
&(A7)\cr
\hfill
Remove horizontal average, normalize:&
$
\Bd_i^* = \BN_d^{-1}
(\Bd_i - \overline{\Bd_i})~.$
&(A8)\cr
\hfill
Compute SVD of normalized tracer residuals:&
$
{\rm SVD} \left({\BS_d^* \doteq
\big[{ \Bd_1^* ,\,
\dots , \, \Bd_s^* }\big]}\right)=
\BE_d^* \BSigma_d^* \BV_d^{*^{\scriptstyle T}}$.
&(A9)\cr
\hfill
Leads to tracer vertical covariance decomposition:&
$
\BC^z_{\trc} = \BE^z_{\trc} \BPi^z_{\trc}
\BE^{z^{\scriptstyle T}}_{\trc}\,$,
&(A10)\cr
&\qquad where $\BE^z_{\trc} = \BN_d \BE_d^*$ and 
$\BPi^z_{\trc} = \bSigma_d^{*^{\scriptstyle 2}}/ s \,$.
\cr
\noalign{\vskip4pt}
%\cr
~Tracer Horizontal Covariance Eigendecomposition\cr
\hfill Specify normalized covariance from analytical fit:&
$
\BC^{{\Br}^{\scriptstyle *}}_{\trc}
$
filled from $R^*_{\phi,\varphi}(\Br_1,\Br_2)$ in (A5).
&(A11)
\cr
\hfill Compute SVD or eigendecompose:&
$
\BC^{{\Br}^{\scriptstyle *}}_{\trc}=
\BE^{{\Br}^{\scriptstyle *}}_{\trc}\,\,
\BPi^{\Br}_{\trc}\,\,
{\BE^{{\Br}^{\scriptstyle *}}_{\trc}}^T .
$
&(A12)\cr
\cr
~Tracer 3D Covariance Eigendecomposition ~~\cr
\hfill
Based on Kronecker product of (A10) and (A12):&
$
\BB_{\trc} = \BC^z_{\trc} \otimes
\BC^{{\Br}^{\scriptstyle *}}_{\trc}~.
$
&(A13)\cr
\hfill Sort eigenproducts $\lambda_z \lambda_r$ and truncate to subspace:&
$p$ dominant $\lambda_z \lambda_r \Rightarrow$
diagonal $\BPi_{\rm trc}\,$.
&(A14)\cr
\hfill Construct corresponding 3D eigenvectors:&
$p$ dominant $\Be_z \otimes \Be_r \Rightarrow$
rank-$p$ orthogonal $\BE_{\rm trc}\,$.
&(A15)\cr
\hfill Result is rank-$p$ decomposition for considered scale:&
$
\BB_{\rm trc}^p = \BE_{\rm trc} \BPi_{\rm trc}
\BE_{\rm trc}^T\,
$.
&(A16)\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\cr
{\capsfb b. Complete Primitive Equation Based Formulation}\cr
\hfill Form ensemble of perturbed {\it a priori} tracer fields:&
$
\xbj_{\rm trc} = \xb_{\rm trc} +  \BE_{\rm trc}
\BPi_{\rm trc}^{1 \over 2} \Bu^j~,~~j=1, ... ,r\, ,
$
&(A17)\cr
&\qquad where here $\Bu^j={\sqrt q} \Be^j\,$.\cr
\hfill
~ Form ensemble of unbalanced {\it a priori} state vectors:&
$
\xtidbj \doteq \big({\widehat \Bu^\b, \widehat \Bv^\b,
\xbj_{\rm trc}, \Bp^\b} \big)\,$,
&(A18)\cr
\hfill Run $r$ dynamical momentum adjustments:&
Integrate in (A1a) with $\Bx_0= \xtidbj$
to yield adjusted $\xbj$.
&(A19)\cr
\hfill
Form differences of PE adjusted fields, normalize:&
$
\BS^\b \doteq
\big[{ \xbone - \xb\,;~\ldots\,;~ \xbq - \xb} \big]
$;~ $\BS^* = \BN^{-1} \BS^\b$.
&(A20)\cr
\hfill
Compute SVD of rank $p \le r$:&
${\rm SVD}_p (\BS^\b) \doteq \BN \; {\rm SVD}_p (\BS^*)=
\BE \BSigma \V^T$;~ $\BE = \BN \BE^*$.
&(A21)\cr
\hfill Restart (A17--21) and evaluate convergence criterion:&
$
\rho
= {
{ \sum_{i=1}^{\tilde p} \, \sigma_{i} (\,
\widetilde {\BPi}^{1\over 2} \widetilde \E^{*^{\scriptstyle T}}
\E^* {\BPi}^{1\over 2}\,) }
\over
{ \sum_{i=1}^{p} \sigma_i ({\BPi})}
}~\ge \alpha\,$,
&(A22)\cr
&\qquad where $\alpha$ is a chosen limit
($1-\epsilon \le \alpha \le 1$) \cr
&\qquad and $\sigma_{i}(\cdot)$ selects the
singular value number $i$.\cr
\hfill
Result is principal error covariance estimate:&
$
\BB^p= \BGamma \,\, \BE \BPi\, \E^T \BGamma^T\,$.
&(A23)\cr
\noalign{\vskip5pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
%\vskip -0.3truecm

{\doublecolumns
In Sect.~3.c of LAL00, algorithms
for estimating the dominant components of 
error covariance matrices based on the functionals (A4b)
are outlined for three types of assumptions involving: a single scale,
multiple but independent scales and general multiple scales.
In each of these assumptions,
for each scale (and scale interactions in the third assumption),
four situations arise (e.g.~Figs.~2--3 in LAL00)
depending on how the vertical and horizontal
decompositions are obtained, either directly from data 
(EOFs of scale-restricted data residuals) or
from an analytical model
fit to data (eigendecomposition of a specified covariance matrix).
For more details
%on these assumptions and situations,
we refer to LAL00. The application dictates the choice among
options; the algorithm chosen here is summarized below and in Table~A.
\bigskip

\vskip 0.3truecm
{\bf A.2 Present construction of the {\it a priori} error subspace}

The discrete state vector
$\Bx$ $\doteq (\widehat \Bu, \widehat \Bv,
\BT, \BS, \Bp)$ $\in \bR^n$ contains gridded values of the
$\widehat u$, $\widehat v$, $T$, $S$ and $\psi$ fields
(see footnote~1).
%3D (dynamical) tracers, the temperature $T$ and salinity $S$;
%3D zonal $\widehat u$ and meridional
%$\widehat v$ internal baroclinic velocities;
%and, 2D barotropic transport stream function $\psi$.
The corresponding {\it a priori} error covariance 
estimate, $\BB$ (superscript $^\b$ omitted), is constructed as follows.
%Sample error or residual matrices are
%denoted by the symbol $\BS$. Their columns
%are field samples, the size of $\S$ depending on the context.

{\capsfn Observed portions}.
For the observed $T$ and $S$,
the assumption of multiple but independent scales is made and (A4b)
becomes (LAL00),
$$
\eqalignno{
C_{\phi \varphi}(\Br_1,\Br_2,z_1,z_2) \simeq \sum_{w=0}^{2} &
R^*_{\phiw,\varphiw}(\Br_1,\Br_2) ~~~~~ &\cr
& Z_{\phiw,\varphiw}(z_1,z_2)\,,
& ({\rm A5})}
$$
for $\phi$ and $\varphi$ any four combinations
of $T$ and $S$. For fair comparisons with the
univariate 2D scheme, for each scale~$w$ in (A5),
$R_{T_{\scriptstyle w} T_{\scriptstyle w}}^*$,
$R_{S_{\scriptstyle w} S_{\scriptstyle w} }^*$
and $R_{T_{\scriptstyle w} S_{\scriptstyle w}}^*$
are also assumed equal. 
The tracer submatrix $\BB_{\rm trc}$ of the error
covariance is then the sum of Kronecker products
of vertical $\BC^z_{\rm trc_{\scriptstyle w}}$
and horizontal $\BC^{{\Br}^{\scriptstyle *}}_{\rm trc_{\scriptstyle w}}$
tracer covariances,
$$
\BB_{\rm trc}= \sum_w \BB_{\rm trc_{\scriptstyle w}}
= \sum_w \BC^z_{\rm trc_{\scriptstyle w}} \otimes
\BC^{{\Br}^{\scriptstyle *}}_{\rm trc_{\scriptstyle w}}~.
\eqno({\rm A6})
$$
The algorithm used here to construct a $\BB_{\rm trc_{\scriptstyle w}}^p$
%(only one scale so that the sum over $w$ can be omitted)
is outlined in Table~A.a (A7--16):
for each scale $w$, $\BC^z_{\rm trc_{\scriptstyle w}}$
is obtained from data (A7--10)
and $\BC^{{\Br}^{\scriptstyle *}}_{\rm trc_{\scriptstyle w}}$
from an analytical model fit to data (A11--12).
%Presently, the parameters of the model for 
%$\BC^{{\Br}^{\scriptstyle *}}_{\rm trc_{\scriptstyle w}}$ 
%were fit to the data residuals
%by a combination of: previous experience (e.g.~Hecht {\it et al.},~1988;
%Milliff and Robinson, 1992; Robinson and Malanotte-Rizzoli, Eds., 1993),
%dynamical considerations (e.g.~Killworth, 1983),
%trial and error, and cross-validation (e.g.~Wahba and Wendelberger, 1980).
Normalization matrices are denoted by $\BN$.
They are block-diagonal: for each field,
the corresponding element of $\BN$ is the sample and spatial
averaged variance.
Further details
%,including specifics of the notation,
are in Sect.~3 and App.~B of LAL00.

{\capsfn Complete multivariate formulation}.
The velocity portions of $\BB$ are computed
through an ensemble of adjustment momentum integrations.
The algorithm for a given scale $w$ is outlined in Table~A.b
(specifics are in Sect.~4 of LAL00).
An ensemble of perturbed tracers
$\xbj_\trc \doteq (\BT^{\b^{\scriptstyle j}}, \BS^{\b^{\scriptstyle j}})$
is created (A17). The corresponding PE adjusted velocities are
obtained by momentum integrations (A19), starting from
unbalanced, e.g.~geostrophic, initial conditions (A18).
Adjusted variability samples can then be formed (A20),
normalized and organized by SVD (A21).
New adjustments (A17--21) are carried out
to increase $p$ until the criterion (A22)
determines that the number of samples is large enough.
Once this occurs, $\BB^p$ is obtained (A23). In (A22),
the pairs ($\E^*, \BPi={1\over q} \BSigma^2$) of rank $p$ and
($\widetilde \E^*, \widetilde {\BPi}={1\over \tilde q}
\widetilde {\BSigma}^2$)
of rank $\tilde p \le p$ correspond to
``new'' and ``previous'' estimates of $\BB^p$
(other criteria are in LER97).
In (A23), $\BGamma$ is a scaling, block diagonal matrix.


\vskip 0.3truecm
{\bf A.3 Error subspace mapping scheme}

With $\xb$, $\BB^p$ (A23)
and (A1a--b), the extremum of (A1c) yields
the {\it a posteriori} estimates
$\xa$, $\BE^\a$ and $\BPi^\a$,
hence $\BB^{\a^{\scriptstyle p}} \doteq \BE^\a\BPi^\a\EaT$.
For a direct inversion,
(A1b) is linearized in the vicinity of $\xb$
and $\xa$ is hypothesized to be a linear
function of $\xb$ and $\By^\o- {\cal H}(\xb)$,
each of which is assumed unbiased.
The extremum is then (LER97),
$$
\eqalignno{
\xa= &\, \xb +\BE \BPi \HpT
( \BH^p \BPi  \HpT+\BR )^{-1} ~~~~~~ &\cr
&~~~~~~~ \big( \By^\o- {\cal H}(\xb) \big),&
\rm{(A24a)}}
$$
%\cr
$$
\eqalignno{
\BU \BPi^\a {\BU}^T = &
\widetilde {\BPi}^\a \doteq
\BPi- \BPi \HpT \cr 
&~~~~~~~ (\BH^p \BPi \HpT+\BR)^{-1} \BH^p \BPi,
\quad ~~~~~~~&\rm{(A24b)}\cr
\BE^\a  =& \, \BE\, \BU~,
&\rm{(A24c)}}
$$
where $\BH$ is the linearization of
${\cal H}$ in the vicinity of $\xb$, $\BH^p \doteq \BH \BE$,
the columns of $\BU$ are ordered orthonormal eigenvectors of
$\widetilde {\BPi}^\a$ and $\BPi^\a$ is diagonal.
Details are in LAL00.
%In (A24a--c), the rank of $\BB^{\a^{\scriptstyle p}}$
%given by (A24b--c) is assumed not larger than that of
%$\BB^p$.
%The data residuals corresponding
%to a state in the complement of the 
%error subspace are neglected (Sect.~A.4). 
The equations (A24a--c) can also be obtained by formally
replacing covariances in the Kalman update step (e.g.~Gelb, 1974)
by their dominant eigendecompositions
$\BB^p$ and $\BB^{\a^{\scriptstyle p}}$. 
With the successive correction approach, (A24a--c)
are repetitively employed for each scale $w$.

\vskip 0.3truecm
{\bf A.4 A few relations with the classic Kalman update}

The difference between the Kalman update and (A24a)
arises because of the different gains:
$\BK \doteq \BP^\b\BH^T\;[\BH \BP^\b\BH^T+\BR]^{-1}$
for the Kalman gain and 
$\K^p \doteq \BP^{\b^{\scriptstyle p}}
\BH^T\, [\BH\BP^{\b^{\scriptstyle p}}\BH^T+\R]^{-1}$ 
for the subspace gain. Three cases are considered next
as a function of the measurement properties $\BH$ and $\R$.
The concept of Data Optimal (LER97) is used.
The rank-$p$ eigendecomposition 
$\BP^{\b^{\scriptstyle p}}=\E \BPi \E^T$ and its complement,
$\BP^{\b^{\scriptstyle c}}
\doteq \BP^\b - \BP^{\b^{\scriptstyle p}} = \BE^c \BPi^c \EcT$,
are assumed known. 

{\capsfn Measurements only sample the error subspace.}
The fixed-time Data Optimal subspace
of rank-$p$ is then observed. Keeping the distinction between normalized
and non-normalized quantities, this
situation can be defined by $\BH=\BGamma\,\E^{*^{\scriptstyle T}}
\BN^{-1}$, where 
$\BGamma \in\bR^{m \times m}$ is a matrix of coefficients
and $\BN$ the normalization matrix (A21).
Since $\BE^{*^{\scriptstyle T}} \BE^{c^{\scriptstyle *}}=0$,
$$\eqalignno{
\K=&\Pb \BH^T \,(\BH\Pb\BH^T+\R)^{-1}\cr
=& \E \BPi \BGamma^T\,
(\BGamma\BPi\BGamma^T+\R)^{-1}
&({\rm A25a})\cr
\K^p=&\Pbp\BH^T\, (\BH\Pbp\BH^T+\R)^{-1}\cr
=& \E \BPi \BGamma^T\,
(\BGamma\BPi\BGamma^T+\R)^{-1}
=\K\,.
~~~~~~ &({\rm A25b})}
$$
Hence, $\K-\K^p=0$, and from (A24b--c),
the error variance reduction achieved,
$$
\eqalignno{
\BP^{\a^{\scriptstyle p}}=&\E 
\big({\BPi - \BPi \BGamma^T\, ( \BGamma\BPi\BGamma^T+\R)^{-1}
\BGamma \BPi}\big ) \E^T ~~~~~~~ & \cr
&+ \BE^c \BPi^c \EcT = \P^\a~,
&({\rm A25c})}
$$
is optimal in the full Kalman update sense. 
For such ideal measurements, the {\it a posteriori} states
and errors of the minimum error variance and present 
estimates (A24a--c) are equal.

{\capsfn Measurements only sample the complement of the error subspace.}
The complement of the fixed-time, rank-$p$ Data Optimal
subspace is then observed. This situation
can be defined by $\BH=\BGamma\, \EcsT \BN^{-1}$, where
$\BGamma \in\bR^{n-m \times n-m}$ is a matrix of coefficients,
$$
\eqalignno{
\K=&\Pb\BH^T\, (\BH\Pb\BH^T+\R)^{-1}\cr
=& \E^c \BPi^c\BGamma^T\, 
(\BGamma\BPi^c\BGamma^T+\R)^{-1} ~~~~
&({\rm A26a})\cr
\K^p=&\Pbp\BH^T\, (\BH\Pbp\BH^T+\R)^{-1}=0\,.~~~~~~
&({\rm A26b})}
$$
Hence, $\K-\K^p=\K$ and from (A24b--c), 
$\BP^{\a^{\scriptstyle p}}=\BP^{\b^{\scriptstyle p}}$.
If the domain of $\BH$
spans the complement of the error subspace,
the data has no influence on (A24a--c).
Nonetheless, by construction, these data are associated
with negligible components $\BP^{\b^{\scriptstyle c}}$ of the {\it a priori} error. They would
thus improve $\xb$ by negligible amounts and can be neglected.
 
{\capsfn Measurements sample the full error space.}
This situation is exemplified by
a complete measurement array $\BH=\I$, which yields
$$\eqalignno{
\K=&\Pb\,(\Pb+\R)^{-1}
=\big(\E\BPi \E^T+\E^c \BPi^c \EcT\big) ~~~~~~&\cr
&~~~~
(\E\BPi \E^T+\E^c \BPi^c \EcT+\R)^{-1} 
&({\rm A27a})\cr
\K^p=&\Pbp\,(\Pbp+\R)^{-1}\cr
=&\E \BPi \E^T\,
(\E\BPi \E^T+\R)^{-1}~.
&({\rm A27b})}
$$
The Frobenius norm of $\K-\K^p$ becomes negligible in comparison to
that of $\K$ when the dominant eigenvalue of $\BPi^c$
becomes negligible in comparison to the norm of $\BR$. 
For example, for $\BR=r \BI$, the weight of the
{\it a priori} error vector $i$ in the Kalman gain (A27a)
is proportional to $\lambda_i \over (\lambda_i + r)$.
For given accuracy, the eigenvectors corresponding to $\lambda_i$ 
negligible in comparison to $r$ are unnecessary.
On the other hand, if the largest eigenvalue of $\BPi^c$
is not negligible in comparison to $r$,
$\K^p$ may not account for most of $\K$. Today,
this situation is very uncommon for lack of good data.
Should it occur, to avoid spurious effects of the truncation,
one must increase $p$ until
all non-negligible data components are considered.
 
\vskip 1.truecm
\centerline{APPENDIX B}
\centerline{\bf Relevant scalings and}
\centerline{\bf non-dimensional numbers}
\vskip 0.3truecm

Table B defines the main non-dimensional numbers that are
relevant for the present mesoscale to subbasin-scale
study. The notation is standard (Gill, 1982; Pedlosky, 1987; 
Cushman-Roisin, 1994): $g$ is the acceleration due
to gravity; $L$ and $D$ the horizontal and vertical length-scales
of motion; $U$ and $W$ the horizontal and vertical
characteristic velocity scales;
$f_0$ and $\beta_0$ the Coriolis parameter 
and its northward gradient at latitude~$\theta_0$;
$\rho_s(z)$ the background density
profile; and, $\rho_0$ a constant reference density.
In addition, $N_s = \sqrt { - { g \over \rho_0}
{\delta \rho_s \over \delta z}}$ is the Brunt-V\"ais\"al\"a frequency,
$R_D = N_s D / f_0$ the internal Rossby
radius of deformation, $R_{\rm ext}= \sqrt {g D} / f_0$ the external Rossby
radius of deformation
and $\nabla_h H$ the horizontal gradient of the 
bottom topography $H$.
Except for $g$, all of these characteristic scales
are assumed regional: they vary with the location and process considered.

For the synoptic and horizontally
isotropic quasi-geostrophic scales in a stratified ocean,
$Ro$, $Bu_{\rm ext}^{-1}$ and $s \ll 1$, but the other numbers
may take a range of values: in particular,
$Bu$ and $\alpha$ increase with the strength of the stratification
and topographic slope, respectively (e.g.~Pedlosky, 1987).
At these scales, using the internal energy (density)
equation to estimate $W$, one can also show that
$Ro\over Bu$ is the ratio of vertical convergence to horizontal
divergence $W L \over U D$, and that
the characteristic horizontal variations
in density $\Delta \rho$ and vertical variations in background density
$\Delta \rho_s$ are related via,
${\Delta \rho \over\rho_0} = {Ro\over Bu} \,
{\Delta \rho_s \over\rho_0}$.
\vfill\eject}
%\bigskip

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf

\centerline{Table B. Main non-dimensional numbers/parameters}
\vskip -0.2truecm
$$\vbox{\halign{
#\hfil&~~#\hfil&\quad#\hfil&\quad#\hfil\cr
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip4pt}
Number/Parameter & Definition & Ratio\cr
\noalign{\vskip4pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
\noalign{\vskip4pt}
\cr
Rossby: & $Ro= {U \over f_0 L}$ & Relative vorticity to planetary vorticity\cr
\noalign{\vskip -1pt}
& & (also nonlinear advection to Coriolis force)\cr
\noalign{\vskip3pt}
External Burger: & $ Bu_{\rm ext}= {g D \over f_0^2 L^2} =
{R_{\rm ext}^2 \over L^2}$ &
Relative vorticity to external vertical\cr
\noalign{\vskip -4pt}
& & vortex-tube stretching (of free-surface)\cr
\noalign{\vskip3pt}
Internal Burger: & $ Bu={N_{\scriptstyle s}^2 D^2 \over  f_0^2 L^2} = 
{R_D^2 \over L^2} $ & Relative vorticity to internal
vertical \cr 
\noalign{\vskip -4pt}
& & vortex-tube stretching (of isopycnal surfaces)\cr
\noalign{\vskip3pt}
Planetary-$\beta$: & $\beta=\beta_0 \cdot {L^2 \over U}$ & 
Planetary (ambient) vorticity gradient \cr
& & to relative vorticity gradient\cr
\noalign{\vskip3pt}
Topographic-$\beta$: & $\alpha = \vert \nabla_h H \vert 
{f_0 \over D} \cdot {L^2 \over U}$ & 
Topographic (ambient) vorticity gradient \cr
& & to relative vorticity gradient   \cr
\noalign{\vskip3pt}
Slope or topographic: & $ s =\vert \nabla_h H \vert \cdot {L \over D}
= \alpha\, Ro $ &  Topographic gradient to motion aspect ratio\cr
\noalign{\vskip4pt}
\noalign{\hrule}
\noalign{\vskip1pt}
\noalign{\hrule}
}}$$}
%\vskip-0.5truecm


{\doublecolumns
\vskip 1.truecm
\centerline{APPENDIX C}
\centerline{\bf Dynamical model and parameters}
\vskip 0.3truecm

The adjustment momentum integrations and forecasts
are carried out using
the numerical primitive equation model
(see footnote~1) of HOPS (e.g.~Robinson, 1996, 1999).
%The state variables are the (dynamical) tracers,
%the temperature $T$ and salinity $S$,
%the barotropic transport stream function $\psi$, and
%the zonal ($x$) and meridional ($y$) zero vertical-average
%internal velocities, $\widehat u$ and $\widehat v$, respectively.
%The corresponding five fields
%are henceforth called the PE fields.
The present values of the main parameters are listed in Table~C.
\bigskip

{\eightpoint
\baselineskip=9pt plus 1pt minus 1pt
\capsf

\centerline{Table C. Main parameters of the primitive equation model}
\vskip 0.2 truecm
%\centerline{\psfig{figure=\tab_src/table1chap8.eps,width=6truein}}
%\centerline{\psfig{figure=\tab_src/Table1.ch8new.eps,width=6truein}}
%\centerline{\psfig{figure=\tab_src/Table1.ch8new2.eps,width=4.5truein}}
\centerline{\psfig{figure=\tab_src/Table1.ch8new2.eps,width=3.5truein}}
}
\noindent
Horizontally, the numerical filtering and parameterization
of subgrid-scale processes is based on a Shapiro filter (e.g.~Shapiro, 1970):
the symbols, $F_u$, $F_v$, $F_T$, $F_S$ and $F_{\,\overline{\omega}_t}$
in Table~C
give the parameters of the filter for each of the five prognostic equations
($\overline{\omega}_t$ is the rate of change of barotropic vorticity) .
The corresponding effective Laplacian horizontal eddy viscosity
$A_h$ and diffusivity $K_h$ (LER97) decay with space scales.
Under 20$\,$km, they are both $3\; 10^8$~\cm2secm1.
At 50$\,$km scales, $A_h=1.6\; 10^5$~\cm2secm1 and
$K_h=3.8\; 10^6$~\cm2secm1.
Beyond 100~km scales, they are both almost null.
The open-boundary conditions are of the Orlanski radiation type
(Orlanski, 1976),
either implicit or partially explicit (ORI/ORE$_{1\over2}$).
The vertical mixing is a Laplacian mixing
(eddy viscosity $A_v$ and diffusivity $K_v$).
Convective adjustment is utilized when the water column
is statically unstable (eddy viscosity
$A_v^{\rm cvct}$ and diffusivity $K_v^{\rm cvct}$).
Across coastlines, the normal flow and
tracer flux are set to zero. Along coastlines,
the tangential flow is dissipated at the coast, using a
Rayleigh friction of relaxation time
$\tau_c$ and Gaussian decay $L_c$.
At the bottom, for the momentum
equations, a dynamic stress balance (drag coefficient $C_d$)
and a Rayleigh friction (relaxation time
$\tau_b$ and Gaussian decay $H_b$) are applied.
%The specifics and implementations of these processes
%and conditions are discussed in LER97.
In the cases illustrated in this work,
stochastic forcings are not utilized.
As discussed in LER97, for the present Levantine application,
such forcings increase the sub-mesoscale variability
but do not modify the main
mesoscale and subbasin-scale results.



\end
