Depth dependence of westward-propagating North Atlantic features diagnosed from altimetry and a numerical 1/6° model

. A 1/6 ◦ numerical simulation is used to investigate the vertical structure of westward propagation between 1993 and 2000 in the North Atlantic ocean. The realism of the simulated westward propagating signals, interpreted princi-pally as the signature of ﬁrst-mode baroclinic Rossby waves (RW), is ﬁrst assessed by comparing the simulated amplitude and zonal phase speeds of Sea Level Anomalies (SLA) against TOPEX/Poseidon-ERS satellite altimeter data. Then, the (unobserved) subsurface signature of RW phase speeds is investigated from model outputs by means of the Radon Transform which was speciﬁcally adapted to focus on ﬁrst-mode baroclinic RW. The analysis is performed on observed and simulated SLA and along 9 simulated isopycnal displacements spanning the 0–3250 m depth range. Simulated RW phase speeds agree well with their observed counterparts at the surface, although with a slight slow bias. Below the surface, the simulated phase speeds exhibit a sys-tematic deceleration with increasing depth, by a factor that appears to vary geographically. Thus, while the reduction factor is about 15–18% on average at 3250 m over the region considered, it appears to be much weaker (


Introduction
Rossby or planetary waves (Rossby, 1939), discussed for the first time by Hough (1897) are instrumental in the western intensification of oceanic gyres (Anderson and Gill, 1975;Anderson and Killworth, 1977), in the large scale adjustment of ocean basins (Pedlosky, 1979;Gill, 1982), and possibly in the variability of the Meridional Overturning Circulation (MOC) (Hirschi et al., 2007). Until recently, the main theoretical framework to describe such waves has been the standard linear theory (referred to as SLT thereafter), which is based on quasi-geostrophic (QG) assumptions, a flat bottom hypothesis, and a decomposition of oceanic motions on a discrete basis of independent normal mode solutions (Gill, 1982). These include a very fast barotropic mode whose associated horizontal velocity profile is independent of depth, and an infinite number of baroclinic modes of increasing vertical complexity. The appealing simplicity of the SLT has made it a widely used interpretative tool by oceanographers, even though the validity, relevance, and accuracy of the dynamical modes it predicts have never been really tested owing to the traditional difficulty of sampling the oceans at the temporal and spatial scales pertaining to planetary wave propagation.
Over the last decade, however, the advent of satellite altimetry has rapidly led to improving our knowledge of the empirical properties of planetary wave propagation by providing global observations of the surface signature of these waves with unprecedented spatial and temporal resolutions. So far, although limited to the surface, such observations have already proved sufficient to question the validity of the SLT on a number of points. Specifically, Chelton and Schlax (1996) showed that observed phase speeds appear to be two to three time faster at mid-and high-latitudes than predicted by the SLT; furthermore, the observed wave amplitude do not appear to remain uniform throughout ocean basins as Published by Copernicus Publications on behalf of the European Geosciences Union. anticipated by the SLT, but to have enhanced variance in the western part of the basins. Chelton and Schlax (1996)'s results greatly stimulated the development of new theories seeking to improve upon the SLT by taking into account such neglected effects as the background mean flow, topography, and nonlinearities, to cite the most important ones. Early theories focused on nondispersive long waves, and initially focused on the effects of a background zonal mean flow over a flat-topography as the most likely explanation to account for the discrepancy, e.g., Killworth et al. (1997). The effect of topography was initially discarded by Killworth and Blundell (1999) as being relevant, but Tailleux and McWilliams (2001) demonstrated that the bottom boundary condition, if properly considered, can have a much greater impact than the background mean flow on Rossby wave phase speeds, see also Tailleux (2003Tailleux ( , 2006 for a further discussion of the topographic effects. Over the past few years, the increase in spatial resolution of satellite altimeter products achieved by merging different altimeter datasets motivated the study of the empirical dispersion relation of Rossby waves. This in turn prompted the most recent theories to include dispersive effects as well. The most complete theory to date is probably that of Killworth and Blundell (2005) which considers the effects of dispersion and background mean flow in the presence of a variable large-scale topography, in the context of WKB theory.
While the extended theories appear to do better than the SLT at reducing the discrepancy with observations, e.g., Maharaj et al. (2007), this does not necessarily mean that they do so for the correct reasons. This is because the extended theories, despite their apparent success, still rely on a number of assumptions whose validity can neither be directly checked against observations, nor rigorously proven mathematically. Indeed, such theories are for the most part based on the WKB approximation whose validity relies in principle on a scale separation between the waves and the medium of propagation, as well as on the non-interaction between the different waves supported by the system. While the idea of a scale separation between the waves and the background large-scale circulation does not appear unreasonable at leading order, it seems much more questionable with respect to addressing the effects of topography since the latter exhibits important variations on both small and large scales. The latter issue is not easily dismissed, because the studies by Killworth et al. (1997) and Tailleux and McWilliams (2001) show that the bottom boundary condition have a much greater impact than the background mean flow on the phase speed when both effects are considered separately. With regard to the noninteraction assumption underlying WKB theory, Tailleux and McWilliams (2002) and Tailleux (2004) suggest that it breaks down in places where the topography has a strong curvature, which is likely to occur widely in the actual ocean.
Another important underlying assumption of most WKBbased theories is that planetary waves are able to retain a normal mode structure similar to that of the standard modes, de-spite the linearised equations in the presence of background mean flow and topography being no longer separable. As a result, these extended theories assume that the propagation remains vertically coherent, as for standard normal modes, with the implication that westward propagation should be observable throughout the vertical with a phase speed independent of depth. Of course, such a particular prediction is difficult to test against observations, although this might become feasible in the future using ARGO floats, as initiated by Chu et al. (2007). From a physical viewpoint, the validity of the normal mode assumption is not obvious. Indeed, vertical normal modes are usually interpreted as a standing wave resulting from the superposition of two modes propagating vertically in opposite directions, with the same amplitude but opposed phase speeds. As a result, the persistence of normal modes requires: 1) that the problem be symmetric with respect to upward and downward propagation; 2) that the amplitude of the reflected waves at the bottom and surface boundaries be the same as that of the incident waves. Obviously, this is technically not true for the linearised equations in presence of a background mean flow and bottom topography, so that the degree to which planetary waves can be regarded to be associated with normal modes should probably be more fully justified than has been the case so far. An alternative is to look at the 3-D propagation of planetary wave packets, as proposed by Yang (2000), but then additional difficulties arise with respect to how to satisfy boundary conditions, while the other difficulties associated with the scale-separation and non-interaction remain.
Finally, the last difficulty in interpreting actual planetary wave propagation stems from the possibility that a large fraction of observed westward propagating signals could actually be associated with nonlinear eddies, rather than with linear waves, as recently pointed out by Chelton et al. (2007). That nonlinearities are likely to be important for the study of westward propagating signals is in fact theoretically motivated by the study of Jones (1979) who was the first to demonstrate that linear Rossby waves are always unstable, an issue reinvestigated by Lacasce and Pedlosky (2004). If linear waves are always unstable, and therefore nonlinear interactions important, the question that naturally arises is how accurate and relevant our linear theories can be?
With so many questions about the validity of extended theories, and the linear or nonlinear nature of actual Rossby waves, it appears urgent and desirable to learn more about the full three-dimensional structure of westward propagation in order to make progress toward checking existing theories and possibly proposing suggestions for improvement. (The terminology Rossby waves is used in this paper as a generic descriptor of westward propagating signals, whether these turn out to be linear waves or nonlinear eddies). So-called "realistic" ocean model simulations can fruitfully complement such idealised (analytical or numerical) Rossby wave studies, to investigate the full 3-D (unobserved) structure of the waves. These numerical simulations take into account high resolution topography and geometry, realistic initial state, forcing and stratification to mimic the observed state and evolution of the real ocean, and can thus be directly compared to observations. However, the unambiguous extraction and identification of individual processes and signals is more difficult in such simulations than in idealised models. Hughes (1995) and Hirschi et al. (2007) made use of "realistic" numerical simulations to investigate Rossby waves over the Southern Ocean surface, and the role of Rossby waves in the variability of the MOC, respectively. The outputs of the ATL6 1/6 • Clipper realistic model (Penduff et al., 2004) are used here to study the vertical structure of phase speeds of the waves. The propagation of Rossby waves has a dominant westward component that appears clearly on longitude/time sections (also known as Hovmöller diagrams). In the literature, the two main signal processing techniques that have been used to extract and study surface Rossby waves from such diagrams are the 2-D Radon Transform (RT) and the 2-D Fourier Transform. The RT (Radon, 1917;Deans, 1983) provides a quantitative estimate of the orientation of lines in Hovmöller plots. Chelton and Schlax (1996) first used the RT to study the characteristics (especially the phase speed) of westward-propagating planetary waves in satellite altimeter data. Since then, several authors have applied the 2-D RT on SST or SLA datasets (Hill et al., 2000;Maharaj et al., 2004Maharaj et al., , 2005Cipollini et al., 2001Cipollini et al., , 2006. Superimposed propagating signals locally extracted by the 2-D RT are generally interpreted as various baroclinic modes of Rossby waves (Maharaj et al., 2004;Cipollini et al., 2006). The 2-D RT is used in the present study to extract from observed and simulated data in the subtropical North Atlantic the signals whose propagation speeds are the closest in magnitude to those of the first baroclinic mode as predicted by the extended theory of Killworth and Blundell (2003). Some limitations of the 2-D RT are pointed out and this paper explains how the method has been adapted. Challenor et al. (2001) applied a 3-D RT on altimeter data in the North Atlantic and showed that the meridional phase speed was negligible in many places; this justifies why only the 2-D RT is used here. Some authors (Subrahmanyam et al., 2000;Osychny and Cornillon, 2004;Maharaj et al., 2007) have used the 2-D Fourier Transform instead of the 2-D RT. This alternative approach highlights the spectral components of the data in time and space and provides a different estimate of superimposed phase speeds. Maharaj et al. (2004); Cipollini et al. (2006) explored the differences between both methods.
The purpose of this paper is to investigate and describe the vertical (unobserved) structure of westward propagation speeds in the CLIPPER ATL6 model. To that end, we first seek to build our confidence in the physical plausibility of the model results by first assessing the ability of the model to simulate the speed and amplitude of westward propagating features as observed by satellite altimetry in the subtropial North Atlantic. The primary method to investigate westward propagation is the 2-D RT applied on simulated isopycnal displacements. The datasets are presented in section 2. Section 3 present how we modified the classical 2-D RT to better select the first baroclinic mode of Rossby waves. Our results are presented in Sect. 4, then summarised and discussed in Sect. 5.

Observations
We make use of the merged Topex/Poseidon-ERS Sea-Level Anomaly (SLA) maps distributed by AVISO. SLA fields as a function of (x, y, t) are available weekly over the world ocean from October 1992 onwards on a 1/3 • MERCATOR grid. Details of the mapping method can be found in Le Traon and Ogor (1998); Le Traon et al. (1998). The reader is also referred to the DUACS handbook (http://www.aviso. oceanobs.com) and to CLS (2004) for further details.

Simulations
Numerical data come from the ATL6-ERS26 1/6 • simulation (Penduff et al., 2004) performed during the French CLIPPER project (Tréguier et al., 1999), a high resolution model study of the Atlantic circulation based on OPA 8.1 numerical code (Madec et al., 1988). The ATL6 model resolves the primitive equations on an isotropic 1/6 • MERCATOR grid spanning the whole Atlantic Ocean (98.5 • W-30 • E, 75 • S-70 • N) with 42 geopotential levels. The model is driven by wind stresses, heat and salt fluxes computed from the ECMWF ERA15 reanalysis and subsequent analyses (Barnier, 1998). The model is limited by four open boundaries located at 70 • N, in the Gulf of Cadiz (8 • W), at the Drake Passage (68 • W) and between Africa and Antarctica (30 • E). The model topography in the North Atlantic is presented in Fig. 1. The temporal resolution of the numerical dataset is 5 days. Simulated data used in this work are the model SLA(x, y, t) and the displacements of 9 isopycnals h (x, y, t) located between 750 m and 3250 m. Positive h values correspond to upward isopycnal displacements.  Figure 2 presents the observed and simulated surface eddy kinetic energy (EKE) in the North Atlantic. The observed EKE is computed over the period October 1992-October 2000 from SLA maps by assuming a geostrophic balance (details can be found in Ducet et al. (2000)). The model EKE is computed for the period 1993-2000 as the velocity vari-

Comparison
below the Ekman layer (100 m) to avoid the ageostrophic Ekman components. The path of the main eddy-active fronts is simulated correctly by the CLIP-PER model in the region of interest. As expected from the relatively modest resolution of the model, the simulated EKE is weaker than in reality. Despite a too zonal North Atlantic Current, the Azores Current is correctly located. This simulation has been compared to various observations and was shown to be realistic in many aspects, including its interannual surface variability (see Penduff et al. (2004) and references therein). This static view of ocean variability is complemented by longitude-time (Hovmöller) diagrams. These diagrams allow a dynamic view along zonal sections of the data showing the westward propagation of signals as oblique structures (Fig. 3). The slope of these structures is inversely proportional to the phase speed of westward-propagating signals. The longitude-time plots of the 9 simulated isopycnal immersion depths (shown along the σ 1 =32.05 and σ 2 =36.95 isopycnals, lying around 1000 and 2000 m respectively) exhibit comparable oblique structures, showing the presence of westward-propagating signals at every depth. The magnitude of isopycnal displacements h are much greater than those of SLA. In Fig. 3, most isopycnal displacement patterns are out of phase with SLA, as would be expected for first baroclinic mode propagating structures.

Data processing
This section presents how westward-propagating features are quantified from Sea Level Anomaly fields (SLA) and simulated isopycnal depths. We modified the classical Radon Transform algorithm to improve the precision of the analysis, extract first baroclinic mode signals, and regularise resulting phase speed fields. The same processing technique is applied independently on observed SLAs, simulated SLAs, and on simulated isopycnals. Figure 4 summarises the whole processing.

Radon transform
Observed (resp. simulated) SLA timeseries are interpolated linearly on a regular 1/3 • ×1/3 • (resp. 1/6 • ×1/6 • ) horizontal grid, and locally detrended in time over 1993-2000. We have verified that the distinct zonal and temporal resolutions of the simulated and observed datasets have no significant impact on the estimated phase speeds. Remaining gaps in observed (resp. simulated) fields are filled in the longitudetime space using a 2-D Gaussian interpolation scheme with 1 • ×7-day (resp. 0.5 • ×5-day) search radii and 2 3 • ×1.4-day (resp. 1 3 • ×1-day) half-maxima fullwidths. Superimposed zonally-propagating features appear as oblique patterns in longitude-time Hovmöller diagrams. The classical 2-D Radon Transform (RT, described in e.g. Hill et al., 2000) projects such diagrams from the Cartesian (x, t) plane onto a line oriented at a given θ, as shown in Fig. 5. When the orientation θ of x is normal to the alignments, the RT thus transforms superimposed westward-propagating signals (distinct oblique alignments of positive or negative anomalies in the x, t plane) into positive or negative peaks along the projected coordinate axis, x . In that case the standard deviation of the RT computed along the x axis (noted S(θ ) hereafter) will exhibit local maxima at θ values that can be readily converted into predominant propagation speeds.   Rossby wave amplitudes are weak and difficult to detect polewards of 50 • N. In the equatorial band, their wavelengths may also become longer than our 20 • analysis window and get aliased by the RT . This study thus focuses on the subtropical North Atlantic from 5 • N to 50 • N between 1993 and 2000.

Link between discrete θ angles and phase speeds
The RT projects local Hovmöller diagrams (with dx and dt denoting the zonal and temporal grid discretisation) along oblique lines whose angles with respect to the vertical span discrete θ values. To each angle θ corresponds the phase speed C p =(dx/dt) tan θ, with dt usually fixed and dx varying with the cosine of latitude φ, so that C p =A cos φ tan θ, with A a constant.
So far, most RT implementations, e.g., Hill et al. (2000), have used a uniform discretisation for θ. However, because of the nonlinear relation between C p and θ, this implies a highly non-uniform distribution for the discretised C p , which may impair an accurate determination of the phase speeds. Indeed, crosses and pluses in Fig. 6 show that in that case, δC p would vary by up to a factor of 50 over the range θ=0−80 • (0-20 cm s −1 ). For this reason, we found it preferable to modify the discretisation of θ to yield a uniform discretisation of C p . In the end, the interval 0-20 cm s −1 was split in constant 0.1 cm s −1 steps. Model and observed surface and isopycnal phase speed fields were determined independently using this modified algorithm to ensure a uniform precision throughout the basin and phase speed range.  (SLA(x, t)) intensity onto a line (x axis) at angle θ . p(x , θ) = t SLA(x, t) dt , with x=x cos θ −t sin(θ ) and t=x sin(θ)+t cos(θ). The θ angle for which the standard deviation of p(x , θ) is maximum is related to the preferential orientation in the Hovmöller diagram and its tangent provides the zonal phase speed of the dominant signal.

Smoothing algorithm -extraction of phase speeds
The standard deviation of these Radon Transforms (noted S(x, y, z, C p ) in the following) can then be computed along x lines (see Fig. 5) at each geographical location to extract the strongest propagating features. Figure 7 shows a local example of S(C p ) from altimeter data; local maxima around 6.5 and 1-1.5 cm s −1 are close to the expected phase speeds of first and second-mode Rossby waves, respectively. However, small-scale structures contaminate S(C p ) there, like in many other locations. This noise affects the automatic detection of significant maxima, and the identification of vertical modes. S fields were smoothed along the C p axis to avoid such ambiguities. A 11-point boxcar filter was chosen after various tests to remove scales shorter than 1 cm s −1 without affecting meaningful peaks (see thick grey line in Fig. 7). This smoothed field will be noted S sm (x, y, z, C p ) in the following; it is defined in the range 0.5-19.5 cm s −1 and is much less contaminated by small-scale noise.
At each (x, y, z), dominant S sm (C p ) maxima are then found using the following criteria: (i) maxima must lie between two increasing and two decreasing values of S sm ; (ii) maxima must lie 10 points (1 cm/s) apart from each other. At each location, these maxima yield values of C m (x, y, z)  (m=1, ..., 4), i.e. the phase speeds of the four dominant propagating signals sorted by decreasing amplitude. The same processing has been applied independently to simulated S fields at the surface (z=0) and along each isopycnal. Figure 8 shows zonal and meridional sections of observed and simulated S sm (x, y, z=0, C p ) fields (colour), the corresponding dominant phase speed C (grey lines), and the C m=1−4 fields computed with unsmoothed S (white circles, white dots for m=1). The westward and southward acceleration of Rossby waves is realistically simulated by the model. Intense westward-propagating structures are found in the western basin (a-f) and between 35 • N and 45 • N (a, b, g, h). The 34-36 • N "waveguide" follows the trajectory of the Azores Current, which might contribute to the generation or amplification of the waves (e.g. Cipollini et al., 1999;Cromwell, 2006). As expected, the S sm maximum (C, grey lines) appears less contaminated by small-scale noise than the first S maximum (C m=1 , white dots), not only in the C p direction but also in the horizontal. This beneficial impact is most clearly seen along the western boundary, and between 35 • N and 45 • N.

Selection of phase speeds close to first baroclinic mode predictions
Several authors have interpreted multiple peaks detected by the 2-D RT as the surface signature of baroclinic Rossby modes (e.g. Subrahmanyam et al., 2000;Maharaj et al., 2004), the strongest signal corresponding to the first mode. Observed and simulated phase speeds are now compared at the surface against extended theory estimates for the first two baroclinic modes (Fig. 8). This theory (Killworth and Blundell, 2003) takes into account the effects of a slowly varying baroclinic mean flow and realistic bathymetry. These theoretical phase speed fields are noisier than the C field; they were thus smoothed horizontally using a 2-D boxcar filter (cuboid of 5 • in longitude, 3 • in latitude). Resulting fields will be noted C t1 and C t2 hereafter, for the first two baroclinic modes respectively. As shown in Fig. 8, the observed and simulated dominant surface phase speeds C (grey lines) correspond very well with C t1 (white lines) over the domain of interest. Over most of the basin, the distance between observed and simulated phase speeds is also smaller than the distance between C t1 and C t2 . This can be seen in Fig. 8 along four transects (33 • N, 25 • N, 21 • N, 42 • W) by comparing C(altimeter)−C(model) (grey lines) to C t1 −C t2 (white lines). Small discrepancies between C(altimeter) and C(model) are thus not linked with any contamination from higher order baroclinic modes.
The Radon analysis also detects observed and simulated phase speeds C 2,3,4 (x, y, z) (white circles) in the range of the second baroclinic mode C t2 (white lines). The agreement with theory of these noisy fields is not as good as for C t1 , as already mentioned in Maharaj et al. (2004). Like Maharaj et al. (2004), we conclude that both the observed and simulated dominant peaks (C) are typical of first baroclinic mode Rossby waves in most regions.
Despite the good agreement between C and C t1 , there remains a few places where the first baroclinic mode does not dominate in S sm (grey lines in Fig. 8). We made use of C t1 and C t2 to remove C values outside the first baroclinic mode range. Cipollini et al. (2006) already used theoretical speeds to extract the first-mode Rossby wave signals: the strongest peak was saved only if its speed was within 1/3 and 3 times the predicted theoretical first baroclinic mode speed (Killworth and Blundell, 2003), otherwise the second largest peak was selected provided it was within that range. Our sorting algorithm uses theoretical speeds too, but the criteria takes into account C t2 to explicitly remove the second baroclinic www.ocean-sci.net/4/99/2008/ Ocean Sci., 4, 99-113, 2008  N (a, b), and x=42 • W (g, h). White open circles represent the strongest four maxima of the unsmoothed field S (absolute maxima shown as white dots). Grey lines represent C, the absolute maxima of S sm . White lines represent C t1 and C t2 (see section 3.4 for definitions), C t1 being always greater than C t2 . mode from all the estimations. The C field is sorted in the following way: at each location (x, y, z), C is conserved only if |C−C t1 |< 2 3 (C t1 −C t2 ), otherwise C is masked. Consequently, C is masked where C t2 is not defined (north of the Gulf Stream region). The approach presented above ensures that the phase speeds discussed in the following are close to first baroclinic mode theoretical estimates, and not contaminated by higher-order modes.

Horizontal smoothing of the phase speed maps
Finally, observed and simulated C(x, y, z) fields were smoothed on the horizontal independently (along the 10 levels for the latter) using the 2-D boxcar filter used for theoretical speeds (5 • ×3 • zonally-elongated cuboids). Smoothed maps are masked where less than 30% unsmoothed C(x, y, z) values are available in a cuboid. The amplitude of the strongest westward-propagating features was computed on the same maps (see Appendix A).
Simulated phase speed fields C(x, y, z) are compared with observations at the surface in Sect. 4.1; their simulated vertical structure is presented in Sect. 4.2.

Model/data intercomparison at the surface
Zones of strong wave amplitudes (Fig. 9c-d) correspond to zones of high EKE (Gulf Stream, Azores current). The simulated Gulf Stream is less energetic than observed and slightly shifted to the north. The similarity of the wave amplitude and EKE fields has already been mentioned by Herrmann and Krauss (1989); Osychny and Cornillon (2004) in the Gulf Stream region and by Cipollini et al. (1999);Cromwell (2001) Fig. 1 and time-averaged depth of selected isopycnals (thin black lines). Bottom topography is shown as thick black lines (m), simulated EKE at 100 m depth is shown as red lines (cm 2 s −2 ).
As expected from the theory and as shown by several authors (e.g. Tokmakian and Challenor, 1993;Polito and Cornillon, 1997;Osychny and Cornillon, 2004), zonal phase speeds increase equatorward and, to a lesser degree, westward ( Fig. 9a-b). This feature is well represented by the model, phase speeds reaching 13 cm/s in the southwest part of the domain. Model-observation mismatches (Fig. 9e) reach their maximum (i.e. 2 to 3 cm s −1 ) in the Gulf Stream region, around 10 • N and close to the boundaries of our domain, but become negligible in the Azores Current. Sim-ulated surface phase speeds exhibit an slight slow bias, i.e. −0.4 cm s −1 on average over the domain.
As mentioned above, simulated first-mode phase speeds are globally closer to observational than theoretical estimates. Both observed (Fig. 9f) and simulated (not shown) phase speeds tend to be slower than predicted by the extended theory, except west of the Mid-Atlantic Ridge (MAR) and along 15 • N-20 • N where both are faster. This suggests that the impact of bathymetry might be underestimated by the extended theory in this area. 4.2 Vertical structure of first-mode phase speeds The processing described above has been applied along each of the 9 selected isopycnals independently, thus providing C(x, y, z), i.e. a 3-D estimate of the zonal phase speeds of first mode baroclinic Rossby waves. This approach might help complement both (surface restricted) observational studies and theoretical investigations of 3-D Rossby waves (Tailleux, 2004(Tailleux, , 2006. Panels g and h in Fig. 9 exhibit C(x, y, z) along the σ 1 =32.05 and σ 2 =36.95 isopycnals (lying around 1000 and 2000 m respectively); these two panels correspond to Fig. 3c-d. These maps show that the main features of surface propagation (i.e. equatorward and westward increase of phase speeds) are found at depth as well. More interestingly, panels b, g, and h in Fig. 9 suggest that phase speeds tend to vary with depth; this feature is investigated in more detail in the following.
The ratio C(z)/C(z=0) will be called the Depth-tosurface speed Ratio, noted DR(x, y, z) hereafter. Figures 10,  11, 12 describe the spatial structure of DR(x, y, z) in complementary ways: along selected sections (Fig. 10); averaged vertically between 1000 and 3000 m (Fig. 11); averaged along isopycnals in two subregions and plotted against the median isopycnal immersion (Fig. 12). Typical values of DR range between 0.6 and 1.1, indicating that zonal phase speeds change with depth. The isopycnal medians of DR throughout the whole domain (Fig. 11) reveal an overall and progressive decrease of phase speeds with depth reaching the median value of 15-18% around 3250 m (Fig. 12). This basin-scale decrease is very similar to that seen south of the Azores Current (thick black line in Fig. 12). This deceleration occurs around 1000 m, i.e. below the thermocline (right panel in Fig. 12). However, distinct regimes seem to prevail in various regions (Fig. 11).
The downward decrease of phase speeds is similar in the whole basin and in the 24-30 • N band. In this latitude band, however, isopycnal DR distributions are much narrower, suggesting some kind of regional homogeneity. Zonal phase speeds clearly decrease downward, reaching around 3000 m 75-85% of surface values ( Fig. 10 (1,2,5), continuous black zone in Fig. 11 and 12). Superimposed on this deceleration pattern, a local influence of the MAR on DR(z) appears along 25 • N ( Fig. 10 (5)). DR reaches a local maximum (0.9) along a line which originates east of the top of the MAR and connects to the surface 1500 km further west. This structure has not been studied in detail, but could illustrate large-scale topographic influences.
The situation is different between 31 • N and 37 • N, i.e. around the Azores Current and the southern part of the North Atlantic Current (see Fig. 10 (1-4), dashed black zone in Fig. 11 and 12). The widths of DR distributions below the thermocline show that the regimes found within and south of the Azores Front are clearly distinct. Zonal phase speeds vary much less with depth (by about 5%) in this frontal region where mesoscale eddies are more energetic, and where dynamical nonlinearities are enhanced. Baroclinic energy is transferred into the barotropic mode at eddy scale, possibly enhancing the vertical coherence of oceanic transients in this eddy-active region (where westward-propagating signals are also stronger). One may thus speculate that westwardpropagating features in eddy-active areas might be more   Fig. 11. Isopycnal medians and quartiles of DR are shown respectively as thick lines and horizontal segments for the 24-30 • N (plain line) and the 31-37 • N (dashed line) latitude bands. Depths correspond to the median immersion of each isopycnal in each region. Thin grey lines labeled "1" show the median (vertical line) and quartiles (horizontal segment) of C t1 /C(z=0). Thin grey lines labeled"2": median and quartiles of C t2 /C(z=0). Right panel: Model mean (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) stratification N (z) (s −1 ) averaged between 24 • N and 30 • N (plain line) and between 31 • N and 37 • N (dashed line).
typical of vertically-coherent eddies than of linear waves, thus moderating the downward deceleration of westward propagation. Further investigations are needed to test this hypothesis. Figure 12 shows a very significant difference between the range of simulated phase speeds C(z) and the range of theoretical second baroclinic mode phase speeds C t2 (C(z) being 2.5 to 3 times greater than C t2 ). The substantial decrease of phase speeds below the surface does not bring simulated phase speeds close to second-mode estimates, but keeps them in the range of first-mode expectations.
North and south of the regions depicted above, results are noisier and more disparate (Fig. 10 (6), Fig. 11). At 20 • N-21 • N, the deceleration is weak in the eastern basin and strong in the western basin, with a transition zone above the MAR, where other processes seem to prevail.
It should be mentioned that this downward deceleration of westward phase speeds and its regional patterns are relatively robust and persist in the following cases: (i) no C psmoothing of the S fields; (ii) no horizontal smoothing of the C(x, y, z) field; (iii) use of median filters instead of boxcar filters.

Discussion and conclusion
The two principal objectives of this work were to [1] assess the ability of a realistic high-resolution OGCM to simulate the surface distribution of Rossby wave phase speeds, and to [2] investigate the (as yet unobserved) vertical variations of these phase speeds. The focus was on westward propagating signals whose phase speeds, estimated by means of the Radon Transform, appeared to be close to that anticipated for the first baroclinic mode. The Radon Transform was applied independently on observed and simulated SLA signals, as well as on 9 simulated isopycnal displacements located between 750 m and 3250 m. Several modifications were made to the classical RT algorithm to deal with the presence of noise, the existence of multiple peaks (possibly associated with higher-order baroclinic modes) and thus better extract first-mode signals: (i) regularisation of the S field in the phase speed direction; (ii) extraction of the first baroclinic mode based on theoretical estimates of first-and second-mode phase speeds from the extended theory by Killworth and Blundell (2003); (iii) horizontal smoothing of the phase speeds.
The model/data intercomparison reveals a good agreement between the simulated and observed phase speeds of westward propagating SSH signals, despite an overall slow bias of about −0.4 cm s −1 . Despite this bias, simulated phase speeds appear to match observed phase speeds better than their theoretical counterparts; moreover, both simulated and observed phase speeds tend to be slightly underestimated by the extended theory of Killworth and Blundell (2003), except west of the MAR and along 15 • N-20 • N, suggesting that the effect of bathymetry might be underestimated in the extended theory.
Interestingly, the main and most unexpected finding of this study is the suggestion that the phase speeds of simulated westward propagating first baroclinic mode do not remain uniform vertically, as anticipated from most existing theories, but exhibit a progressive downward deceleration. Typically, phase speeds appear to be slower by about 15-18% (the median value over the domain of interest) around 3250 m compared to their surface value, with the region 24 • N-30 • N being typical of this behavior. Our results also indicate, however, that the degree of downward deceleration is not independent of the geographical location, and that it could be influenced by the degree of nonlinearity of the westward propagating signals. For instance, the downward slow-down of the waves appears to be only about 5-8% between 31 • N and 37 • N where westward propagating signals are perhaps more representative of nonlinear vertically-coherent eddies generated along the Azores current than of linear Rossby waves.
Ocean Sci., 4, 99-113, 2008 www.ocean-sci.net/4/99/2008/ If further confirmed, the present results would likely constitute an important advance in our understanding of westward propagation in the oceans. Indeed, in the linear context, the existence of a vertical dependence of the phase speed of such signals questions the often made normal-mode assumption of many WKB-based theories that the phase speed should be independent of depth. So far, the only study not relying on the normal-mode solutions in a continuously stratified ocean is that of Yang (2000), which focused on the three-dimensional propagation of a wave packet in a background zonal mean flow. Interestingly, the latter study predicts a downward slow-down of the horizontal phase speed in the case of a background zonal mean flow exponentially decaying with depth. Although this is in qualitative agreement with the present results, a more detailed comparison with Yang (2000)'s theory is left for a further study because of its very idealised character. Alternatively, the vertical dependence of the phase speed found in this paper could also potentially imply the inadequacy of linear theories to describe the vertical structure of westward propagation; it could perhaps also result from interference effects resulting from the superposition of various interacting normal modes, as would be expected from a strongly nonlinear signal.
Although the general tendency for downward deceleration appears to be robust, in the sense that none of the modifications made to our analysis algorithm altered the fact that it remained clear throughout the basin, it would probably be useful to carry out further work to better understand the properties of the Radon Transform when applied to situations it was not originally designed for. Indeed, while the Radon Transform is expected to work well for signals whose phase speeds and amplitudes remain relatively uniform over the domain considered, the meaning of its predictions when this is not the case is less clear. Although the issue of non-uniform propagation could be important, there is little evidence so far to support the idea that it would significantly alter the present results.
Long observational time series are not available yet to sample precisely the subsurface signature of Rossby waves. Comparable investigations could be performed in various basins from ensembles of global ocean simulations at various resolutions, such as the one built by the Drakkar Group (2007). This would help better understand the origin of the features revealed in the present study, in particular concerning the effect of different factors (turbulence, topography, stratification, mean flow, etc.) on the vertical structure of Rossby waves.