# Assessment of observational networks with the Representer Matrix Spectra method-application to a 3D coastal model of the Bay of Biscay

Matthieu Le Hénaff, Pierre de Mey, Patrick Marsaleix

► To cite this version:

Matthieu Le Hénaff, Pierre de Mey, Patrick Marsaleix. Assessment of observational networks with the Representer Matrix Spectra method-application to a 3D coastal model of the Bay of Biscay. Ocean Dynamics, 2009, 59 (1), pp.3-20. 10.1007/s10236-008-0144-7. hal-02110195

HAL Id: hal-02110195

https://hal.science/hal-02110195

Submitted on 5 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessment of observational networks with the Representer Matrix Spectra method—application to a 3D coastal model of the Bay of Biscay

Matthieu Le Hénaff · Pierre De Mey · Patrick Marsaleix

**Abstract** The development of coastal ocean modeling in the recent years has allowed an improved representation of the associated complex physics. Such models have become more realistic, to the point that they can now be used to design observation networks in coastal areas, with the idea that a “good” network is a network that controls model state error. To test this ability without performing data assimilation, we set up a technique called Representer Matrix Spectra (RMS) technique that combines the model state and observation error covariance matrices into a single scaled representer matrix. Examination of the spectrum and the eigenvectors of that matrix informs us on which model state error modes a network can detect and constrain amidst the observation error background. We applied our technique to a 3D coastal model in the Bay of Biscay, with a focus on mesoscale activity, and tested the performance of various altimetry networks and an in situ array deployment strategy. It appears that a single nadir altimeter is not efficient enough at capturing coastal mesoscale physics, while a wide swath altimeter would do a much better job. Testing various local in situ array configurations confirms that adding a current meter to a vertical temperature measurement array improves the detection of secondary variability

modes, while shifting the array higher on the shelf break would obviously enhance the model constraint along the coast. The RMS technique is easily set up and used as a “black box,” but the utility of its results is maximized by previous knowledge of model state error physics. The technique provides both quantitative (eigenvalues) and qualitative (eigenvectors) tools to study and compare various network options. The qualitative approach is essential to discard possibly inconsistent modes.

**Keywords** Stochastic modeling · Model errors · Data assimilation · Array design · Coastal ocean

## 1 Introduction

The development of complex ocean models, thanks to the increase in computational facilities along with the increasing number of oceanic observations, has helped to get a better representation of the ocean physics and variability. Data assimilation techniques, first developed in meteorology, have been adapted to oceanography with the main objective of forecasting the ocean state, at least at basin scale. The goal of observation networks design is nowadays not so much about discovering unknown features of the ocean physics but rather about correctly detecting some specific aspects of it, such as large-scale currents, thermocline depths, eddies, etc. These relevant observations would in return help models better represent the associated physics, leading to a more precise analysis of the ocean state. One satisfactory way to identify a “good” network would therefore be to test its efficiency in (1) detecting and (2) controlling the errors of the state of a model suitable for whatever physics the network has been deployed to measure.

M. Le Hénaff (✉) · P. De Mey
LEGOS, Université de Toulouse/CNRS,
14 avenue Edouard Belin,
31400 Toulouse, France
e-mail: Matthieu.Le.Henaff@legos.obs-mip.fr

P. Marsaleix
Laboratoire d'Aérologie,
14 avenue Edouard Belin,
31400 Toulouse, France
The idea of model state error control is the base of many studies about array design in recent years, most of them accompanying the development of data assimilation in oceanography. Bennett (1985) and McIntosh (1987) first used a variational data assimilation technique to study tide gauge deployment. Within that approach, Bennett (1985, 1990) defined the notion of *array modes*, which gives an insight of the impact an observation array has on a model. Barth and Wunsch (1990) and Barth (1991) carried out optimized array design by the minimization of a cost function, through simulated annealing or genetic algorithms. The latter approach has been followed by Hernandez et al. (1995) to optimize a drifter cast strategy. Many observation network studies rely on observing system experiments, assimilating data from true (OSEs) or simulated observations (OSSEs) to quantify the performances of an array and compare various networks. Hackert et al. (1998) used OSEs involving a reduced order Kalman filter to define the best locations for the Pilot Research Moored Array in the Tropical Atlantic (PIRATA) network. Guinehut et al. (2002) used an optimal interpolation (OI) technique to address the float deployment problem in the North Atlantic. Schiller et al. (2004) used OSSEs to design an Argo float deployment in the Indian Ocean. In the coastal ocean, different altimeter configurations have been compared using OSSEs to control the model error on the continental shelf (Mourre et al. 2006). In the Mediterranean Sea, Raicich (2006) used OI to compare the various available in situ arrays. The majority of these studies rely on data assimilation and comparison of the different array performances through data assimilation, in error reduction, or correlation to a “true” state. Performing data assimilation is computationally expensive for the complex physics and fine scales characterizing coastal processes. Using simplified data assimilation schemes or interpolation techniques to perform array design is relevant to study open ocean processes of relatively simple variability but is probably not ideal to study highly nonlinear, non-isotropic physics, for which the model error budget is much more complex. Advanced data assimilation techniques such as 4D-Var or ensemble Kalman filter (EnKF) are probably more adapted when assimilating data in a coastal model, but the computation cost becomes very high if many network options are to be tested. We would certainly be interested by an approach assessing the relative performances of various observation networks on a complex coastal zone at a relatively low computational cost, while keeping the idea that a “good” network must be able to detect and provide an efficient constraint of model state error.

Observations are usually sparser in coastal ocean than in the open ocean. For instance, routine altimetry observation data are filtered to keep large space and time scales variability and are not satisfactory in coastal zones (Anzenhofer et al. 1999). Improving altimetric data close

to the coast is a challenge currently addressed by several studies, for example, Vignudelli et al. (2005) and Bouffard et al. (2008). Recent massive Argo floats deployments have mainly been designed to target basin open ocean circulation and mesoscale, although some of them eventually end in coastal zones. Kurapov et al. (2005) on the Oregon Shelf or Auclair et al. (2003) in the Gulf of Lions studied the complexity of coastal dynamics in a data assimilation framework. The need for observation data in coastal zones is essential, and it is important to make sure that future observation arrays are optimally designed.

Ensemble simulations provide an efficient way to estimate the complex model state (or *forecast*) error covariance matrix, as in the EnKF (Evensen 1994). In this paper, we use the ensemble approach to calculate several important quantities characterizing how a given observation network is able to detect state errors. By means of studying the *scaled representer matrix* associated to a network, it is possible to examine the dominant error modes a network can measure. The associated eigenmodes are used to define *modal representers*, which give an insight of the potential a network has in terms of model state error control. These are close to the “array modes” defined in Bennett (1985, 1990). The Representer Matrix Spectra of various networks can quantitatively be compared to one another, while a more qualitative analysis can be carried out, thanks to the associated eigenvectors and modal representers. This technique, which we coin “Representer Matrix Spectra” technique in this paper, is explained and illustrated in Section 2. The application on a coastal ocean model of the Bay of Biscay is depicted in Section 3, with two examples illustrating how useful this technique can be, provided that the dominant modes are robustly associated with known dynamical processes in the model state error subspace. Conclusions are addressed in Section 4.

## 2 The RMS technique

### 2.1 Theory

Let us denote $\mathbf{x}$ as the state vector under study. It is composed of a set of $n$ variables of interest, distributed in space and time to adequately define the state of the system and its evolution. More precisely, $\mathbf{x}$ is an *augmented state* vector, to distinguish it from the “classical” definition of a state vector, which is usually instantaneous.

We consider now a set of $p$ observations of the augmented state vector, denoted $\mathbf{y}^o$. If $\mathbf{x}^t$ is the true state of the ocean, $\mathbf{y}^o$ can be defined from this true state as

$$ \mathbf{y}^o = H(\mathbf{x}^t) + \varepsilon \quad (1) $$

where $H$ is the observation operator, which projects the state vector onto the observational space, and $\varepsilon$ is the
observation error, assumed to have zero mean (no bias) and a known covariance equal to $\mathbf{R}$.

Our purposes in this study are to examine whether and under which circumstances one particular array ($H, \mathbf{R}$) can be said to be “objectively satisfactory,” and whether, given two arrays ($H_1, \mathbf{R}_1$) and ($H_2, \mathbf{R}_2$), one of them can be said to be more efficient than the other.

We assume we have an a priori estimate—a guess—of the state and its associated error statistics. The guess can be any prior estimate such as a previous model simulation or forecast. The case where we would have no a priori knowledge whatsoever is not so interesting because, in that case, any observational array would bring valuable information proportionally to its cost. Our a priori guess writes:

$$ \mathbf{x}^g = \mathbf{x}^t + \eta \quad (2) $$

where $\eta$ is the prior model state (guess) error, assumed to be of zero mean, with covariance equal to $\mathbf{P}^g$, the guess error covariance matrix.

We assume here that $\eta$ and $\varepsilon$ are uncorrelated with each other, which is not necessarily true. The information increment brought by the observation data on top of our guess is called the innovation vector $\mathbf{d}$:

$$ \mathbf{d} = \mathbf{y}^o - \mathbf{y}^g = \mathbf{y}^o - H(\mathbf{x}^g) \approx \varepsilon - \mathbf{H}\eta \quad (3) $$

where $\mathbf{H}$ is the linearized observation operator at point $\mathbf{x}^t$. The second-order statistics of $\mathbf{d}$ can be used to characterize the amount of discrepancy brought in by the observational array. From Eq. 3, we have

$$ \langle \mathbf{d}\mathbf{d}^T \rangle = \mathbf{H}\mathbf{P}^g\mathbf{H}^T + \mathbf{R} \quad (4) $$

The right-hand side part of Eq. 4 is a sum of two matrices. The second one is the observational error covariance matrix; the first one, known as the *representer matrix*, is the projection, onto the observational space, of the guess error covariance matrix. By construction, these two matrices are symmetric definite positive. The sum of these two matrices is referred to as the *stabilized representer matrix*. If $\mathbf{R}$ dominates, the discrepancies between the observations and the model are mostly due to observational errors, and the array is pretty much useless. If $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ dominates, then the discrepancies are mostly due to guess, and the observational array would probably be useful to assess and then possibly correct the model in a data assimilation step. The Representer Matrix Spectra technique is based on the comparison of both the model and observations covariance matrices, without performing data assimilation. It is inspired by the *array modes* technique described by Bennett (1985, 1990). In a variational data

assimilation framework, Bennett diagonalizes the stabilized representer matrix:

$$ \mathbf{H}\mathbf{P}^g\mathbf{H}^T + \mathbf{R} = \mathbf{V}\mathbf{S}\mathbf{V}^T \quad (5) $$

In his case, the largest eigenvalues of this matrix define the dominant modes for the correction of the model through the data assimilation step, and this correction spreads to the model space through the matrix $\mathbf{\Psi} = \mathbf{P}^g\mathbf{H}^T\mathbf{V}$. $\mathbf{\Psi}$ is a matrix whose columns the author calls *array modes*.

Here, we simply wish to compare $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ and $\mathbf{R}$ without performing data assimilation. The study of the spectrum of $\mathbf{H}\mathbf{P}^g\mathbf{H}^T + \mathbf{R}$ does not tell us in itself which prior error modes are detected by the array and which ones are not seen amidst the observational noise. A solution would be to compare the spectra of both $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ and $\mathbf{R}$ to determine the modes for which the eigenvalues of $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ are larger than those of $\mathbf{R}$, which is a rather simple problem in the case of independent, homogeneous observational errors. Usually, in data assimilation, observation errors are considered independent, with standard deviation consistent with the assumed measurement error $\varepsilon^o$. In that simple case, $\mathbf{R}$ is diagonal with constant values equal to $\langle(\varepsilon^o)^2\rangle$. If an eigenvalue of $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ is larger than $\langle(\varepsilon^o)^2\rangle$, then the corresponding prior error mode would be detected amidst the measurement noise, and its information would be valuable for further assessment and assimilation.

In this paper, we tackle the more complex case of a non-diagonal, non-constant $\mathbf{R}$ matrix. To that purpose, we define a new representer matrix $\chi$, which is scaled by the observation error covariance and which we call the *scaled representer matrix*:

$$ \chi = \mathbf{R}^{-1/2}\mathbf{H}\mathbf{P}^g\mathbf{H}^T\mathbf{R}^{-1/2} \quad (6) $$

The problem of comparing the spectra of $\mathbf{H}\mathbf{P}^g\mathbf{H}^T$ and $\mathbf{R}$ is then replaced here by the problem of comparing the spectra of $\chi$ and of the identity matrix, as we multiply both matrices by $\mathbf{R}^{-1/2}$ on the left and right sides.

Diagonalizing $\chi$, we have access to the *scaled error modes* associated with the array under study:

$$ \chi = \mu\sigma\mu^T \quad (7) $$

$\sigma$ is a diagonal matrix containing the ordered eigenvalues $\lambda_i$ of the $\chi$ matrix. This time, we are interested in eigenvalues $\lambda_i$ that are larger than one; each eigenvalue is associated with a distinct eigenmode in observational space, contained in the columns of $\mu$. These modes exhibit the prior error patterns associated with each eigenvalue in observation space.

From Eq. 7, and the orthonormality of the scaled modes, we get

$$ \mu^T\chi\mu = \mu^T\mathbf{R}^{-1/2}\mathbf{H}\mathbf{P}^g\mathbf{H}^T\mathbf{R}^{-1/2}\mu = \sigma \quad (8) $$

$\sigma$ appears here as a rotated scaled representer matrix, in the new basis defined by $\mu^T$. In the Kalman filtering formalism,
$P^gH^T$, whose columns are vectors of covariance between the model at an observation location and the rest of the model state space is named *matrix of representers*. Here $\rho = P^gH^T R^{-1/2}\mu$ can be seen as a matrix of representers for the modes, i.e., a matrix of modal representers. Indeed, each column of $\rho$ contains the patterns associated with a scaled mode, but in model state space. This is a very valuable piece of information giving insight into the impact an observation network would have on the model state. As this analysis is performed altogether in space and time, the scaled modes may give information about the relevance of the spatial resolution of a network and of its time sampling. The modal representers have their signature in the full model space and time.

Modal representers obviously remind us of Bennett's array modes. They are different in that modal representers are homogeneous to the state vector variables, while array modes are homogeneous to the product between the considered state and observed variables and in that the eigenmodes from which they are built are different. The two approaches are similar only if $R$ is diagonal and homogeneous. In this case, we write $R = rI$, and therefore, $HP^gH^T$, $\chi$, and $HP^gH^T + R$ all have the same eigenvectors. We simply have $\rho = r^{-1/2}\Psi$.

If one wants to compare two networks defined by ($H_1$, $R_1$) and ($H_2$, $R_2$), one has to calculate the scaled representer matrix and the associated spectra. The network that has the largest number of eigenvalues above unity can be defined as the *objectively* best network, in the sense that it is able to detect more of the scaled modes. Because of the scaling by the $R$ matrix, the direct comparison of eigenvalue spectra is relevant even if the number of observations is not the same for both networks. Moreover, the scaling by $R$ theoretically allows comparing networks that do not measure the same physical variables, provided that such a comparison makes physical sense. As a complement to spectral analysis, visual inspection of the modes associated with the largest eigenvalues of the spectrum is probably necessary, since it provides some form of verification that the appropriate error-space physics are captured by the observational array.

Our definition of the *objectively* best network is not universal. One could, for instance, consider that the most efficient network is the one that explains the largest scaled error variance, i.e., which has the largest sum of eigenvalues of $\chi$ or that it is the one explaining the largest part of the scaled variance associated to some physical or dynamical features one is interested in. The former criterion is an acceptable alternate quantitative criterion, but would not facilitate physical interpretation in terms of the modal representers. The second criterion can use our approach as long as the state space has been redefined as a projection onto the desired physics.

Practical aspects of the technique, e.g., how calculations are carried out and how $P^g$ is accessed, are reviewed later in the paper.

## 2.2 Example

Let us consider a single-variable 2D model at a time $t$, at which time the prior error covariance matrix $P^g$ is assumed to be known, and two observation networks, one triplet in the North–South direction ($H_1$), the other triplet in the East–West direction ($H_2$), both observing at that time $t$ only. We consider the observation errors to be independent and homogeneous so that the observational error covariance matrix $R$ is diagonal and the same for both networks. Figure 1a shows a map of the standard deviation of the prior error at observation time and the location of the two triplets. As the standard deviation patch appears to be isotropic, one may not look further and infer that both networks would be equivalent in constraining the model state error subspace. However, from the spatial correlation of prior errors (Fig. 1b), we can see that error correlation scales are larger in the E–W direction. This is akin to anisotropy found in shelf seas, due to the neighboring presence of the coast or shelf break. Anisotropy in the model state error subspace directly affects the performances of the two networks and translates into distinct scaled Representer Matrix Spectra (Fig. 1c). The N–S triplet spectrum is characterized by three eigenvalues higher than 1, while the E–W triplet spectrum has only one. In effect, the N–S triplet is able to detect more degrees of freedom

<figure><img src="image_2.png" /><figcaption>**Fig. 1** a Model sea level error standard deviation. Two observation networks are plotted: east–west (E–W) network with circles, north–south (N–S) network with squares with a central cross. The central point belongs to both. b Correlation in sea level error between the central point and the rest of the domain. c Scaled Representer Matrix Spectra of E–W (thin line with squares) and N–S networks (thick line with circles). The dot-dashed line is the spectrum of the identity matrix, equal to 1</figcaption></figure>
(df s), while the E–W triplet (redundantly) constrains only 1 df of state error. This interpretation is confirmed by Fig. 2, showing, for each network, the three modal representers contained in ρ, as defined in the previous section, associated with eigenvalues of the scaled representer matrix. While the first modal representer exhibits the same pattern for both networks, consisting of a large central coherent zone in the E–W direction, the second and third modes differ. In the case of the N–S network, modal representers 2 and 3 do bring fresh information with respect to mode 1. In the case of the E–W network 2, modes 2 and 3, of low amplitude, are only modifiers of mode 1 with little new information, given the longer correlation scales in that direction.

This simple example illustrates how two networks with
equal cost and seemingly equivalent return, in an attempt to
measure an apparently isotropic patch of error variance, are
discriminated by the RMS technique in terms of detection,
and potential control of, model state error structures.

3 Application to a 3D coastal model of the Bay of Biscay

We now aim at illustrating the RMS technique in a realistic
case of coastal ocean array design. We are particularly
interested in observing short-scale, low-frequency mesoscale
structures in the vicinity of the coastal transition zone
between the deep ocean and the shelf. We first set up a
coastal ocean model of the Bay of Biscay. The guess error
covariance matrix is determined by stochastic modeling with
perturbations of wind forcing. We then use the RMS
technique to test several remote sensing (altimetry) scenarios
and in situ arrays on a subdomain.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Modal representers in sea level for the E–W (top) and N–S (bottom) networks. *From left to right* First, second, and third modes. *Continuous lines* are for positive values, *dashed lines* for negative values. Note that the colorbar scaling of modes 2 and 3 is different for both networks</figcaption></figure>

3.1 Model implementation

The 3D, primitive equation free surface Symphonie model
is set up in the Bay of Biscay, from the southwest of Cape
Finisterre to the middle of the English Channel (cf. Fig. 3),
during July and August 2004. We use the Bay of Biscay
configuration of the model (Marsaleix et al. 2008), with
2-day averaged outputs. The model uses bulk formulae to
calculate heat fluxes and a sigma-step vertical scheme. The
grid has a 3-km horizontal resolution and 45 vertical levels,
with higher density toward the surface. Meteorological
forcings are the Meteo-France Aladin model outputs, with a
10-km resolution over the domain and a 3-h temporal
resolution. The model initial state and open-boundary fields
are taken from the Mercator Psy2v1 model weekly data and
have been extrapolated on our grid using the VIFOP tool
(Auclair et al. 2006). River outflows from the Loire, the
Gironde, and the Adour rivers are also modeled from in situ
data, even if the study time is a bit short to fully study their
influence.

These modeling parameters have been designed to
correctly represent the coastal area mesoscale. The 2-day
output period has been chosen to smooth both the inertial
effects (period of about 16 h at this latitude) and the diurnal
and semi-diurnal tide effects (tides were eventually not
modeled but had at first been intended to be).

Figure 3 shows the bathymetry of the domain and the
main dominant dynamical features of the zone. A slope
current in the upper layer follows the continental shelf,
from North of Spain to South of Ireland (Pingree and Le
Cann 1989, 1990). Along the Cantabrian and Armorican
Slopes, variations along the steep bathymetry line are also
places of mesoscale activity with formations of meanders
<figure><img src="image_4.png" /><figcaption>**Fig. 3** a Bathymetry of the Bay of Biscay (m) with contours indicated on the color bar. Arrows are dominant dynamics features: slope currents (continuous), large-scale anticyclonic circulation (dot), and eddies (circles). b Zoom on the zone framed in a</figcaption></figure>

and eddies, which eventually get into the open ocean (Pingree and Le Cann 1992a, b; Van Aken 2002). The Landes Plateau region with a smoother slope at the shelf edge, framed by two steep canyons, the Cap Ferret Canyon in the North and the Cap-Breton Canyon in the South, is typical in that it generates high mesoscale activity (Bardey et al. 1999; Garcia-Soto et al. 2002). The open ocean zone is influenced by the large-scale anticyclonic circulation, associated with the North Atlantic subtropical gyre (Pingree 1993; Van Aken 2001). Typical water masses of the zone are masses from the North Atlantic advected by the large scale currents from the subtropical gyre, and salty intermediate water formed in the Mediterranean Sea advected along Portugal and spreading in the Bay of Biscay (Arhan et al. 1994; Van Aken 2000).

The summer period has been chosen because it is supposedly a period with lower, but significant, meteorological variability. Actually, summer 2004 was characterized by bad weather, with several low-pressure systems crossing the zone. On July 7th, a severe storm forms on the Bay of Biscay and strikes Brittany, as seen on Fig. 4a. Then, a calmer period follows until mid-August, when two

major low-pressure systems cross the zone from West to East, as shown in Fig. 4b,c. These two events are very similar in shape and affected the entire modeling domain. Figure 4d shows the time evolution of the mean atmospheric pressure over the domain. We clearly see the three aforementioned events, with relative high pressure between mid-July and mid-August.

## 3.2 Validation

The model is to be validated in terms of hydrology and dynamics in order to get sufficient trust into the way the mesoscale features are reproduced. Vertical distribution of density drives the mesoscale activity, as its features, meanders, or eddies have a horizontal extension typical of the first Rossby radius of deformation. This radius depends on the Brunt–Väisälä frequency $N(z)$:

$$N^2(z) = - \frac{g}{\rho_0} \frac{d\rho}{dz} \qquad (9)$$

and thus on the density profile of the water column, driven by salinity and temperature. On Fig. 5a are shown six

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Ten-meter wind and atmospheric pressure over the zone on a July 7, b August 12, and c August 19. Wind reference arrow is 10 m s⁻¹, P is in hPa. d Averaged pressure over the domain during the study period. The dashed line is the time mean of the series, 1,016 hPa. Dates of a, b, and c are in dot-dashed lines</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5** a **Squares** Available observations from in situ drifters during the study period (see Table 1). **Thin lines** 200, 2,000, and 4,000 m isobaths. **Circles** Point 1 (on the shelf) and point 2 (in the interior) used in Fig. 9. **b** θ–S diagram from the observation (**thin line with crosses**), in **white square** on **a**, and the corresponding model output (**thick line with circles**) on July 12. The **dotted lines** are constant linearized σ₀ values</figcaption></figure>

every 0.5 g m⁻³. Depths of local extrema indicated for each curve, as well as the depth of the junction of the two. **c** Surface to 200-m depth Brunt–Väisälä frequency N(z) (10⁻³ s⁻¹) for the observation (**thin line with crosses**) and the corresponding model output (**thick line with circles**) on the same day. **d** Same as **c** but from the surface to 2,000-m depth

observations by Argo floats, available on the Coriolis website (http://www.coriolis.eu.org/cdc/) for summer 2004, with enough data to get temperature and salinity profiles. Figure 5b shows the θ–S diagrams from both in situ data and the corresponding Symphonie output, for the first of the observation points (white square on Fig. 5a). We see major differences on this plot, mainly in salinity: a salinity bias in the surface layer of about 0.05 psu or more. Deeper in the water column, we clearly see on the in situ data the signature of Mediterranean deep water, with a much higher salinity. This salinity signature is very weak, although present in the model. We also see from this plot that typical water masses, namely intermediate and Mediterranean waters, are found at different depths in the observations and the model. These remarks make us think that the initial conditions from Mercator Psy2v1 suffer from errors and biases in the representation of the deep water, especially as far as the Mediterranean water mass is concerned.

These discrepancies have surprisingly little effect on the N(z) profiles, which have been calculated using the *neutral density gradient* method from Chelton et al. (1998). We see on Fig. 5c the upper layer N(z) profile from the same observation as in Fig. 5b. and the corresponding model output. The maxima of N(z) are of the same order of magnitude, between 20 and 25 10⁻³ s⁻¹, but the depth of the maximum is deeper in the observations. Below this high gradient at the base of the mixed layer, the minimum of N(z) from the observation is not well represented by the

model, as seen on Fig. 5d; but deeper in the water column, both profiles have comparable amplitude.

From the N(z) profiles, we can calculate the first Rossby radius using the formula from Chelton et al. (1998):

$$R_1 = \frac{1}{|f|\pi} \int_{-H}^{0} N(z) dz \qquad (10)$$

As the N(z) profiles from observations and model outputs are very close, corresponding first Rossby radii are very close, too. Table 1 gives a comparison of the calculated first Rossby radii for the available points in the zone during our study period. Discrepancies between observations and model outputs are of about 2%. The agreement between available observations and model outputs makes us reasonably confident in the realism of our simulation in terms of typical

<table><caption>**Table 1** Observed and modeled values of the internal radius of deformation R1 (km) at locations indicated as squares on Fig. 5a on the indicated dates</caption><thead><tr><td></td><td>Observed</td><td>Modeled</td></tr></thead><tbody><tr><td>Point a (07/09)</td><td>15.01</td><td>15.19</td></tr><tr><td>Point b (07/19)</td><td>15.11</td><td>15.45</td></tr><tr><td>Point c (07/29)</td><td>16.09</td><td>16.46</td></tr><tr><td>Point d (08/08)</td><td>15.39</td><td>15.41</td></tr><tr><td>Point e (08/18)</td><td>15.91</td><td>15.88</td></tr><tr><td>Point f (08/28)</td><td>15.84</td><td>16.12</td></tr></tbody></table>
horizontal scales of mesoscale features in the region covered by the floats, i.e., the southwest part of the Bay of Biscay abyssal plain.

We also check typical dynamical scales from our model by comparing sea surface variations with observations along altimeter tracks. These tracks have been processed using the X-track tool (Roblou et al. 2007), which allows getting valuable information closer to the coast as compared to usual gridded altimetric datasets and spatially filtered to remove very short scale features. Figure 6 shows Jason track 137, seen on the map on the left. That track crosses the open ocean and the shelf areas. The first pass on July 10 shows overall good agreement, probably due to the fact that the Mercator Psy2v1 model, from which our initial condition was taken, assimilates altimetry data. A similar level of agreement remains throughout the run, at least for the spatial scales of open ocean features. Major discrepancies appear however on top of the wide shelf in the north and close to the Spanish coast. Using 2-h model outputs in a test experiment (not shown) led to similar discrepancies, while comparisons of the model with tide gauge sea level, from which tides have been removed by harmonic analysis, showed fair agreement. Those discrepancies of altimetric data on continental shelves are not unexpected, partly because of the inadequate sampling of altimetry for fast dynamical processes and partly because of the specific limitations of altimetry in coastal areas (Anzenhofer et al. 1999). However, far from land, we found a fair agreement in the size and amplitude of the dominant sea surface elevation features, i.e., the type of dynamics we are interested in.

### 3.3 Ensemble strategy

The model is used to test various observation networks, ultimately using the technique described in Section 2. We used a Monte Carlo technique to generate the scaled representer matrix $\chi$ from ensemble statistics. The model was run several times (a few tens, depending on the case), each time randomly perturbing the same forcing parameter,

here the wind vector. The dispersion of the ensemble is considered, for a large enough number of members in the ensemble, as a useable proxy of the guess error. The guess error covariance matrix is approximated by:

$$ \mathbf{P}^g = E[\boldsymbol{\epsilon}^g \boldsymbol{\epsilon}^g \mathrm{T}] \approx (\mathbf{x}_i^g - \mathbf{x}_i^g)(\mathbf{x}_i^g - \mathbf{x}_i^g)^\mathrm{T} \quad (11) $$

where $\epsilon^g$ is the guess error and $\mathbf{x}_i^g$ is the “prior” model state for the ith member of the ensemble.

We decided to perturb the wind components, as the wind is a major driver of the dynamics of the ocean and, in particular, of our coastal current variability. Alternative sources of perturbations, considered in separate studies, are the open boundary conditions from the large-scale forcing model, the bathymetry, and the major river run-offs. The open boundary conditions have a major impact on the deep ocean dynamics as it drives the inflows and outflows of the various water masses, but a smaller one on the shelf, as the shallow water column there is more influenced by the atmosphere and the bottom stress. The bathymetry, on the opposite, would be more likely to affect the shelf than the open ocean, as it directly affects trapped coastal barotropic waves forced by tides and winds that are major dynamical features there. It would also affect barotropic dynamics in the open ocean, but not as predominantly as on the shelf. The river run-offs would predominantly affect the shelf water masses. Here, we only consider errors due to winds, as it is not the scope of this paper to test the various sources for errors.

In order to generate as realistic as possible wind perturbations, we first calculated the first ten bivariate empirical orthogonal functions of the zonal and meridional wind components in the modeling period and domain, using a Lanczos decomposition (Toumazou and Crétaux 2001):

$$ \mathbf{u}_{\mathrm{ref}}(x, y, t) - \overline{\mathbf{u}_{\mathrm{ref}}(x, y, t)} = \sum_{k=1}^{10} \lambda_k \mathbf{u}'_k(x, y) \zeta_k(t) \quad (12) $$

where the overbar denotes the time mean over the whole vector. $\mathbf{u}_{\mathrm{ref}}$ is the wind vector used for the unperturbed

<figure><img src="image_7.png" /><figcaption>Fig. 6 Comparison of local sea level variations (cm) along the track (km) seen on the map (200, 2,000, and 4,000 m isobaths shown), for the observations (thin lines with crosses) and the corresponding model outputs (thick lines with circles). The spatial mean along the track has</figcaption></figure>

been removed for both the observations and the model outputs. Observations are spatially smoothed by a Hanning filter with a 40-km frame
reference simulation, $\lambda_k$ is the kth eigenvalue of the decomposition, $u'_k(x,y)$ is the kth bivariate EOFs in space and $\zeta_k(t)$ is the kth EOF amplitude time series. The perturbed wind vector is then calculated as the sum of the original wind vector plus a linear combination of perturbations of the EOFs, which writes for any member m of the ensemble:

$$u_m(x, y, t) = u_{\text{ref}}(x, y, t) + \alpha \sum_{k=1}^{10} \delta_k^m \lambda_k u'_k(x, y) \zeta_k(t) \quad (13)$$

The method is similar to the one put forward by Auclair et al. (2003). We used a random Gaussian coefficient $\delta_k^m$ of zero mean and standard deviation of 1 for each term of the sum. The ten $\delta_k^m$ coefficients are changed every 5 days, corresponding to the peak of highest variability in the spectrum associated to the time series of the first EOF mode, in order to enlarge the variability of the ensemble. This period of 5 days is longer than typical low pressure events, which generally last 2 or 3 days, but allows the generated error patterns to develop in time. $\alpha$ represents the amount of uncertainty of the wind field. This factor has been set to 0.3, meaning 30% of uncertainty, slightly larger than the 20% uncertainty assumed for Aladin atmospheric fields (Loik Berre, Meteo-France personal communication). This overestimation of errors has been decided in order to artificially enhance ensemble dispersion, as the simulation is relatively short. Figure 7 shows the first three bivariate EOFs of the wind components with the associated time series and explained variance. The first two modes are dominant in explained variance. They both seem associated with the low-pressure systems of mid-August, since we see on their time series that they reach maximum amplitude at that time. The first mode represents cyclonic winds associated to low pressure over Great Britain, typical of the pressure system of the area. The second one is composed of Northward wind that is involved in the two low-pressure events of August. The third mode depicts strong cyclonic winds around Brittany and is clearly associated with the low-pressure event of the beginning of July. The time series show that this mode is barely involved

later on. Using combinations of EOFs ensures a realistic perturbation, since it will enhance or damp patterns composing the reference wind at any particular time.

Our ensemble had to be limited to 82 members, for practical reasons linked to the computer memory required for each simulation. However, we conducted a rough verification of the statistical robustness with respect to ensemble size. Figure 8 shows the dominant ensemble EOFs for instantaneous sea surface elevation, and their relative explained variance, for 50, 60, 70, and 80 members. The spatial patterns illustrated for the first three modes are similar for each tested ensemble size. The corresponding variances explained by the first three EOFs modes in sea surface elevation, explaining together about 86% of the total variance, are very close to each other. The relative robustness in the ensemble sea surface elevation statistics gives us some confidence in the relevance of our results, even though a larger number of simulations would probably be required to formally ensure the stability of the ensemble statistics.

## 3.4 The error subspace

As the RMS technique aims at assessing how an observational network is able to detect the model state error subspace, it is necessary to first get a good grasp of the dynamical processes present in that subspace. We want to know what physical processes are at work during the birth and growth of state error features and what the typical patterns of the error structures and their time evolution are.

Figure 9 illustrates the time evolution of the ensemble variability associated with the 82 model runs, in wind stress, sea surface elevation, and temperature on the vertical, as well as the reference wind stress at two locations shown on Fig. 5. Point 1 is on the Armorican Shelf; point 2 is in the middle of the basin. These points illustrate two types of model behavior in the error subspace associated to the wind forcing errors. The reference wind stress appears to be more intense on the open ocean (Fig. 9b) than closer

<figure><img src="image_8.png" /><figcaption>Fig. 7 a First EOF mode for the wind (amplitude and direction, no dimension) on the zone during the study period, with its associated time series (no dimension). Percentage of variance explained in the title. b, c Second and third modes. a-c Dashed lines Zero values</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8** *Left* First three EOFs modes in model sea surface elevation on July 29th and 30th, 2004, as a function of ensemble size (no dimension). *Right* Explained variance (%) by the first ten EOFs modes, for the four ensemble sizes tested</figcaption></figure>

to the coast (Fig. 9a). The wind stress ensemble standard
deviation, on Figs. 9c and d, follows the wind stress
amplitude itself, due to the perturbation protocol. The error
regime on the shelf can be characterized by short-term, high-
amplitude processes forced by meteorological events with
rather constant “background” error, whereas in the open
ocean, the error regime shows lower amplitudes, rising with
each new meteorological event. On the shelf, the sea level is

highly sensitive to winds that excite high-amplitude coastal
trapped waves in the error subspace (Fig. 9e). These
barotropic waves (Carrère and Lyard 2003) have a short life
span of a few hours to a few days, which is the time scale of
the ensemble variance peaks. In temperature, on Fig. 9g, the
maximum ensemble variance is located just below the mixed
layer, as seen on the isotherm contours. As the thermocline is
thinner on the shelf, perturbations in winds that affect

<figure><img src="image_10.png" /><figcaption>**Fig. 9** *Left* Model outputs from point 1 (cf. Fig. 5). *Right* Same for point 2. *a, b* Evolution of the reference wind stress intensity (N m⁻²) during the simulation. *c, d* Wind stress intensity ensemble standard deviation (N m⁻²). *e, f* Sea surface elevation ensemble standard deviation (cm). *g, h* Upper 100 m temperature ensemble standard deviation (°C) on the vertical (y-axis, in m). *Thin white lines* Reference model potential temperature (contour every 1°C). Bottom is indicated (*black line*) on *g*. For all: reference dates from Fig. 4 in *dashed lines*</figcaption></figure>
vertical mixing lead to large temperature ensemble variance at the thermocline depth. In the mixed layer, temperature variability is large after a major wind event but is otherwise rather constant, as is the sea level variance. In the open ocean, the sea level and temperature ensemble variances are overall smoother than on the shelf. The sea level ensemble variance (Fig. 9f) reflects the first blast of wind at the beginning of July and rises at each new wind event during the simulation. The depth of the maximum temperature ensemble variance (Fig. 9h) is located at the thermocline as it was on the shelf. Its amplitude is lower, probably because the vertical temperature gradient is weaker. Closer to the surface, the ensemble variance rises at each new high wind event as does the sea level ensemble variance and affects the whole mixed layer.

As we are focusing on mesoscale activity, mostly in relation with coastal currents, we now check that these features are visible in the model state error subspace. Figure 10a shows the vortical structure at 30-m depth in the reference model run on August 6. We clearly see major mesoscale patterns in the open ocean area and at the shelf edge. Two major cyclonic eddies neighbor the Armorican Shelf edge, forced as they are by the shear between the northward slope current and the southward interior current. A large cyclonic meander spreads inwards from the Spanish coast at 5.5° W. Small cyclonic patches appear to be stuck between this latter structure and the two major eddies along the Armorican Shelf slope. Close to the model western boundary, two small eddies are found, a cyclonic one in the North at the shelf edge and an anticyclonic one in the South, along the Cantabrian shelf. In the Southeast corner of the Bay of Biscay, a mushroom-like meander is developing on the Landes Plateau toward the interior, with a cyclonic part in the North and the anticyclonic counterpart in the South in which

currents are following the curve of the shelf slope. Along the shelf edge and on the shelf itself, cyclonic and anticyclonic features have much smaller scale and appear to be less coherent during the simulation. Figure 10b shows the ensemble potential temperature standard deviation at the same depth and on the same day. The largest amplitude for temperature ensemble variance at this depth is found on the Cantabrian shelf and slope, on the Landes Plateau and the adjacent shelf, and at the edge of the Brittany peninsula. Temperature ensemble variance on the shelf can be locally very large when reaching the bottom of the mixed layer where the temperature vertical gradient is intense. On the Landes Plateau, the large amplitude temperature ensemble variance meander follows the mushroom-like vorticity pattern seen on Fig. 10a, and the anticyclonic part of that structure, North of it, extends inside the Cap-Breton Canyon within the anticyclonic part. In the interior of the basin, the temperature ensemble variance is larger at the edges of dominant vortical features. It is particularly true around the cyclonic meander spreading from the Spanish coast, the shape of which appears clearly in temperature ensemble variance and around the anticyclonic patch north of it. Meanders of temperature ensemble variance also follow vortical dynamics west of these features. Temperature ensemble variance thus appears intensely influenced by mesoscale activity in the error subspace in that region. Figure 10c is a map of the instantaneous sea level ensemble standard deviation on the same day. Sea level ensemble variance is higher in shallow waters than in the open ocean as expected. In the interior of the basin, the ensemble variance partly follows mesoscale features, like the meander spreading from the Spanish coast. Two patches of high-amplitude ensemble variance are present at the edge of the Landes Plateau and in the Cap-Breton Canyon north of it,

**Fig. 10** **a** Model relative vorticity ($\partial\nu/\partial y - \partial u/\partial x$, $10^{-4} \text{ s}^{-1}$) at 30-m depth on August 6.

**b** Model potential temperature ensemble standard deviation (°C) at 30-m depth on the same day.

The zonal line at 44.8° N is the location for transect used in **d**.

**c** Model sea surface elevation ensemble standard deviation (cm) on the same day. **a–c Thin lines** 200, 2,000, and 4,000 m isobaths. **d Vertical transect** along 44.8° N of model potential temperature ensemble standard deviation (°C) in the upper layer on the same day. **Thin white lines** Reference model potential temperature on the same day (contour every 1°C).

**Thin black line** 30-m depth reference

<figure><img src="image_11.png" /></figure>
which are places of high mesoscale activity. Smaller patches
are present along the Armorican Shelf slope in the vicinity of
the two large cyclonic eddies. Figure 10d shows the
ensemble standard deviation of model temperature on the
vertical, along the transect shown on Fig. 10b, at the same
date. In the open ocean, from 10.5° W to 4° W, the depth of
the temperature ensemble variance maximum and its
amplitude appear to be correlated to the depth and amplitude
of the maximum vertical temperature gradient, as inferred
from the isotherms. The amplitude of the temperature
ensemble standard deviation rises when approaching the
Landes Plateau, around 4° W. There, it strengthens and
deepens. Closer to the coast, the temperature ensemble
variance seems more complex with smaller variance patches
at several depths. Short space scale variations of the isotherm
depths there indicate strong local dynamics, which undoubt-
edly play a major role in the temperature error pattern found.
This illustrates the complexity of the Landes Plateau zone.

3.5 Comparison of altimetry networks

We now use the previously described approach (RMS
technique and stochastic modeling) to estimate the perfor-
mance of a few observational networks. The Appendix
gives a brief overview on how the calculations are handled
in practice, as well as some cost-related considerations.

We first compare the performances of two altimeters.
One is the classical Jason 1 altimeter; the other is a wide
swath altimeter, such as the instrument envisaged for the
Surface Water Ocean Topography (SWOT) projected
mission, inspired by the Wide Swath Ocean Altimeter
(WSOA) project (Alsdorf et al. 2007). This “imaging”
altimeter would be composed of a classical nadir radar

altimeter, plus an interferometer measuring sea surface
elevation in the cross-track direction and therefore provid-
ing much wider space coverage than a classical radar, nadir-
only altimeter. We are interested not only in comparing the
relative efficiency of each instrument at detecting model
variability modes but also in estimating what physical
processes these instruments are able to measure. Two
experiments are performed with the wide-swath instrument.
The first experiment accounts for all the available measure-
ments points; the second accounts for only those which are
over the deep ocean (more than 150 m depth). The second
experiment is designed to test the influence of the space
locations of altimetric data. In Fig. 11, we see the
observation points of the various networks. Jason 1
measurement locations have been taken from real observa-
tions, with a datum every 5.8 km on average, corresponding
to the 1 Hz data post-processing. The tracks of the wide
swath altimeter are centered on those of Jason 1, which is
one of the realistic options considered for the SWOT
mission. The measurement points are 15 km distant in both
along-track and cross-track directions, as proposed in the
WSOA mission definition (Fu 2003). The error budget is,
in both cases, rather simple. Observation errors are
considered independent in space and time, equal to a
Gaussian white noise of zero mean. The observation error
covariance matrices **R** are thus diagonal. In the case of
Jason 1, the standard deviation is taken equal to 3.9 cm.
This value is larger than the 2-cm error achieved for large-
scale variability in datasets available on basin-wide areas
(Le Traon and Ogor 1998). The same error of 3.9 cm in sea
level for coastal ocean has already been used by Mourre et
al. (2006) and is relevant for the achieved error along the
satellite tracks using improved data processing in coastal

<figure><img src="image_12.png" /><figcaption>**Fig. 11** *a* Observation locations and dates for the Jason 1 experiment, with 200, 2,000 and 4,000 m isobaths in *thin black lines*. *b* Same as *a* for the wide swath experiment. *c* Same as *a* for the wide swath deep ocean-only experiment. *a–c* Dates for each measurement are in color levels. *d* Representer matrices spectra for the Jason 1 experiment (*thin line with crosses*), the wide swath experiment (*thick line with circles*), and the wide swath deep ocean-only experiment (*dot-dashed line with squares*). *Dotted line* Spectrum of the identity matrix, equal to 1</figcaption></figure>
areas (Bouffard et al. 2008). The wide swath altimeter error budget depends on the distance to the nadir. On the nadir, the error is one of the radar altimeter and is then set to 3.9 cm as well, then it increases from 3.9 cm at nadir to 5.3 cm at the far end of the swath (Fu 2003). This does not take into account other sources of errors, such as the satellite roll that induces correlated errors along the satellite track [the impact of such errors has been studied by Le Hénaff et al. (2008) in the case of barotropic shelf dynamics]. The observational sampling period spans the period from July 27 to August 5, corresponding to a full Jason 1 orbit and does not include a storm (thus simplifying error propagation analysis). Figure 11a–c shows the observation points and the associated days for the three networks. Jason 1 has 732 observations, the full wide swath altimeter has 3,910, and the wide swath altimeter over deep ocean only has 2,779. For each network, the scaled representer matrix χ has been calculated, then diagonalized using a Lanczos method (Toumazou and Crétaux 2001). Figure 11d provides the spectrum of each network χ matrix. Considering the number of eigenvalues larger than 1, as defined in Section 2, the most efficient network is the full wide swath altimeter, which detects 4 *df* of model state error, then the wide swath deep ocean-only altimeter, detecting 2*df*, and finally Jason 1 of only 1 *df*.

By itself, the spectral analysis does not tell which model state error physics are gained or lost when choosing one

network over another. To access that information, we have to examine the corresponding eigenvectors. Figure 12 shows the corresponding dominant 3 eigenvectors of χ for each altimetric network. The first mode is a basin mode consisting of an exchange of water between the open ocean surface elevation and the shelf surface elevation. This large-scale mode appears to be detected by the three networks. The second and third modes, seen clearly on the wide swath eigenvectors, have shorter space scales, closer to the order of magnitude of open ocean mesoscale features. The second mode exhibits positive patterns of large amplitude on the edge of the Landes Plateau, north of Spain at 5° W and at the northern edge of the Armorican Shelf, where the shelf slope enlarges. Negative patterns appear at the extreme northwest of the domain and in a small patch at the base of the shelf slope just north of the Landes Plateau. Mode 2 eigenvectors for all three networks have similar locations for their relative maxima and minima, considering the differences in the data sampling density. The third mode shows high-amplitude patterns of both signs mainly on the shelf, a major positive signals on the Landes Plateau, and smaller positive signals along the Cantabrian slope. The high-amplitude, double-signed feature on the shelf lets us think of short time-scale error-space events being detected by the wide-swath network there. Modes 2 and 3 from the full wide swath altimeter show many of the mesoscale sea level features of the interior of the basin as seen on Fig. 10:

<figure><img src="image_13.png" /><figcaption>**Fig. 12** First three dominant eigenmodes (no dimension) of the scaled representer matrices, for the three altimetry experiments from Fig. 11: Jason 1 (top), full wide swath (middle), and white swath deep-ocean only (bottom). *From left to right* First, second, and third modes. *Thick line* 200, 2,000, and 4,000 m isobaths</figcaption></figure>
the meander North of Spain, the patches on the Landes Plateau and in the Cap-Breton Canyon, and the same patches along the Armorican Shelf.

The combined examination of the number of error modes detected by a network and of the patterns associated to these modes allows characterizing the performances of the various networks in a qualitative way, considering the physics associated to the various modes. The previous discussion seems to indicate that the dominant three modes are related across networks in terms of those physical processes. Jason 1 can detect the water sloshing mode but not the higher error-space modes associated with the mesoscale variability. A wide swath altimeter appears to be able to detect that last type of error-space variability. Another interesting feature is the opposite-sign swaths on the Armorican shelf on wide-swath mode 3. Both swaths were measured at different times and led to the detection of error features with opposite signs, hinting at the presence of fast error-space processes on the shelf. These two patches of opposite signs reveal an underlying error-space anti-correlation, on the shelf, between the two swath passes, which occurred at different times. This is consistent with high-frequency, wind-forced shelf dynamics being present in the error subspace there. In order to enhance separability between the high-frequency shelf processes and the abyssal plain mesoscale dynamics over such a short time span, we introduced the third network, which only covers the deepest areas with the wide swath. With that configuration, the wide swath instrument is able to detect only two modes, the second mode being as before associated with mesoscale features. Overall, a wide swath altimeter seems to be able to detect error processes linked to mesoscale activity in a coastal zone, which is a requirement from the coastal modeling community.

This example shows that the RMS method is not to be used as a black box by just counting eigenvalues. The

definition of the observation network has to be designed, keeping in mind what the model is able to represent and what error processes we are interested in detecting.

## 3.6 In situ array design

The RMS technique may also be used to study several options for a cruise or for the deployment of a local ad hoc array. To illustrate this, we now examine various arrays of in situ measurements in one particular subregion of the Bay of Biscay. Several array options are being considered to study the mean currents and variability of the Landes Plateau, which is a zone of complex dynamics with high mesoscale activity on top of an irregular sloping bottom. As before, we wish to know whether the array options considered would help us detect local (and possibly remote) model state errors and whether we can hierarchize the options given the criteria illustrated in this paper. Three networks are tested, as shown on Fig. 13a. The first two networks, A and B, only measure temperature along the water column every 10 m down to 2,000 m (e.g., coarsely mimicking CTD or glider measurements). Network A, composed of 326 temperature measurements on 12 vertical profiles, is centered at the edge of the shelf along the Landes Plateau, while network B, with 316 measurements again on 12 vertical profiles, is shifted 35 km to the south. It is closer to the Cap-Breton Canyon and reaches the shelf south of it. Finally, network C is located as A, but current meters are added along the central transect. Counting the additional velocity measurements, network C is composed of 538 measurements, on the same profiles as network A. Our temporal sampling scheme is simplified in this example: we consider observations to be concomitant on July 30 only. The observation errors are once again chosen to be independent. Observation error standard deviations have been set to 0.2°C in temperature, 1 cm s⁻¹ in zonal and meridional

<figure><img src="image_14.png" /><figcaption>**Fig. 13** a Location of the observations in network A (white squares), network B (black circles), and network C (same as A, dots in central squares indicate current-meter locations). b Representer matrices</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>spectra for network A (thin line with crosses), network B (thick line with circles), and network C (dot-dashed line with squares). Dotted line Spectrum of the identity matrix, equal to 1</figcaption></figure>
<table><caption>**Table 2** Total variance of the scaled errors (Global Detected Error Variance, GDEV) for the three networks (no dimension) (first column) and part of the variance explained by the first four modes (%) (columns 2 to 5)</caption><thead><tr><th></th><th>Total scaled errors variance</th><th>Mode 1</th><th>Mode 2</th><th>Mode 3</th><th>Mode 4</th></tr></thead><tbody><tr><td>Network A</td><td>427</td><td>76.2</td><td>8.9</td><td>4.3</td><td>2.3</td></tr><tr><td>Network B</td><td>555</td><td>77.3</td><td>8.5</td><td>5.5</td><td>2.2</td></tr><tr><td>Network C</td><td>530</td><td>65.3</td><td>12.4</td><td>7.5</td><td>5</td></tr></tbody></table>

currents, for any measurement point. **R** matrices are thus diagonal, homogeneous over each measurement type.

Spectra for the scaled representer matrices $\chi$ are shown on Fig. 13b. As a first remark, there are much more modes with eigenvalues larger than one than in the previous example. This is due to the high-amplitude model guess errors in temperature, especially at the base of the mixed layer, which are twice as large as the observational error. In contrast, errors on sea surface elevation in the altimetric examples were generally smaller than the observational error, at least on the open ocean. Spectra of networks A and B appear to be very close to each other for the dominant modes, though one must be careful since the scale is logarithmic. Table 2 shows that the global detected error variance (GDEV), measured by the total scaled error variance and equal to the trace of $\chi$, is much larger for network B because of the contribution of the first eigenvalue. These two spectra also differ in that network B has 14 eigenvalues larger than one, one more than network A. Shifting the array network helps better detect the dominant mode and enriches the range of model guess error *dfs* the network is able to detect. Network C has 15

eigenvalues larger than one and appears to be more efficient than networks A and B for measuring the dominant modes other than the first, for which network B is the most efficient. Table 2 indicates that network C overall measures less GDEV, whereas its entire spectrum is above network B's, except for the first eigenvalue. Network C has relative explained variances for modes 2, 3, and 4 higher than the other two networks and a relatively less dominant first mode. Adding current meters thus improves the ability of the network to detect a palette of secondary modes.

Another point of interest is to compare the impact each network would have in state space, i.e., on the model fields if such a network was assimilated, at least for the dominant modes. This can be achieved through the study of the modal representers defined in Section 2. These state-space objects are more convenient here than the sole eigenvectors of $\chi$, as groups of vertical measurements are difficult to analyze together. Modal representers of the three networks in terms of model sea level are shown on Fig. 14 for the first three modes. The first-mode modal representer patterns are very similar for the three networks. We notice on this mode the large sea level pattern already seen in sea level

<figure><img src="image_16.png" /><figcaption>**Fig. 14** First three modal representers in sea level (cm), for the three in situ networks from Fig. 13: network A (top), network B (middle), and network C (bottom). From left to right First, second, and third modes. Thick dots 200, 2,000, and 4,000 m isobaths</figcaption></figure>
ensemble standard deviation (Fig. 10c), which was also seen in mode 3 of the altimetry networks analysis (Fig. 12). In this particular mode, network B appears to show the largest amplitude, while the other two seem rather equiv- alent. This network has a larger potential to impact the French coast model sea level due to its position closer to the shelf and reaching the along coast variability patch of mode 1. Modes 2 and 3 differ across the networks studied. Indeed, mode 2 for network C has the same pattern as mode 3 for networks A and B. This modal representer has larger amplitude for network C: adding current meters to the initial array would give network C the potential to correct the model sea level along the coast and at the western edge of the Landes Plateau, even more so than shifting the array on the Cap-Breton Canyon. This ability had not been found for the first mode, for which networks A and C had comparable eigenvalues and modal representer amplitudes. Similarly, network B modal representer has a larger amplitude than A: once again, shifting the array to the south of the Plateau allows a larger impact on the model sea level along the French coast. For the remaining modes, it is not possible to assert that they detect the same error mode.

Here again, a wealth of useful information can be
obtained beyond the mere representation of Representer
Matrix Spectra, thanks to the other tools introduced in this
paper. In situ temperature data appear to provide a
potentially strong constraint on the model sea level, through
the covariance of model state errors. Shifting the array
South on the Cap-Breton Canyon allows reaching a region
of higher model sea level variability, thus better impacting
the sea level along the coast. Adding current meters to the
array modifies the order of the error modes caught by the
network and enhances the capability to detect secondary
modes.

4 Conclusion

A relatively simple methodology to assess the performance
of any observational network at detecting model state errors
without having to perform data assimilation has been
developed and tested on realistic cases. Given a model
prior state estimate, our technique, called the Representer
Matrix Spectra (RMS) technique, aims at analyzing how
observation networks are potentially able to detect and
constrain degrees of freedom of the model state errors. In
this paper, the methodology was implemented with the help
of stochastic modeling in order to provide the required error
estimates. Using that approach in the case of wind forcing
errors, we have illustrated the potential of wide swath
altimetry to detect errors associated with mesoscale vari-
ability in the deepest part of the Bay of Biscay and with the
slope current variability along the shelf break, as well as

with high-frequency wind error-induced processes on the
shelf. In a second example, we focused on instantaneous
measurements on the Landes Plateau. As the prior model
state errors are very large in temperature on the vertical,
measurements appear very valuable to the model. Shifting
the array closer to the south edge of the Landes Plateau
allowed measuring more of the along-coast variability,
while adding current measurements better detects secondary
variability modes.

The configuration adopted in these two examples was
based on a single ensemble implementing a unique source
of state errors, here those of the forcing winds. The way the
perturbed wind field is generated, using a set of random
coefficients every 5 days, could be refined by using a
smoother temporal scheme. Other sources could be consid-
ered, for example, the model initial state. Alternative
ensembles, mixed ensembles, or ensembles drawn from
data assimilation schemes, such as the EnKF, would
probably be needed to fully assess network performances
in realistic cases.

These two examples illustrate the potential of the method
in both quantitative and qualitative ways. The RMS
technique is not able to define an optimum observation
network by itself, but it gives a way to assess and compare
array options at low computational cost, as one has to
calculate only one ensemble of model runs and then
perform the various analyses. Moreover, the analysis of
array modes and modal representers gives valuable insight
into the impact a network would have on the model state
error subspace physics, thus helping choose a network that
would focus on specified physics. This qualitative part of
the RMS technique is undoubtedly the most interesting and
promising.

Results from the RMS technique are indicative of the
potential impact an observation network may have in a data
assimilation context. We are thus interested in confronting
our approach with other observation network assessment
techniques, such as OSSEs involving the EnKF (Mourre et
al. 2006). This will be one of the next steps of this work.
We also wish to test the use of on-the-fly RMS analysis
during EnKF integration, with the objectives of studying
the regime dependency of the array mode spectra and of
confronting observability criteria (such as in this paper) to
controllability criteria. Another point of interest would be to
test the RMS method with more complex observation error
budgets, involving correlation in observation errors, which
may affect the perception of the model state error modes by
the array. Finally, this approach should be complemented
by a model-testing approach once the “optimized” network
data are finally available.

**Acknowledgments** We thank MERCATOR, a GODAE-participating system, for providing the model boundary data used in the nesting.
Most of the work was conducted during Matthieu Le Hénaff's Ph.D. thesis, for which a grant from DGA is gratefully acknowledged. This project also benefited from interactions within the MERSEA and ECOOP European projects, and (in the case of nadir and wide-swath altimetry) within the CNES/NASA OST and SWOT communities. We also thank two anonymous reviewers for their constructive remarks.

## Appendix: practical calculations and cost considerations

The ensemble of 82 simulations provides us with a proxy of model guess error statistics. Manipulating such an ensemble of high-resolution, 3D simulations is not straightforward. Our model state has a full dimension (space, time, and variables) of $n = O(10^9)$, which is large. Calculating the various terms explained in Section 2.1 can thus rather be heavy depending on how this is done. Here, the matrix $\chi$ is calculated in three steps. First, we directly calculate $\textbf{HP}^g\textbf{H}^T$, and not $\textbf{P}^g$, which is too costly to determine, and too heavy to store, as its size is $O(n^2)$. The projection onto the observation space through $\textbf{H}$ is a spatial trilinear interpolation (horizontal and vertical) of the model values on the observation points, for each member; we consider in our study that the observation times correspond to some of the model output dates, as each model output is already a 2-day average. If $p$ is the dimension of the observation space and $N$ the size of the ensemble, the computation cost to calculate the representer matrix $\textbf{HP}^g\textbf{H}^T$ is $O(p^2N)$. The second step is the calculation of $\textbf{R}^{-1/2}$. It is straightforward when we consider uncorrelated errors: $\textbf{R}$ is diagonal, so $\textbf{R}^{-1/2}$ is also diagonal, with values equal to the inverse of the square root of the diagonal terms of $\textbf{R}$. When we consider a more complex error budget (not the case in our tests), the computation cost of $\textbf{R}^{-1/2}$ increases. The calculation of $\textbf{R}^{-1}$ can be performed for instance using a Cholesky decomposition, e.g., using Lapack libraries (DPPTRF/DPPTRI functions), since $\textbf{R}$ is symmetric definite positive by construction. A Cholesky decomposition of $\textbf{R}^{-1}$, which is symmetric definite positive too, can also be used to calculate $\textbf{R}^{-1/2}$. Using $\textbf{R}^{-1} = \textbf{LL}^*$ ($\textbf{L}$ indicating the conjugate) and the SVD decomposition of $\textbf{L}$ as $\textbf{L} = \textbf{USV}^*$, then $\textbf{LL}^* = (\textbf{USU}^*)(\textbf{USU}^*) = \textbf{R}^{-1}$, so $\textbf{R}^{-1/2} = \textbf{USU}^*$. Both Cholesky and SVD decompositions cost $O(p^3)$ calculations, which is more costly than the calculation of $\textbf{HP}^g\textbf{H}^T$ when $p > N$, which is the general case. The last step is the calculation of $\chi$ from $\textbf{HP}^g\textbf{H}^T$ and $\textbf{R}^{-1/2}$; it requires $O(p^2)$ operations when $\textbf{R}^{-1/2}$ is diagonal, $O(p^3)$ when it is not, which is of the same order as the cost of the calculation of $\textbf{R}^{-1/2}$ itself. The final eigenvalue decomposition of $\chi$ requires $O(p^2)$ operations. Overall, in our case where $\textbf{R}^{-1/2}$ is diagonal, the most expensive part is the calculation of $\textbf{HP}^g\textbf{H}^T$, with a cost of $O(p^2N)$ calculations.

The matrix of modal representers $\rho = \textbf{P}^g\textbf{H}^T\textbf{R}^{-1/2}\mu$ is also costly to compute. As we are interested in the

dominant modes of $\chi$, we select an arbitrary number $q \ll p$ of modes we wish to examine, using a block $\mu'$ of dimension $p*q$. If $\textbf{R}^{-1/2}$ is diagonal, the final matrix $\rho'$ is calculated in one step, requiring $O(N\textit{npq})$ calculations, which is larger than calculating $\textbf{P}^g\textbf{H}^T$ separately in $O(N\textit{np})$ calculations, but avoiding having to store this matrix of large dimension $n*p$. If $\textbf{R}^{-1/2}$ is not diagonal, calculating $\rho'$ in one single step requires $O(N\textit{np}^2q)$ operations, which is very large, so one may prefer to calculate and store $\textbf{P}^g\textbf{H}^T$ and $\textbf{R}^{-1/2}\mu'$ separately. $\textbf{P}^g\textbf{H}^T$ costs $O(N\textit{np})$ calculations and has a size of $n*p$. This induces a limit in the size of the observation space: if $p = 10^3$ observations, $\textbf{P}^g\textbf{H}^T$ has a size of $O(10^{12})! \textbf{R}^{-1/2}\mu'$, of dimension $p*q$, costs $O(p^2q)$ operations, which is expectedly lower than the calculation of $\textbf{P}^g\textbf{H}^T$, since $n > p$ and in the general case, $N > q$. The final calculation of $\rho'$ from the two matrices then requires $O(np^2q)$ calculations, which is the most expensive step in terms of calculations. In this study, we only calculate $\rho'$ in the case of diagonal $\textbf{R}$ matrices, in a single step of $O(N\textit{npq})$ calculations.

## References

* Alsdorf DE, Rodríguez E, Lettenmaier DP (2007) Measuring surface water from space. Rev Geophys 45:RG2002, doi:10.1029/2006RG000197

* Anzenhofer M, Shum CK, Rentsh M (1999) Coastal altimetry and applications. Tech. Rep. n. 464, Geodetic Science and Surveying, The Ohio State University Columbus, USA

* Arhan M, Colin de Verdière A, Mémery L (1994) The eastern boundary of the subtropical North Atlantic. J Phys Oceanogr 24:1295–1316

* Auclair F, Marsaleix P, De Mey P (2003) Space-time structure and dynamics of the forecast error in a coastal circulation model of the Gulf of Lions. Dyn Atmos Ocean 36:309–346

* Auclair F, Estournel C, Marsaleix P, Pairaud I (2006) On coastal ocean embedded modeling. Geophys Res Lett 33:L14602, doi:10.1029/2006GL026099

* Bardey P, Garnesson P, Moussu G, Wald L (1999) Joint analysis of temperature and ocean colour satellite images for mesoscale activities in the Gulf of Biscay. Int J Remote Sensing 20:1329–1341

* Barth N (1991) Oceanic experiment design II: genetic algorithms. J Atmos Oceanic Technol 9:434–443

* Barth N, Wunsch C (1990) Oceanographic experiment design by simulated annealing. J Phys Oceanogr 20:1249–1263

* Bennett AF (1985) Array Design by Inverse Methods. Prog Oceanogr 15:129–151

* Bennett AF (1990) Inverse methods for assessing ship-of-opportunity networks and estimating circulation and winds from tropical expendable bathythermograph data. J Geophys Res 95:16111–16148

* Bouffard J, Vignudelli S, Cipollini P, Menard Y (2008) Exploiting the potential of an improved multimission altimetric data set over the coastal ocean. Geophys Res Lett 35:L10601, doi:10.1029/2008GL033488

* Carrère L, Lyard F (2003) Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing—comparisons with observations. Geophys Res Lett 30(6):1275, doi:10.1029/2002GL016473
Chelton DB, DeSzoeke RA, Schlax MG, Naggar KE, Siwertz N (1998) Geographical variability of the first baroclinic Rossby radius of deformation. J Phys Oceanogr 28:433–460

Evensen G (1994) Sequential data assimilation with a nonlinear quasi-
geostrophic model using Monte Carlo methods to forecast error
statistics. J Geophys Res 99(C5):10143–10162

Fu L (2003) Wide-swath altimetric measurement of ocean surface topography. Report Jet Propul. Lab, Pasadena, CA, 67

Garcia-Soto C, Pingree RD, Valdés L (2002) Navidad development in the southern Bay of Biscay: Climate change and swoddy structure from remote sensing and in situ measurements. J Geophys Res 107(C8):3118, doi:10.1029/2001JC001012

Guinehut S, Larnicol G, Le Traon P-Y (2002) Design of an array of profiling floats in the North Atlantic from model simulations. J Mar Syst 35:1–9

Hackert EC, Miller RN, Busalacchi AJ (1998) An optimized design for a moored instrument array in the tropical Atlantic Ocean. J Geophys Res 103(C4):7491–7510

Hernandez F, Le Traon P-Y, Barth NH (1995) Optimizing a drifter cast strategy with a Genetic Algorithm. J Atmos Oceanic Technol 12:330–345

Kurapov AL, Allen JS, Egbert GD, Miller RN, Kosro PM, Levine MD, Boyd T, Barth JA (2005) Assimilation of moored velocity data in a model of coastal wind-driven circulation off Oregon: Multivariate capabilities. J Geophys Res 110:C10S08, doi:10.1029/2004JC002493

Le Hénaff M, De Mey P, Mourre B, Le Traon P-Y (2008) Contribution of a wide-swath altimeter in a shelf seas assimilation system – Impact of the satellite roll errors. J Atmos Oceanic Technol (in press)

Le Traon P-Y, Ogor F (1998) ERS-1/2 orbit improvement using TOPEX/POSEIDON: The 2 cm challenge. J Geophys Res 103 (C4):8045–8058

Marsaleix P, Auclair F, Floor JW, Herrmann MJ, Estournel C, Pairaud I, Ulses C (2008) Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Model 20:61–89

McIntosh PC (1987) Systematic Design of Observational Arrays. J Phys Oceanogr 17:885–902

Mourre B, De Mey P, Menard Y, Lyard F, Le Provost C (2006) Relative performance of future altimeter systems and tide gauges in constraining a model of North Sea high-frequency barotropic dynamics. Ocean Dynamics 56:473–486

Pingree RD (1993) Flow of surface waters to the west of the British Isles and in the Bay of Biscay. Deep-Sea Res 40B:369–388

Pingree RD, Le Cann B (1989) Celtic and Armorican slope and shelf residual currents. Prog Oceanogr 23:303–338

Pingree RD, Le Cann B (1990) Structure, strength and seasonality of the slope currents in the vicinity of the Bay of Biscay region. J Mar Biol Assoc UK 70:857–885

Pingree R, Le Cann B (1992a) Three anticyclonic slope water oceanic eddies (SWODDIES) in the southern Bay of Biscay in 1990. Deep-Sea Res 39:1147–1175

Pingree R, Le Cann B (1992b) Anticyclonic eddy X91 in the southern Bay of Biscay, May 1991 to February 1992. J Geophys Res 97:14353–14367

Raicich F (2006) The assessment of temperature and salinity sampling strategies in the Mediterranean Sea: idealized and real cases. Ocean Sci 2:97–112

Roblou L, Lyard F, Le Hénaff M, Maraldi C (2007) X-track, a new processing tool for altimetry in coastal oceans, *ESA ENVISAT Symposium*, Montreux, Switzerland, April 23–27, 2007, ESA SP-636, July 2007

Schiller A, Wijffels SE, Meyers GA (2004) Design requirements for an Argo float array in the Indian Ocean inferred from observing system simulation experiments. J Atmos Oceanic Technol 21:1598–1620

Toumazou V, Crétaux JF (2001) Using a Langczos eigensolver in the computation of empirical orthogonal functions. Mon Weather Rev 129:1243–1250

van Aken HM (2000) The hydrography of the mid-latitude Northeast Atlantic Ocean: II, The intermediate water masses. Deep-Sea Res 47A:789–824

van Aken HM (2001) The hydrography of the mid-latitude Northeast Atlantic Ocean: Part III. The subducted thermocline water mass. Deep-Sea Res 48A:237–267

van Aken HM (2002) Surface currents in the Bay of Biscay as observed with drifters between 1995 and 1999. Deep-Sea Res 49A:1071–1086

Vignudelli S, Cipollini P, Roblou L, Lyard F, Gasparini GP, Manzella GRM, Astraldi M (2005) Improved satellite altimetry in coastal systems: Case study of the Corsica Channel (Mediterranean Sea). Geophys Res Let 32:L07608, doi:1029/2005GL22602