## Abstract

Three strategies were compared for extrapolating surface geostrophic velocities to below the surface: S1, using only the barotropic or first baroclinic mode; S2, using a fixed or “phase locked” linear combination of the first baroclinic mode and the barotropic mode; and S3, a strategy similar to S2 but using a new set of basis functions. For S2 and S3, the phase locking allows one to impose zero velocity at the seafloor. The new basis functions start from zero at the surface, are not degenerate with respect to the free-surface boundary condition, and represent the adjustment of the pressure at a given depth from density surfaces responding to sea surface height undulations. In idealized primitive equation simulations, strategy S3 had the least error and allowed extrapolation to deeper levels, suggesting the new basis functions performed significantly better than the traditional baroclinic modes. In contrast, strategies S1 and S2 made poor predictions by 400-m depth. Large temporal fluctuations in the fraction of energy in the barotropic and first baroclinic modes could explain the poor predictions by strategies S1 and S2. This brings into question the interpretation of the sea surface height gradients measured by satellite altimetry in terms of first baroclinic mode motions.

## 1. Introduction

With over 17 yr of continuous satellite altimetry, we now have much longer time series of near-surface geostrophic velocities than the typical 1- or 2-yr time series from subsurface moored current meters. What is the best method to extrapolate the surface geostrophic velocity below the surface? To what depth can we obtain a reasonably accurate estimate? Given the importance of satellite altimetry in developing our knowledge of the dynamics of the ocean (Scott et al. 2010a, and references within), these are critical questions. The answers depend upon the dynamical interpretation and vertical structure of the geostrophic ocean currents revealed by the sea surface height (SSH). There are currently two theories of what sets the vertical structure. The first theory suggests that the SSH gradient reflects mostly geostrophic motions of the first baroclinic mode. In a timely and strenuous effort using a global archive of 107 moorings, Wunsch (1997) suggested that the daily averaged currents were dominated by the barotropic and first baroclinic mode. Despite concerns about the grave inadequacy of 107 moorings of limited duration to describe the inhomogeneous World Ocean, this result has remained a corner stone of interpreting altimeter data. For given the surface intensified density stratification typical almost everywhere in the World Ocean, the upshot is that the SSH gradient revealed by altimeter data should predominately reflect first baroclinic mode signals, with the corollary that the abyss is dominated by the barotropic mode. Indeed, this has been found to be consistent with results from quasigeostrophic simulations (e.g., Smith and Vallis 2001; Scott and Arbic 2007).

Recently, Lapeyre and Klein (2006) proposed an alternative dynamical interpretation of altimeter data based upon surface quasigeostrophic theory (SQG) (Held et al. 1995). Lapeyre and Klein (2006) posed the problem in terms of inverting the potential vorticity (PV) (e.g., Hoskins and McIntyre 1985) to obtain the three-dimensional (3D), geostrophically balanced circulation in the upper ocean. Although the interior PV remains unobservable, the surface temperature anomalies are observable and provide information on the inhomogeneous surface boundary. Matching the surface buoyancy boundary condition (BC) to a reference solution (“surface solution” in their terminology) provides an approximation to the upper ocean PV and therefore the complete 3D circulation. The resulting approximate circulation is associated with the so-called buoyancy modes of SQG flow. Lapeyre (2009) found that the SSH gradient signal in an eddying simulation of the North Atlantic Parallel Ocean Program (POP) model could be interpreted as mostly associated with buoyancy modes. LaCasce and Mahadevan (2006) applied the SQG theory to SST data in several regions.

The mathematical motivation for including the buoyancy modes is that the surface buoyancy anomaly boundary condition can be met. However, the price to pay is that the modes are no longer orthogonal in the depth integral. This complicates the representation of the vertical structure because the depth-integrated energy includes an interaction term between the traditional baroclinic and surface buoyancy modes, which can be large and negative (Lapeyre 2009, Fig. 5, where the decomposition of the surface spectral kinetic energy density leads to “components” in the barotropic and SQG mode that are over 100 times larger than the total). If this is deemed a tolerable complication, then one might ask what other formulations that meet inhomogeneous boundary conditions can be considered.

Here, we put forward a complete set of basis functions that meet a specified sea surface height boundary condition rather than sea surface buoyancy condition. This seems a more natural choice because the height boundary condition is directly related to the signal provided by the altimeter data. We call these “Dirichlet basis functions” to distinguish them from the traditional baroclinic modes for reasons that will become clear. We investigate whether these basis functions provide a better representation of the vertical structure of mesoscale (14-day period and longer) oceanic currents.

In representing the vertical structure of ocean currents, two different goals can be sought. Most ambitiously, given the sea surface height anomaly (SSHA), one can hope to predict time series of current anomalies at depth *u _{p}*(

*z*,

*t*) via vertical extrapolation in real time of the surface geostrophic current anomalies to the oceanic interior. For example,

where *u _{g}* = −(

*g*/

*f*)(∂/∂

*y*SSHA) is the surface geostrophic current anomaly and

*F*(

*z*) is the much desired function required for the extrapolation to depth

*z*. In principle

*F*(

*z*) could have a complicated geographical variation, but in practice the function would only be of use if it had some degree of universality: for example, a single baroclinic mode or a fixed (phase locked in time) linear combination of modes. Less ambitiously, one might only require statistics of the currents at depth. The latter goal is fundamentally easier because

*F*(

*z*) can be represented with multiple modes with possibly nonunity correlation (i.e., phase locking is not required). One could then describe, for instance, the vertical structure of eddy kinetic energy from maps of altimeter data, without attempting to predict the currents at a particular place, time, and depth.

Herein we consider the more ambitious type of prediction problem by assessing the error statistics of time series of currents extrapolated from the surface. We compare three different strategies for building the extrapolation function, *F*(*z*): S1, the traditional barotropic or first baroclinic mode only; S2, a phase-locked linear combination of the barotropic and traditional first baroclinic mode; and S3, a phase-locked linear combination of the barotropic, the first, and possibly the second Dirichlet eigenfunction. In particular, we address the following questions: 1) How well can the first (traditional) baroclinic mode predict the observed currents in the upper ocean? The simplistic extrapolation with the barotropic mode provides a rudimentary benchmark. 2) Does the new basis set provided by the Dirichlet eigenfunctions work better? 3) At what depth do we lose predictive skill? In section 2, we briefly review the theory behind the traditional baroclinic modes and then provide the mathematical and physical foundation of the new set of basis functions. We assess the skill of all three strategies in extrapolating the surface geostrophic currents to depth using primitive equation model simulations in section 3. In section 4, we consider a different but related problem: given complete knowledge of the currents at all depths, how many modes are required to describe the vertical structure of the currents? Addressing this related question sheds light on the problem of extrapolating the currents below the surface given only knowledge of the surface currents. We conclude in section 5 by summarizing the answers to the questions posed above, and we recap how the main results from other sections help explain these answers.

## 2. Basis functions for extrapolating surface geostrophic currents below the surface

De Mey and Robinson (1987) considered the extrapolation of altimeter data to below the surface using empirical orthogonal functions (EOFs). A limitation of using EOFs is that one must assume that historical knowledge of the local current statistics is sufficient to construct the EOFs. Although EOFs could be constructed where moorings have been deployed, it is far from obvious how to extend these EOFs to locations where moorings have not been deployed. Watts et al. (2001) have developed an empirical method based upon statistical correlation between the SSH and the observed interior density structure that appears to work well in the Southern Ocean. In principle, such a method might ultimately be extended globally (exploiting Argo observations), but first considerable work needs to be done to build the correlation functions for each region. Given the near-global coverage of altimeter data, we instead seek a method that can be immediately employed almost everywhere. Here, we explore methods based upon basis functions that can be constructed from the climatological density data such as that contained in the *World Ocean Atlas 2005 (WOA2005*; Locarnini et al. 2006; Antonov et al. 2006). Our first two strategies, S1 and S2, used the traditional baroclinic modes. In section 2a, we present some of the background on these modes and their relation to the forced and freely propagating modes of the linear system. In section 2b, we describe an alternative set of basis functions.

### a. Traditional baroclinic modes

The traditional baroclinic modes arise from separation of variables solutions to the governing equations of oceanic motions. Gill (1982, chapter 6) derives the traditional baroclinic modes as the long wave limit of the internal wave modes of a nonrotating, uniformly stratified, flat-bottom ocean. The baroclinic modes are actually more general and can be derived with appropriate assumptions from the primitive equations on a rotating sphere, as presented pedagogically by Wunsch and Stammer (1997). Linearizing the primitive equations, solutions will then separate between time-dependent horizontal functions and time-independent vertical functions as long as the buoyancy frequency and surface and bottom boundary conditions imply consistent horizontal structure. Wunsch and Stammer achieve this with a simple surface forcing term {surface pressure *P _{a}* =

*A*

_{0}(

*k*,

*l*,

*σ*) exp[

*i*(

*kx*+

*ly*−

*σt*]} that separates, in space and time, buoyancy frequency that depends only upon the vertical coordinate and a flat bottom. One then obtains a horizontal partial differential equation (PDE) and separate vertical ordinary differential equations that are coupled through the separation constant

*γ*that depends upon the horizontal scales and the frequency of the surface forcing. For horizontal scales (

*k*,

*l*) falling within the circle defining the dispersion relationship for the nondivergent Rossby waves for a given frequency

*σ*≪ |

*f*| (retaining their notation to simplify consultation),

one obtains vertical structure functions that oscillate in the vertical. If the vertical structure function of the horizontal currents is *F*(*z*), then

with bottom boundary condition

The upper boundary condition is obtained by matching the prescribed pressure and reduces to the homogeneous condition for a free surface when surface pressure forcing is zero, *A*_{0} = 0,

[see also Gill 1982, Eqs. (6.11.15) and (6.11.7)]. With the free-surface boundary condition, Eq. (4), and a homogeneous bottom boundary condition, one obtains the free modes (or resonantly forced modes) of the linear system, which for a flat bottom are the barotropic and baroclinic Rossby waves. The corresponding countably infinite vertical eigenfunctions *F _{m}*(

*z*) are the traditional baroclinic modes, with

*m*= 0 for the barotropic mode and

*m*= 1 for the first baroclinic mode. Imposing a rigid lid the upper boundary condition becomes a Neumann boundary condition [see Eq. (7)],

which in practice gives almost identical modes to the more general homogeneous condition of Eq. (4), as noticed, for example, by Gill (1982), and a point that we will recall below. We will call these modes the traditional or “Neumann” baroclinic modes and denote them with to distinguish them from a set of basis functions derived in section 2b. They are the eigenfunctions of a regular Sturm–Liouville problem and thus form a complete orthogonal basis (Haberman 1998). For *γ*^{2} < 0 and for small scales, the solutions of Eq. (2) approach the buoyancy modes described by Lapeyre and Klein (2006) (see appendix A).

Anticipating our results below, note that, even if the motions of the ocean were completely described by free modes, the motion at a given location can simultaneously contain contributions from different modes. Also, we will see in section 4 that the relative contributions from the barotropic and first baroclinic mode were found to vary widely and rapidly in time in the Modular Ocean Model version 4 (MOM4) simulation.

### b. Basis functions resulting from a homogenous surface Dirichlet boundary condition

In this section, we derive another complete, nonorthogonal set of basis functions. An advantage of this basis set is that the physical interpretation of the basis functions is immediately clear. The barotropic motions are due completely to the free-surface undulations and are immediately associated with the altimeter signal. The remaining eigenfunctions result from vertical displacements of the density surfaces. We refrain from calling this basis set “baroclinic modes,” even though they are eigenfunctions of the same operator as the traditional baroclinic modes, because they do not generally represent free modes of the system, as discussed in the last paragraph of this subsection.

Our starting point and development follows almost identically that of Lapeyre and Klein (2006), except for the imposed surface boundary condition. Lapeyre and Klein posed the problem in terms of inverting the quasigeostrophic potential vorticity anomaly *Q* for the geostrophic streamfunction anomaly *ψ*,

where *f*_{0} is the Coriolis parameter; *N* is the buoyancy frequency; and *z* is depth measured positive upward, with origin (*z* = 0) at the mean sea surface. Before we introduce the new basis functions resulting from the Dirichlet boundary conditions, we first recall the traditional boundary condition solution. With the traditional homogenous Neumann boundary conditions imposed at the sea surface and seafloor, one can solve this standard linear PDE using the method of expanding *Q* in terms of eigenfunctions of the corresponding homogeneous problem *Q* = 0, which separates into horizontal and vertical parts, as described in elementary textbooks on PDEs (e.g., Haberman 1998). The eigenfunctions are the standard baroclinic modes introduced in section 2a and satisfy the regular Sturm–Liouville eigenvalue problem for the geostrophic streamfunction anomaly,

where the eigenvalue *λ* is the inverse of the corresponding Rossby deformation radius and serves as the separation constant between the horizontal and vertical problems. The boundary conditions [Eqs. (7) and (8)] assume that density perturbations at the sea surface and seafloor are small [Gill 1982, cf. (6.11.16) and (6.11.7)]. Lapeyre and Klein (2006) emphasize that although the traditional baroclinic modes provide a complete, orthogonal basis set, surface buoyancy anomalies cannot be represented by them. The implication is that they might not provide an efficient basis set because many modes might be required to expand the vertical structure near the surface. Instead, Lapeyre and Klein used the standard technique for PDE boundary value problems of introducing a reference solution (“surface solution” *ψ*_{sur} in their nomenclature, which is also sometimes called a “particular solution”) that meets the boundary condition. That is, they separate the streamfunction *ψ* into a surface streamfunction *ψ*_{sur} that meets the inhomogeneous surface boundary condition and the remainder *ψ*_{int},

The reference solution is arbitrary other than meeting the boundary conditions (Haberman 1998). Obvious choices are to make it zero everywhere, in which case an infinite second derivative is required at the surface, leading to a delta distribution potential vorticity sheet at *z* = 0 modifying *Q* on the RHS of the resulting equation for *ψ*_{int} (we prove this in appendix B). Lapeyre and Klein (2006) instead chose not to modify the equation for *ψ*_{int} by defining a reference (or surface) solution that satisfies a homogeneous PDE. We belabor these elementary methods of PDEs to emphasize the flexibility and meaning of the buoyancy modes.

Here, we depart from Lapeyre and Klein (2006) by choosing to impose a Dirichlet surface boundary condition. In particular, in keeping with the information provided by the satellite altimeter, we specify that the surface flow match the geostrophic flow implied by the altimeter and choose the reference solution to match this condition. This seems more natural than specifying a surface buoyancy condition that may not be directly related to the SSHA. If we retain the traditional homogeneous Neumann boundary condition at the seafloor, we immediately see that a candidate reference solution can be the *z*-independent function

The resulting problem for *ψ*_{int} satisfies an inhomogeneous PDE,

with homogeneous boundary conditions,

This can be solved by the method of eigenfunction expansion, by first solving the related homogeneous problem with *Q* − ∇^{2}*ψ*_{sur} = 0. The horizontal and vertical problems separate, leading to the regular Sturm–Liouville eigenvalue problem,

The PDE Eq. (12) is identical to that for the traditional baroclinic modes [Eq. (6)], but the different boundary conditions [Eq. (11)] ensure different baroclinic modes or, more properly, different eigenfunctions. For example, a nontrivial constant function cannot meet both boundary conditions. However, recall *ψ*_{sur} was chosen to be barotropic, so a barotropic interior eigenfunction is not needed. The barotropic particular solution is nonorthogonal to the interior eigenfunctions (orthogonality being defined by the unweighted integral of the product of the two functions from *z* = −*h* to *z* = 0). The distinguishing characteristic was the homogeneous Dirichlet surface boundary condition, so we will refer to the eigenfunctions as the Dirichlet basis functions to distinguish them from the traditional “Neumann baroclinic modes.” These are contrasted in Fig. 1.

Because the above eigenproblem is a regular Sturm–Liouville problem, the resulting solutions are countable, infinite, and form a complete set of orthogonal basis functions (Haberman 1998). In general, the first derivative is nonzero at the top boundary so that surface buoyancy anomalies can be represented.

The physical interpretation of the Dirichlet basis functions is immediately clear. The surface solution is barotropic and represents the pressure signal associated with the SSH anomalies. In the absence of interior adjustment, the associated horizontal pressure gradient would be felt throughout the water column. However, the interior of course adjusts by tilting density surfaces so as to reduce the abyssal horizontal pressure gradient. The adjustment process is in general quite complicated (Philander 1978; Müller and Frankignoul 1981; Wunsch and Stammer 1997). In practice, the details of the adjustment processes are not required to make progress. Spatial variations in *f* allow for large-scale mass field adjustment via baroclinic Rossby waves in a process analogous to the inverse barometer response to atmospheric pressure loading described in detail by Wunsch and Stammer (1997).

If a Dirichlet baroclinic eigenfunction was excited in the unforced problem, could it propagate on its own? The Dirichlet eigenfunctions, unlike the Neumann eigenfunctions, were not approximately consistent with the free-surface boundary condition. They are not free modes of the system and should be thought of as only useful basis functions. This is discussed more fully in appendix C.

## 3. Extrapolation of surface geostrophic currents to depth in MOM4 simulations

We consider the simple setting of idealized primitive equation simulations where we are not distracted by measurement noise. The Boussinesq, primitive equation, ocean general circulation model MOM4 was configured in an idealized midlatitude basin and run at a variety of resolutions. The configuration and simulations were described in detail by Cole (2010) and are a minor variation of the runs described by Scott et al. (2010b). The bathymetry includes continental slopes, a large seamount, and a ridge. The vertical resolution was held fixed with 40 vertical levels with spacing varying from 15 m near the surface to 160 m at the deepest point of 3500 m. Two horizontal resolutions were used: ⅓° and latitude and longitude. The equation of state was linear with *α* = 0.0015 K^{−1} and with no salinity dependence. The idealized double gyre zonal wind stress was held constant in time. The simulations were spun up for several decades from rest and with initial potential temperature that decreased exponentially with 700-m *e*-folding depth. We analyzed a year of daily averaged time series in a subregion in the western boundary current extension. All our analysis presented below is for a single column of grid points [fixed (*x*, *y*)] in this subregion.

The traditional Neumann and the new Dirichlet basis functions were calculated from the annual average buoyancy frequency of the run by solving Eq. (6) with boundary conditions Eqs. (7) and (8) for the Neumann modes and boundary conditions Eq. (11) for the Dirichlet basis function. The results are presented in Fig. 1. The model geostrophic currents *υ _{g}*(

*z*,

*t*) were identified from the daily averaged horizontal pressure gradient. Because the MOM4 runs had a residual in the momentum budget resulting from the barotropic equations being solved with shorter time step (barotropic split of 30 was used in these simulations; i.e., there were 30 barotropic time steps for every single baroclinic time step), we diagnosed the effective horizontal pressure gradient from the residual of all other terms as done by Cole (2010), who explored this issue in detail. The geostrophic currents agreed almost exactly with the total current. A 14-day low-pass filter was applied but had negligible effect.

From the surface value of the geostrophic velocity, *υ _{g}*(0,

*t*), we made extrapolations to all depths of the general form,

In Table 1, we give short names to the various schemes we employed. The first letter denotes the basis functions used, either N for Neumann or D for Dirichlet. The numbers following N or D refer to the eigenfunctions included. Determining the constants *a*, *b*, and *c* required imposing various criteria. First, we require the current to match the surface current at *z* = 0, i.e., *υ _{p}*(0,

*t*) =

*υ*(0,

_{g}*t*). For strategy S1, schemes N0 and N1 used only the barotropic or first traditional baroclinic mode, respectively, and that was sufficient to define

*a*or

*b*. Because both S1 schemes gave unrealistically large currents in the abyss, we also tried the S2 scheme, N01, constructed from the barotropic and first traditional baroclinic mode phase locked to match the surface flow and give zero flow at the seafloor,

*z*= −

*h*. We emphasize that we do not anticipate for some physical reason that the modes should be perfectly phase locked but rather that, if we wish to construct a vertical profile of the form Eq. (13) to make a deterministic extrapolation of the surface currents to below the surface, we must make the pragmatic assumption of phase-locked modes. Choosing

*b*to minimize the kinetic energy at the bottom leads to the same value of

*b*as simply setting the bottom velocity to zero.

The strategy S3 made extrapolations using the Dirichlet eigenfunction basis set with the same form, Eq. (14), but we had additional criteria to impose and therefore the possibility of including three terms. In all cases, the only parameters we can solve for are the *a*, *b*, and *c* of Eq. (14). The D0 scheme is identical to N0. The D1 scheme was not considered. For the schemes with two or three terms, additional criteria are required. Again, we imposed zero bottom velocity. Recall that the free-surface boundary condition, Eq. (4), was not degenerate for the Dirichlet eigenfunctions, and thus we could impose that D01 and D012 also meet this boundary condition. The linear system for the parameters *a* and *b* in the case of D01 was overdetermined, and we used least squares to estimate the parameters, with the most weight given to matching the surface velocity, intermediate weight given to meeting the free-surface boundary condition, and the smallest weight given to the zero bottom velocity. We chose a default weighting of (10, 5, and 1) for matching the surface velocity, meeting the free-surface boundary condition, and zero bottom velocity. For sensitivity experiments we used (1, 1, and 1); (2, 1, and 1); (2, 2, and 1); or (2, 1, and 2). The depth where the mean squared (MS) error, , first reached half the variance of the model velocity, VAR(*υ _{g}*) (the “signal-to-noise ratio”), gives a measure of the results. We found varying the weighting as described caused the depth where the signal-to-noise ratio first reaches one-half to vary between 653 and 831 m for the D01 scheme. The worst result was for a weighting of (2, 1, and 2) for matching the surface velocity, meeting the free-surface boundary condition and zero bottom velocity (the result shown in Fig. 3). Removing the free-surface boundary condition constraint reduced the performance of D01, though it was still better than N0, N1, and N01 in the upper 700 m.

A ⅓° simulation was analyzed but the results were not presented because the modeled currents were so weak at depth, the extrapolations all overestimated the strength of the currents below about 800 m. The simulation was found to be much more energetic especially at depth, which is almost certainly more realistic than the ⅓° simulation but probably still has a weak bias in the abyss (Hurlburt and Hogan 2000; Penduff et al. 2006; Arbic et al. 2009; Scott et al. 2010b). A snapshot of the profile of geostrophic currents reveals nonnegligible speeds in the abyss (Fig. 2). The corresponding plot of MS errors in the extrapolated currents (Fig. 3) summarizes the success of the various schemes. The worst case for the D01 scheme was plotted in Fig. 3 and reported in Table 1.

An unexpected result was that the barotropic (N0) and first baroclinic mode (N1) and their phase-locked linear combination (N01) all had similar skill down to almost 400 m. Less surprisingly, below 400 m N0 was the worst prediction scheme. The N01 was better than the N1 scheme below 1000 m, but this does not represent real skill because it only reflects that the currents were weak in the abyss and the N01 scheme exploited the extra degree of freedom to bring the velocities to zero at the bottom. The N01 scheme was slightly worse than simply using N1 in the 500–1000-m depth range. Unlike the other schemes, D01 had nonzero error at the surface, but it was a small fraction of the variance. In the depth range of about 250–800 m, the D01 scheme was the best overall. Only the D012 scheme might be considered superior for the upper ocean because its error went to zero at the surface and gave similar skill as D01 down to about 450 m.

For depth ranges where the MS error in Fig. 3 was much larger than the variance, the prediction scheme was, of course, predicting currents that were systematically too large (e.g., the N0 scheme below about 500 or 600 m or the D012 scheme around 1000 m). On the other hand, the schemes that were forced to zero at the seafloor of course underestimated the currents in the abyss. To reveal these systematic biases, we plotted in Fig. 4 the time mean of the errors,

for each prediction scheme. The sgn(*υ _{g}*) function ensures that the quantity is positive when the prediction is systematically too weak or in the wrong direction and negative when the prediction is in the correct direction but too strong. The D01 scheme was the only prediction with significant positive mean error in the upper few hundred meters. Below about 600 m, the N0 scheme mean error was negative, suggesting the currents were generally equivalent barotropic but, not surprisingly, weaker than the surface current. In contrast, the N1 scheme mean error was positive and larger than the mean |

*υ*| current below about 1100 m, because the predicted currents tended to be in the wrong direction below the zero crossing of the first baroclinic mode (about 1100 m; cf. Fig. 1a). Of course, for the remaining schemes the mean error approached the mean |

_{g}*υ*| in the deep ocean because the prediction approached zero at the seafloor.

_{g}Overall, the Dirichlet eigenfunctions appear to be more useful than the traditional baroclinic modes with Neumann boundary conditions at describing the vertical structure of the MOM4 simulated currents. In the next section, we try to understand these results in more detail by using the currents at all depths to do a decomposition at each point in the time series in terms of both the traditional and new basis functions.

## 4. Modal decomposition of currents given complete knowledge of currents at all depths

With the model data analyzed in the previous section, we have the luxury of information on the currents at all depths. This allows one to decompose the daily currents in terms of either the traditional or new basis functions, in a manner that we define precisely below in Eqs. (16) and (17). This decomposition allowed us to explore how many modes were required to represent the currents to a given accuracy. This is a very different problem than the main problem addressed in the paper. We include this section because it may help shed light on understanding why the new basis functions performed better than the traditional baroclinic modes in the extrapolation of the MOM4 simulated currents. The results of this section also reveal the need to be cautious in interpreting the altimeter data in terms of modes, either baroclinic or buoyancy modes.

The decomposition is defined as

where *z _{i}* with

*i*∈ {1, 2, …, 39} are the discrete vertical levels where the basis functions or were defined. Because we obtained the eigenfunctions from derivatives of the density field, these levels were midway between the 40 model grid levels and we had to linearly interpolate the currents to these levels. Thus, we had 39 eigenfunctions (

*m*= 0, 1, … 38), and in this limit (

*M*= 38) the decomposition becomes exact to within machine precision. The modal coefficients at each time point

*t*where

_{n}*n*∈ {1, 2, …, 355} were obtained with simple least squares with

*M*= 38. We also tried ridge regression with trade-off parameter between 0 and 0.2 (Marquardt and Snee 1975) and obtained almost identical results, so only simple least squares results are presented.

For the decomposition problem posed in Eqs. (16) and (17), a more efficient basis set would allow for more accurate representation of the currents at all depths with fewer terms. To give some measure of this efficiency, we have plotted in Fig. 5 the MS discrepancy between the RHS and LHS of Eqs. (16) and (17) as a function of *M*. Note that *M* = 38 was always used in the fitting procedure, but the reconstruction was explored with various numbers of terms included in the expansion. The MS discrepancy at each point in time was divided by the depth-integrated variance at the corresponding time. The mean is over all 355 days of data, and the error bars represent the statistical uncertainty obtained from boot strapping with 2000 subsamples. This plot clearly shows that the traditional modes were a more efficient basis set by this measure. At first sight, this might seem surprising. How can the traditional modes be a more efficient basis set for the decomposition problem and yet do a poorer job at the related problem of extrapolating the surface velocity to below the surface? To understand this, we must consider the coupling between the modes.

Recall that in the extrapolation problem considered in section 3 we had much less information: we were restricted to knowledge of only the surface currents as opposed to currents at all depths. Going beyond a single-mode fit introduced another degree of freedom that also required making additional assumptions including the assumption of zero velocity at the bottom and phase-locked modes. However, it turns out that the traditional baroclinic modes were only weakly phase locked. A corollary was that the relative portion of energy in the barotropic and baroclinic modes was not constrained but actually varied wildly. In Fig. 6, we plotted the fraction of energy in the first four traditional modes. As expected from Fig. 5, most of the energy was in the barotropic and first baroclinic modes, and only in rare spikes did the higher baroclinic modes surpass 25%. (Modes *m* ≥ 4 contributed a much smaller fraction of the variance and were not shown for clarity. Modes *m* = 4 and *m* = 5 had peaks reaching 11% and 13%, respectively, whereas even higher modes had monotonically decreasing maximum peaks.) As a result, the fraction of energy in the barotropic and first baroclinic mode tended to compensate each other: when one was larger the other was smaller. The difficulty for the extrapolation problem was that the fraction of energy in the barotropic and first baroclinic modes varies wildly, approaching 95% one day and then plunging to near 0% a few days later. Although the fraction of energy in these two dominant modes must be strongly anticorrelated (Pearson correlation coefficient was *r* = −0.82), the modal coefficients and were more free to vary independently. The latter point is emphasized in the scatterplot of versus in Fig. 7. The correlation coefficient between and was only *r* = 0.67, albeit statistically significant at the 99% level assuming at least 12 degrees of freedom or more, which was rather weak. In contrast, the corresponding coefficients for the new basis functions, and , were very strongly anticorrelated with *r* = −0.98. Thus, the phase-locking assumption that was exploited in the extrapolation problem was much more applicable for the Dirichlet basis. The extrapolations with the Dirichlet basis also benefitted from imposing the criterion of meeting the free-surface boundary condition.

## 5. Conclusions

In searching for a function to extrapolate the surface geostrophic currents below the surface, we started from the problem of inverting PV, as suggested by Lapeyre and Klein (2006). Instead of choosing a particular solution that matched the sea surface density, we chose one that matched the sea surface height because this can be directly observed with satellite altimetry. This changed the boundary conditions for the Sturm–Liouville eigenfunctions from Neumann to Dirichlet, which has a first order effect on the eigenfunctions (cf. Fig. 1).

Returning to the questions posed in the introduction, we conclude the following: 1) The first (traditional) baroclinic mode had some skill in extrapolating the surface currents in the upper ocean, though not much more than the simple extrapolation with the barotropic mode. Using both the first baroclinic and barotropic modes by assuming complete coupling between them (i.e., phase locked) did not degrade the skill in the upper ocean and improved it below about 1200-m depth. 2) The prediction scheme using a phase-locked linear combination of the first two terms of a new set of basis functions (the “Dirichlet eigenfunctions”), scheme D01, had similar skill to the traditional Neumann baroclinic modes in the upper 200 m and was much better between about 300 and 700 m. The scheme employing the first three eigenfunctions, D012, was similar to D01 in the upper 500 m or so but predicted currents far too strong below that. 3) All schemes predicted currents that remained significantly correlated to the actual currents throughout most of the water column. How deep one can extrapolate the surface currents, however, depends upon one’s tolerance for expected error. Conservatively, one might require a signal-to-noise ratio greater than 2. By this measure the schemes with traditional baroclinic modes, scheme N0, N1, or N01 lost predictive skill by 380-, 412-, or 400-m depth, respectively. The schemes with Dirichlet eigenfunctions allowed for deeper extrapolations. The D01 and D012 schemes had signal-to-noise ratio above 2 in the upper 653 and 535 m, respectively. Tuning the choice of weighting in solving the overdetermined system to obtain the coefficients *a*, *b*, and *c* in Eq. (14) improved these results.

The other parts of this paper were aimed at helping to make sense of the above results. From the background on linear modes given in section 2a, one can anticipate several problems with using the traditional baroclinic modes to extrapolate surface currents to below the surface. First, the assumptions that one had to introduce to obtain separation of variables solutions to the governing equations do not strictly apply to the ocean. Of more concern is the fact that much of the variability in the ocean lies well outside the region where the separation constant *γ*^{2} > 0, so that the so-called negative-depth modes might be more appropriate than the baroclinic modes. We showed that these negative-depth modes reduce to the buoyancy modes of SQG well outside the *γ*^{2} = 0 circle that defines the dispersion relation for nondivergent Rossby waves; for example, in the limit where . [So for motions with 10-day period at midlatitudes and wavelength smaller than 200 km, the negative-depth modes look like the SQG modes. These are on the small end of what can be resolved with present altimetry but are well within the scope of planned future missions (Scott et al. 2010a).] Even within the region where *γ*^{2} < 0, where the baroclinic modes apply, for subinertial motions the discrete baroclinic modes correspond to the freely propagating or resonantly forced Rossby waves. It is not clear why these should dominate geostrophic turbulence or nonresonantly forced linear modes.

The modal decomposition in section 4 performed for the MOM4 simulation provided further insight. The traditional basis functions were a more efficient basis, with the first two modes capturing almost 85% of the depth-integrated variance compared to only about 70% for the new basis functions (cf. Fig. 5). Despite the traditional modes being a more efficient basis set, they performed worse for the extrapolation problem, which can be understood in terms of the degree of modal coupling. Using a linear combination of two eigenfunctions allowed for better extrapolations (scheme N01 and D01) but required the additional assumption of perfect temporal correlation or “phase locking” between the coefficients multiplying the eigenfunctions. In reality, the coefficients are only partially correlated. (Using realistic correlations would allow for extrapolation to depth of statistical properties only, such as moments of the velocity time series.) Indeed, the lower temporal correlation (*r*^{2} = 0.45) between the traditional barotropic and first baroclinic mode coefficients accounted for the lower ability to extrapolate below the surface. In contrast the barotropic and first Dirichlet eigenfunctions were highly temporally correlated (*r*^{2} = 0.96), which allowed more accurate extrapolation below the surface.

In exploring the modal decomposition in the idealized MOM4 simulations we found the fraction of energy in the barotropic and first baroclinic mode was highly variable in time (cf. Fig. 6). This variability has grave implications for the interpretation of the altimeter data in terms of the first baroclinic mode. Although these model results were consistent with the idea that the barotropic and first baroclinic modes account for much of the kinetic energy and this energy is approximately equally distributed between the two modes on average, the ocean model at any given instance departed widely from this average. The currents could change from being over 95% barotropic on one day to less than 5% barotropic a few days later. Wunsch (1997) presented only mean statistics on the modal decomposition. Also of great relevance for interpreting the altimeter data is the variability about the mean. This variability about the mean is presently unknown for the real ocean and apparently highly significant in this model run.

## Acknowledgments

The MOM4 code was obtained from NOAA/GFDL and we especially thank Bill Merryfield for help in configuration and Matt Harrison and Steven Griffies for help running MOM4. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high-performance computing (HPC) resources that have contributed to the research results reported within this paper (http://www.tacc.utexas.edu). RBS and DF were supported by NSF Grants OCE-0526412 and OCE-0851457; a contract with the NOC, Southampton; and a NASA subcontract to Boston University. This work was motivated by conversation with Guillaume Lapeyre in 2006. We benefitted especially from helpful conversations with Anand Gnanedesikan (on the adjustment problem) and Glenn Flierl (regarding Sturm–Liouville theory). The manuscript also benefited from conversations with Patrice Klein, Xavier Carton, and Guillaume Routlet. An anonymous reviewer requested section 4.

### APPENDIX A

#### Negative-Depth Modes and SQG Buoyancy Modes

For *γ*^{2} < 0 and for small scales, the solutions of Eq. (2) are qualitatively different than with *γ*^{2} > 0; for example, with *N*(*z*) a monotonically increasing function, then the vertical structure functions *F*(*z*) are also monotonically increasing, becoming exponential functions when *N* is independent of *z*, so *F* = *F _{s}* exp(|

*γ*|

*Nz*). These are the so-called negative-depth eigenfunctions discussed more generally by Philander (1978). Interestingly, for motions well outside the circle

*γ*

^{2}= 0 (i.e., in the limit where ), we find from Eq. (1) that

and the negative-depth eigenfunction solutions have precisely the vertical structure of the buoyancy modes described by Lapeyre and Klein (2006).

### APPENDIX B

#### The Singular PV Sheet at the Ocean Surface

Bretherton (1966) considers quasigeostrophic flow with inhomogeneous surface boundary conditions. The problem is framed in terms of the meridional PV flux but the important point is the treatment of the PV itself and especially near the boundary. He pointed out the following:

The boundary contribution may be thought of as due to eddy flux in a sheet of potential vorticity along the boundary. Any flow with potential temperature variations over a horizontal rigid plane boundary may be considered equivalent to a flow without such variations, but with a concentration of potential vorticity very close to the boundary. The dynamical consistency of this picture is assured by the persistence with time of Eq. (3), together with Eq. (5) describing conservation of potential vorticity in the interior.

His Eq. (3) expressed the boundary condition,

which can of course be made homogeneous through a trivial variable change *ψ*′ = *ψ* − *ψ*^{0} where satisfies his Eq. (3).

Here, we show that this infinite sheet of PV can be obtained from taking the limit of a particular solution. This complements the two-layer interpretation with the upper layer approaching zero thickness (Boss et al. 1996). Recall in section 2b that we suggested one obvious choice for the reference solution was to make *ψ*_{ref} zero everywhere, in which case an infinite second derivative is required at the surface, leading to a delta distribution potential vorticity sheet at *z* = 0 modifying *Q* on the RHS of the resulting equation for *ψ*_{hom}. The homogeneous boundary conditions are still imposed on *ψ*_{hom}, but *ψ* = *ψ*_{hom} almost everywhere because *ψ*_{ref} → 0.

Is it possible to choose a particular solution with these properties? Can we make *ψ*_{ref} → 0 while still meeting the boundary condition? Is the singularity confined to the boundary? Let us look in more detail at the nature of the particular solution in this limit. Let us start with a specific functional form

where *a* and *z*_{0} are free parameters that we will vary to obtain the solution with the desired properties. First, we require that the particular solution meets an inhomogeneous boundary condition at the surface. Without loss of generality, in suitably chosen units,

We are interested in the limit *a* → +∞, for then we have

This is of course very convenient because then we do not have to explicitly solve for *ψ*_{ref}. However, what happens to the PV in this limit? Is the singularity confined to the boundary and is it integrable? Near the boundary, we anticipate the stretching term dominates,

At the boundary *z* = 0, we have in the limit *a* → +∞,

However, immediately inside the domain, say at *z* = −*ε* where 0 < *ε* ≪ 1, we have

The stretching term of the particular solution remains finite within the domain, even arbitrarily close to the surface because the exponential factor overpowers the factor linear in *a*. Of course the singularity is integrable because ∂*ψ*_{ref}/∂*z* = 0 for *z* < 0. This justifies the representation of the inhomogeneous boundary condition as a singular PV sheet of the form

for then *Q*_{ref} has the desired properties: (i) singularity is confined to *z* = 0, and (ii) integrating from *z* < 0 to *z* = 0 gives the expected result

Reiterating, because *ψ*_{ref} = 0 almost everywhere, we do not have to solve for this particular solution explicitly.

A few subtleties are worth emphasizing. Note that, with the method outlined above, we actually solve the PDE with homogeneous boundary conditions [Eq. (10)] with RHS,

and conveniently obtain the solution to the inhomogeneous BC problem almost everywhere. Only at the boundary where ∂*ψ*_{hom}/∂*z*|_{z}_{=0} = 0 is the homogeneous solution defective. However, this deficiency is of no practical consequence because the value of ∂*ψ*/∂*z*|_{z}_{=0} is of course known. This subtlety has escaped at least one pair of authors, who in a well-known paper (Rhines and Holland 1979) wrote

This is wrong because they have confused the solution for the homogeneous problem with that of the full solution precisely at the point where it matters (*z* = 0). Another subtlety is that the singular PV sheet arises as an artifact of using this trick to solve for the complete solution while imposing homogeneous boundary conditions. If we were to choose a different reference solution, then we do not necessarily have a singularity in the reference solution PV at the boundary. Green (1987) referred to the singular PV sheet as B3, which is short for “Bretherton bogus boundary condition”; see Green (1987) for complementary comments on the potential for confusion arising from the Bretherton (1966) device.

### APPENDIX C

#### Free Modes versus Eigenfunctions

For the Dirichlet eigenvalue problem, we have imposed an unusual surface boundary condition and so it is necessary to check that the solutions also meet the boundary conditions of the free surface [Eq. (4)]. For constant *N*, they will not because the reduce to sine functions (e.g., the first mode is the sine function on the interval [0, *π*/2]), so the first derivative is maximum at the surface where a free mode would have it to be zero. One might note from Fig. 1 that it appears the first derivative of our is small near the surface. This can be understood from the *N*^{2} profile and the ODE that the eigenfunctions obey. Everywhere, the eigenfunctions satisfy the ODE,

At the surface, we have imposed *F* = 0, so for large then the low eigenfunctions will have *dF*/*dz* small at the surface. In our solution, we indeed found small *N* near the surface associated with the well-mixed layer and large *N* below associated with the base of this layer, a generic feature of the ocean. The diagnosed at the surface, although small, is multiplied by a large constant *g*/*N*^{2} in Eq. (4) and then implies , or about 25%–33% of the maximum value, where the maximum was taken over depth for the given mode. The Dirichlet eigenfunctions, unlike the Neumann eigenfunctions, were not approximately consistent with the free-surface boundary condition.

This should be contrasted to the situation when one imposes a rigid lid, which forces *dF*/*dz* = 0 (and produces the Neumann eigenfunctions), which then gives a value of *F*(0) that when substituted into the free-surface condition Eq. (4) implies a nonzero but small derivative. The approximation is much better in that case. For comparison, for our *N* profiles, the rigid lid approximation implied a surface density anomaly that was typically about 4% of the maximum density anomaly, where the maximum was taken over depth for the given mode. In some sense, we have turned to our advantage the fact that the Dirichlet eigenfunctions do not redundantly approximate the free-surface boundary condition. In section 3, we could impose the free-surface boundary condition to further constrain the extrapolation of the surface currents to depth.

## REFERENCES

## Footnotes

Current affiliation: Laboratoire de Physique des Océans UMR6523, CNRS/UBO/IFREMER/IRD, Brest, France.

^{+}

Current affiliation: National Oceanography Centre Southampton, Southampton, United Kingdom.