## Abstract

Using observations from the Argo array of profiling floats, the large-scale circulation of the upper 2000 decibars (db) of the global ocean is computed for the period from December 2004 to November 2010. The geostrophic velocity relative to a reference level of 900 db is estimated from temperature and salinity profiles, and the absolute geostrophic velocity at the reference level is estimated from the trajectory data provided by the floats. Combining the two gives the absolute geostrophic velocity on 29 pressure surfaces spanning the upper 2000 db of the global ocean. These velocities, together with satellite observations of wind stress, are then used to evaluate Sverdrup balance, the simple canonical theory relating meridional geostrophic transport to wind forcing. Observed transports agree well with predictions based on the wind field over large areas, primarily in the tropics and subtropics. Elsewhere, especially at higher latitudes and in boundary regions, Sverdrup balance does not accurately describe meridional geostrophic transports, possibly due to the increased importance of the barotropic flow, nonlinear dynamics, and topographic influence. Thus, while it provides an effective framework for understanding the zero-order wind-driven circulation in much of the global ocean, Sverdrup balance should not be regarded as axiomatic.

## 1. Introduction

The general circulation of the ocean plays a central role in the global transport of heat, freshwater, carbon, oxygen, nutrients, and other constituents. The ocean’s large-scale circulation strongly influences the distribution of these important quantities and thus shapes both the global climate and patterns of marine biological production. A fundamental zero-order theory of the large-scale circulation, Sverdrup balance (Sverdrup 1947) directly relates meridional geostrophic transport to wind stress. Despite its key role in our understanding of ocean dynamics, the extent to which this tenet quantitatively describes the observed global ocean circulation has received surprisingly little investigation.

One of the most likely causes for the lack of analyses of Sverdrup balance is the dearth of subsurface velocity observations. Subsurface measurements have traditionally been limited in both space and time because of the resource-intensive nature of shipboard observations and the sheer magnitude of the world’s oceans. In his seminal series, Reid (1989, 1994, 1997, 2003) compiled over six decades of hydrographic and tracer data to estimate the mean geostrophic circulation of the Atlantic, Pacific, and Indian Oceans. His analysis relied on somewhat subjective choices for constraints on the total transport and reference velocities and also assumed that both the flow and mixing take place along isopycnal surfaces. Several other studies (Macdonald and Wunsch 1996; Macdonald 1998; Ganachaud and Wunsch 2000; Ganachaud 2003) have used hydrographic sections from the World Ocean Circulation Experiment (WOCE) in conjunction with inverse box models to produce quantitative estimates of global ocean circulation. Such estimates, however, are sensitive to particular assumptions concerning the accuracy of the conservation equations used to constrain the flow. Although these analyses represent considerable progress in our ability to determine subsurface geostrophic flow, the spatial and temporal coverage of the data used in these estimates remains relatively narrow.

In contrast, the Argo array of more than 3000 autonomous profiling floats (Roemmich et al. 2004) offers an abundance of hydrographic and velocity observations with unprecedented resolution in both space and time. These instruments provide velocity estimates and measure temperature and salinity in the upper 2000 decibars (db) of the ocean, relaying the data to shore-based data centers in real time. Argo floats have been deployed throughout the global ocean with a nominal horizontal spacing of 300 km. Each float collects hydrographic profiles at regular intervals, typically 10 days. In this study, absolute geostrophic velocities are determined by combining these hydrographic profiles with the direct estimates of velocity provided by the float trajectories. The resulting quantitative estimate of the large-scale circulation has significant advantages, including vastly improved spatial and temporal data coverage and no ad hoc assumption regarding the reference velocity.

The goals of this work are to 1) compute the absolute geostrophic velocity in the upper 2000 db of the global ocean during the period from 1 December 2004 to 30 November 2010 using observations from the Argo array of profiling floats and 2) use the resulting velocity estimate, together with satellite observations of wind stress, to examine the validity of Sverdrup balance on a global point-by-point basis. Given that this highly simplified theory of the wind-driven circulation is not expected to correspond exactly to the observed ocean everywhere, determining if and when it quantitatively predicts geostrophic circulation is valuable. Investigating where Sverdrup balance succeeds and where it fails, and understanding why, will also help further our knowledge of the dynamics that dominate different ocean regimes. Section 2 gives a brief derivation of Sverdrup balance and reviews relevant results from previous work. The data and methods used to compute the gridded velocities and evaluate Sverdrup balance are discussed in sections 3 and 4, respectively. Section 5 presents the absolute geostrophic velocity fields followed by an analysis of the validity of Sverdrup balance and the implications of these findings. Last, section 6 reviews the conclusions of this study.

## 2. Background

Since Sverdrup (1947) first outlined the elegant formulation of the wind-driven circulation now known as Sverdrup balance (hereafter SB), it has become a canonical part of circulation theory. The derivation of SB is based on the linear vorticity equation connecting meridional geostrophic velocity *υ*_{g} to the vertical gradient of vertical velocity *w*:

Following standard notation, *f* is the Coriolis parameter and *β* is its meridional gradient. The validity of (1) requires steady flow as well as small Rossby and Ekman numbers. Integrating vertically from a depth *z* = *h* to the base of the Ekman layer (denoted as *z* = 0) gives

where *V*_{g} is the meridional component of the geostrophic transport. Basic Ekman theory provides that the vertical velocity at the base of the Ekman layer is related to the wind stress vector ** τ** = (

*τ*

^{x},

*τ*

^{y}) through

where *ρ* is the density and **k** is the unit vector in the vertical direction. Next we make the critical assumption that *w*(*z* = *h*) = 0. Given this, (2) and (3) can be combined to produce

Expanding the right-hand side of (4) gives a useful partition in terms of the meridional component of Ekman transport *V*_{E} and what is commonly thought of as the Sverdrup transport *V*_{Sv}:

This statement of SB separates *V*_{g} from the wind-derived transport *V*_{Sv} − *V*_{E}, while also showing that the total meridional transport *V* = *V*_{g} + *V*_{E} equals *V*_{Sv}. The form of SB provided by (5) is directly tested in this study, using observations provided by the Argo array.

Relatively few studies have attempted to assess the validity of SB using observational data, and those that have are mostly limited to fairly narrow regions. Leetmaa et al. (1977) first explored a zonally integrated form of SB in the North Atlantic. Using a hydrographic transect across 25°N and assuming a level of no motion, they found that the zonally integrated southward transport was consistent with SB. Wunsch and Roemmich (1985), however, identified several issues with this study that called this conclusion into question. These included the extrapolation of a level of no motion through the Gulf Stream recirculation, contrary to the evidence there; inconsistency with the meridional heat transport needed to produce a substantial quantity of deep water at high latitudes; and a contradiction between the sign of the bottom *w* obtained by integrating (1) downward from a surface Ekman pumping velocity and the bottom *w* expected based on continuity and the slope of the seafloor. At the same latitude, Schmitz et al. (1992) subsequently concluded that the simple balance between southward interior geostrophic transport and northward Gulf Stream transport assumed in the Leetmaa study was not supported by water mass characteristics.

In the North Pacific, Qiu and Joyce (1992) determined that interannual variations in Sverdrup transport were strongly correlated with the transport of the Kuroshio and the North Equatorial Current. Hautala et al. (1994) found good agreement between the zonally integrated meridional geostrophic transport inferred from a hydrographic section across 24°N and that computed from wind data according to SB. Using altimetric data, Vivier et al. (1999) examined a time-dependent topographic SB, in which the integral in (2) is computed from the seafloor at *z* = *H* to the sea surface and *w*(*z* = *H*) is given by continuity [*w*(*z* = *H*) = **u** · **∇***H*]. They found that this version of SB explained variations in sea surface height in much of the North Pacific and the subtropical South Pacific.

More recently, Wunsch (2011) carried out a global quantitative analysis of SB using the ocean state estimate produced by the Estimating the Circulation & Climate of the Ocean-Global Ocean Data Assimilation Experiment (ECCO–GODAE) project from a least squares fit of general circulation models to global-scale datasets available after 1992. On a point-by-point basis, the form of SB given by (2) was tested by integrating the meridional velocity, averaged over 16 years, from the depth where |*w*| < 10^{−8} m s^{−1} to the surface. This transport was then compared to *f*/*β*[*w*(*z* = −117 m)], where *w*(*z* = −117 m) was assumed to be the Ekman pumping velocity. Using that formulation, SB was found to be reasonably accurate over much of the interior subtropical and tropical oceans. In this study, we provide another global assessment of the validity of SB, based solely on data obtained as part of the ongoing Argo program together with satellite wind stress observations.

## 3. Data

### a. Oceanographic data

The oceanographic data in this study were collected by the Argo array of autonomous profiling floats and consist of upper ocean profiles of temperature *T* and salinity *S* as a function of pressure *p*, along with float trajectories determined by satellite positioning. The Argo array is composed of several different types of quasi-Lagrangian profiling floats. While the individual float types vary in terms of their specific components, all the floats used here were equipped with a buoyancy engine, a satellite antenna for positioning and transmission, and an instrument package containing a conductivity–temperature–depth (CTD) sensor. To ensure high quality, nearly all of the data included in our study (99.9% of the total profiles) were collected by floats with SeaBird Electronics CTDs, with the remaining small quantity of data from a selected subset of floats with Falmouth Scientific, Inc. (FSI) CTD sensors.

Each float repeatedly executed a mission cycle that consisted of four phases: 1) descent to the predefined “parking pressure” *p*_{park}; 2) subsurface drift at *p*_{park} for a predetermined length of time; 3) ascent to the surface with concurrent measurements of *T*, *S*, and *p*; and 4) surface drift with transmission of the hydrographic data and position fixes to a satellite. The data transmitted from the floats were compiled at two onshore Global Data Assembly Centers (GDACs) and were made publicly available usually within 24 h. More details regarding the technical aspects and operation of the floats as well as the processing of the data are given by Roemmich et al. (2004) and the Argo Data Management site (http://www.argodatamgt.org/Documentation).

All data used in this study were obtained from the Coriolis GDAC at the French Research Institute for Exploitation of the Sea (IFREMER), France (www.coriolis.eu.org). Only data that had undergone thorough quality control testing and were designated as delayed mode were included. In addition, data from floats with known problems (i.e., those on the “greylist” compiled by the GDACs) were discarded. Data from marginal seas, such as the Gulf of Mexico and the Indonesian seas, were also excluded. The time period for this analysis spanned from 1 December 2004 to 30 November 2010. Prior to this time period, global Argo data coverage was sparse, and complete global delayed-mode data were not available for later years at the time of preparation.

#### 1) Profiles

A total of 485 704 *T* and *S* profiles from 5027 Argo floats were used in this study. The profiles had variable vertical resolution, with observations typically at 5- or 10-db increments near the surface and gradually increasing to 100-db increments below 1500 db. Only profiles extending to 870 db or deeper were included in this dataset. The data are for the most part evenly distributed in space (Fig. 1a), although the edges of the basins are not as well covered as the interior. The temporal distribution of the data is skewed toward the later part of the study period (Fig. 1b) because the Argo array did not reach full coverage until early 2007.

At each profile location, *T*, *S*, and *p* were used to compute 1) a profile of dynamic height *D* as a function of *p*; 2) profiles of potential temperature *θ*, *S*, and *p* as a function of potential density *σ*_{θ}; and 3) maximum *p*, average *θ*, and average *S* of the mixed layer. To determine *D*, specific volume anomalies *δ*(*S*, *T*, *p*) = 1/*ρ*(*S*, *T*, *p*) − 1/*ρ*(35, 0, *p*) were calculated and then integrated with respect to pressure according to

The reference pressure *p*_{0} was set to 900 db to allow easy combination with the velocity data derived from the float trajectories. Dynamic height was computed at 29 regular pressure levels spanning the upper 2000 db. Extrapolation of the *D* profile beyond the pressure range of the original *T* and *S* profiles was limited to less than 28 db at the bottom of the profile and less than 4 db at the top.

The integral in (6) was calculated by interpolating the *δ* profile with a piecewise polynomial composed of cubic Hermite basis functions (Fritsch and Carlson 1980) and then analytically integrating the interpolant. This form of interpolating polynomial was chosen because it efficiently produces a smooth function that is monotonic on each interval and is therefore shape preserving with no superfluous “bumps” or “wiggles.” The resulting interpolated profiles were visually appealing and had minimal introduced density inversions. To further ensure the quality of the profile data, we manually inspected all interpolated profiles. Profiles that contained gaps of more than 400 db, density inversions ≥0.02 kg m^{−3} db^{−1}, or unrealistic interpolations (as determined by comparison to nearby profiles) were either corrected by removing a portion of the data or discarded entirely if correction was not possible.

The deepest *p* of the mixed layer along with the average mixed layer *θ* and *S* were computed for each profile using the density algorithm of Holte and Talley (2009). Profiles for which the algorithm failed to produce a mixed layer *p*, such as those that did not include any *p* < 20 db, were excluded, resulting in a dataset of 476 319 mixed layer data points. For each profile, *σ*_{θ} was computed, and then *θ*, *S*, and *p* were interpolated onto regular *σ*_{θ} levels, again using a piecewise cubic Hermite polynomial without extrapolation. To span the full range of the Argo profiles used here, the *σ*_{θ} levels ranged from 22.4 to 27.95 kg m^{−3} and were linearly spaced with decreasing intervals from 0.25 to 0.01 kg m^{−3}.

#### 2) Trajectories

Trajectory data for each float were used to estimate the drift velocity at the parking pressure **u**(*p*_{park}) for each cycle. Only floats with a nominal *p*_{park} ≥ 800 db were included in this study. Drift velocities were computed using the time the float began its descent *t*_{des} and the time it reached the surface at the end of the subsequent ascent *t*_{asc}, together with the corresponding diving location **x**_{des} and surfacing location **x**_{asc}:

The *t* and **x** of the velocity estimate were taken as the average of *t*_{asc} and *t*_{des} and the average of **x**_{asc} and **x**_{des}, respectively. This analysis only used subsurface drifts for which 5 ≤ (*t*_{asc} − *t*_{des}) ≤ 25 days. Shorter cycles do not adequately measure geostrophic velocity due to a stronger influence from high-frequency fluctuations, while longer cycles result in excessive smoothing in space and time.

For each cycle along the float trajectory, **x**_{des} and **x**_{asc} were determined with one of two methods depending on the type of satellite positioning system used by the float, the Argos system (96.6% of the data used here) or global positioning system (GPS; 3.4%). For Argos floats, which remain at the surface for 6–12 h while obtaining multiple position fixes of varying accuracy, the procedure followed that detailed by Park et al. (2005). Both **x**_{asc} and **x**_{des} are generally unknown because a variable amount of time elapses between *t*_{asc} and the first position fix and between the last position fix and *t*_{des}. The float experiences much stronger currents at the surface than at *p*_{park}, and thus the accuracy of **u**(*p*_{park}) is greatly improved if **x**_{asc} and **x**_{des} are estimated taking into account the surface drift and the elapsed time (Park et al. 2005). To accomplish this, *t*_{des} and *t*_{asc} are needed. These times were either obtained directly from the GDAC when available or estimated from the float metadata and/or surface fixes using the protocols given in Park et al. (2005). When it was not possible to calculate *t*_{asc} and *t*_{des} from the available data for a particular cycle (37.6% of the data), they were approximated by computing the average Δ*t* between *t*_{asc} and the first surface fix and the average Δ*t* between *t*_{des} and the last surface fix for all other cycles of the float and then adding these to the first and last fixes, respectively. In the small number of cases where we were unable to estimate *t*_{asc} and/or *t*_{des} by any of these means (3.8%), the times of the first and last surfaces fixes were used instead.

Next **x**_{asc} and **x**_{des} were estimated via a least squares fit of the fixes from each surface drift to a theoretical trajectory **x**(*t*), derived from a linear velocity **U**_{L} and an inertial oscillation of magnitude *U*_{I} and phase *ϕ*:

where *t*_{0} is the time of the first fix. Each surface fix was weighted by 1/*ϵ*^{2}, where *ϵ* is the standard deviation of the fix (150, 350, and 1000 m for Argos fixes of level 3, 2, and 1, respectively). The fitted trajectory was then extrapolated to *t*_{des} and *t*_{asc} to give **x**_{des} and **x**_{asc}. Park et al. (2005) showed that this method produces velocities that are more accurate and less biased than simple linear extrapolation of the surface fixes. The total uncertainty of each velocity estimate computed in this manner is on the order of 0.2 cm s^{−1} (Park et al. 2005).

Floats using GPS positioning obtain a single surface fix with a much higher accuracy (10 m) and remain on the surface for a significantly shorter period of time (typically <30 min). In these cases, **x**_{asc} and **x**_{des} for each surface drift were both taken as the sole position fix, and *t*_{des} and *t*_{asc} were set to the corresponding time.

The mean of the *p* measurements recorded along the drift track gave the measured *p*_{park} for each cycle. If no in situ *p* values were reported, *p*_{park} was assumed to be the preprogrammed nominal *p*_{park}. Only cycles for which 780 ≤ *p*_{park} ≤ 2020 db were included in this dataset. Data were manually inspected if the difference between the nominal *p*_{park} and the measured *p*_{park} exceeded 100 db. When it appeared that the float had shoaled, the pressure sensor had failed, or the float had come into contact with the seafloor (as determined by comparing the float trajectory with available bathymetric maps), the affected data were discarded. Any cycles flagged by the GDACs for likely grounding of the float were also removed from the dataset.

After estimating *p*_{park} and **u**(*p*_{park}), any outlying velocity estimates more than 10 standard deviations from the mean over the whole dataset were excluded from further analysis. To combine the velocity data with the dynamic height data, the **u**(*p*_{park}) estimates were adjusted to the common reference pressure *p*_{0} (900 db). The geostrophic shear between *p*_{park} and *p*_{0} was calculated by interpolating the mapped *D* fields (described below in section 4) to the location of **u**(*p*_{park}). The **u**(*p*_{0}) was then computed by adding this shear to **u**(*p*_{park}). Data for which the magnitude of the velocity adjustment was greater than 20 cm s^{−1} were eliminated from the dataset. The discarded data were all located within 5° of the equator, where the assumption of geostrophy becomes unreliable as the Rossby number grows as 1/*f*.

The resulting velocity dataset, consisting of 442 211 estimates of **u**(*p*_{0}) from 4799 floats, provides a representative sample of the geostrophic velocity at *p*_{0}. The distribution in space (not shown) and time (Fig. 1b) is comparable to that of the profile data, although there are slightly fewer estimates due to the differences in the quality and requirements of each dataset.

### b. Wind stress data

Wind stress data needed to compute the Sverdrup and Ekman transports in (5) were derived from observations from the National Aeronautics and Space Administration (NASA) SeaWinds scatterometer onboard the Quick Scatterometer (QuikSCAT). These data were obtained from the Centre de Recherche et d'Exploitation Satellitaire (CERSAT) at IFREMER in Plouzané (France). A gridded product was used that gave global monthly 10-m wind fields at 0.5° × 0.5° resolution, spanning from August 1999 to October 2009. Wind stress magnitude and direction as well as wind stress curl were provided.

## 4. Methods

### a. Geostrophic velocity

Using the Argo dynamic height and velocity data, monthly absolute geostrophic velocities in the Atlantic, Pacific, Indian, and Southern Oceans were computed on 29 pressure surfaces spanning the upper 2000 db for the period from December 2004 to November 2010. The procedure consisted of three basic steps. First, the *D* data were gridded on each *p* surface to produce 1° × 1° maps of both *D* and velocity relative to *p*_{0 }**u**_{rel}, which is derived from *D* according to

Second, the **u**(*p*_{0}) estimates were mapped onto the same spatial and temporal grid. A geostrophic streamfunction *ψ*_{0} was computed simultaneously, where *ψ*_{0} is defined here as

Third, the simple addition of the two fields gave monthly gridded absolute geostrophic velocities at each *p* level, as shown by summing (7) and (8). Similarly, the streamfunction for the absolute geostrophic flow *ψ*(*p*) was found by combining *ψ*_{0} and *D*(*p*) at each level.

All mapping was carried out using a variation of the method put forth by LeTraon (1990). Essentially this technique maps a large-scale signal using the objective function fitting procedure developed by Davis (1985) and a small-scale signal using the conventional objective analysis described by Bretherton et al. (1976). The two signals are characterized by different scales in both space and time. This procedure has the advantage of requiring no a priori knowledge of the statistics of the large-scale field; however, the statistics of the small-scale signal must be specified. Rather than simply assuming a decorrelation length scale, as is typically done, the present study uses a novel iterative approach to determine the statistics of the small-scale signal directly from the data. Full details of the mapping procedure including the iterative scheme are given in the appendix. For all *p* levels, mapping was restricted to areas where the depth of the seafloor *H* ≥ 1000 m, as determined from a box average of 1 arc-minute global relief model of Earth's surface (ETOPO1) bathymetric data (Amante and Eakins 2009). Additional grid points were excluded for *p* levels >1000 db if *p* ≥ *H*.

Uncertainties in the gridded velocities arise from a number of sources including array bias, measurement errors, and sampling errors (Katsumata and Yoshinari 2010). The mapping procedure used here has the distinct advantage of producing at each grid point a quantitative estimate of the mean squared error *e*^{2}. This error estimate combines the usual errors in the large- and small-scale signals with the error arising in the small-scale field due to uncertainty in the large-scale signal (LeTraon 1990). The total error in the absolute velocity was given by the sum of the *e*^{2} associated with the relative and reference velocities.

The gridded absolute velocities described here [Absolute Geostrophic Velocities from Argo (AGVA)], together with accompanying error estimates, are available for use by the oceanographic community and can be obtained online (http://flux.ocean.washington.edu/agva).

### b. Density

Assuming hydrostatic balance, volume transport can be written as the integral of velocity between two pressure surfaces *p*_{1} and *p*_{2} (with *p*_{1} < *p*_{2}):

where *g* is the acceleration due to gravity. Because this expression requires estimates of in situ *ρ* in addition to **u**, monthly maps of *θ*, *S*, and *p* on *σ*_{θ} surfaces were produced using the iterative two-step mapping technique presented above. These quantities were mapped on *σ*_{θ} surfaces instead of *p* surfaces due to the significant spurious anomalies that are introduced by averaging *θ* and *S* on *p* surfaces (Lozier et al. 1994). Mixed layer *p*, *θ*, and *S* were also mapped with the same method. The gridded *θ*, *S*, and *p* on *σ*_{θ} surfaces, together with the mixed layer maps, were used to compute monthly 1° × 1° maps of *ρ* and *σ*_{θ} as a function of *p*. Mean squared errors were also calculated for each mapped quantity. The gridded *θ* and *S*, vertically interpolated onto the same *p* levels as the **u** fields, are also available at the website listed above.

### c. Sverdrup balance

To assess the extent to which SB accounts for global meridional geostrophic transport, the two sides of (5) were computed. There is some question as to what time period is appropriate for testing SB. In other words, does SB hold over synoptic time scales or does it only apply when the wind stress and geostrophic velocities are averaged over a longer time scale? Sverdrup (1947) and most oceanographic textbooks imply that SB should emerge after averaging over a long interval but do not specify a particular length. Instead of averaging, some authors (e.g., Niiler and Koblinksy 1985; Vivier et al. 1999) chose to discuss a time-dependent SB. Similar to Wunsch (2011), this study will examine the longest averaging periods permitted with the datasets used here, namely 6 years for the Argo data and 10 years for the QuikSCAT data, with an overlap between the two of almost 5 years.

The wind-derived transport in (5) was calculated using the zonal wind stress and wind stress curl derived from QuikSCAT together with the mean mixed layer *ρ* computed from the gridded mixed layer *θ* and *S*. Monthly maps of wind-derived transport were produced for the entire QuikSCAT era (from August 1999 to October 2009) and then interpolated onto the same 1° × 1° horizontal grid as the mapped Argo data. These were then averaged over the whole time period and over the period contemporaneous with the Argo data (from December 2004 to October 2009).

To evaluate the geostrophic transport in (5), the mapped *υ*_{g} and *ρ* were used to compute *V*_{g} according to (9). In practice the integration was carried out by first finding the mean over the study period (from December 2004 to November 2010) of *υ*_{g}(*p*), *θ*(*σ*_{θ}), *S*(*σ*_{θ}), *p*(*σ*_{θ}), and the mixed layer properties. The mean was also calculated for the period overlapping with the QuikSCAT data. Using these mean values, was constructed for each grid point (where the angle brackets denote the mean). Next and were interpolated to common *p* levels with piecewise cubic Hermite polynomials. The resulting interpolants were used to calculate , which was then integrated analytically. While in principal is not exactly equal to , the difference here is considered negligible as .

Of critical importance to any analysis of SB is the choice of *h*, the bottom boundary of the transport integral in (2). One approach is to inquire as to the expected depth of the wind-driven circulation. The ventilated thermocline theory of Luyten et al. (1983) suggests that the wind-driven transport extends to the densest outcropping isopycnal in a given basin. Rhines and Young (1982) identify the depth of potential vorticity homogenization as the base of the wind-driven circulation. Motivated by these considerations and the observations of Talley (1988), Hautala et al. (1994) selected the depth of the 26.5 and 27.0 isopycnals in the eastern and central subtropical North Pacific, respectively. Other authors have integrated velocity from the surface to a fixed depth (e.g., Leetmaa et al. 1977). The most straightforward idea, namely, that *h* should be the depth where *w* goes to zero, arises directly from the derivation of SB. Although historical observations of *w* are generally unsuitable for practical application of this concept, the ocean state estimate used by Wunsch (2011) did allow *h* to be defined in this way.

Here, we have evaluated SB using two different methods of determining *h*. The first choice for *h* is the depth of the maximum mean isopycnal outcropping in wintertime along the line of zero Sverdrup transport, which was prompted by the ventilated thermocline theory. In the second approach, instead of selecting an a priori value for *h*, at each grid point the depth and corresponding pressure were found that gave the best fit to SB. If the meridional geostrophic transport agrees with the transport predicted by SB, then this depth should also be the location in the water column where *w* is approximately zero. Although the existence of such a depth does not by itself guarantee that SB is satisfied, this value will be used as a starting point for further analysis. In both cases, the mean maps of *θ*, *S*, *p*, and mixed layer properties were used to determine the desired limits of integration.

The uncertainty associated with the computed geostrophic transport was derived from the error estimates for *υ*_{g}, *θ*, *S*, and *p* using a Monte Carlo simulation with 500 iterations. The standard error around the mean was calculated for each quantity using the uncertainties on the monthly maps. In each iteration, , , , , and the mean mixed layer properties were individually perturbed by an amount randomly drawn from a uniform distribution on the interval (−*e*, *e*), and the transport was calculated as described above. The error in was decomposed into the error associated with and that associated with , and these were varied separately. At each grid point, the errors in each variable were assumed to be correlated vertically but independent of the errors in the other quantities. The variance of the resulting series of 500 transport values provided an uncertainty on the computed geostrophic transport.

## 5. Results and discussion

### a. Absolute geostrophic velocities

We first present the mean absolute velocity fields at four representative pressure levels: 5, 200, 1000, and 1500 db (Figs. 2–5). In each of these plots, *ψ*(*p*) is contoured at the intervals specified in the captions. Colors indicate highs in *ψ*(*p*) in red and lows in blue; as (8) shows, these represent anticyclonic and cyclonic flows, respectively. All estimates at latitudes ≤5°N/S are discarded because the geostrophic assumption becomes problematic in this region. In addition, results are only plotted where the uncertainty in the geostrophic velocity *e*_{|u|} < *σ*_{|u|}/8, where *σ*_{|u|} is the standard deviation of |**u**|.

At the shallow, near-surface level (Fig. 2), the dominant features are the anticyclonic subtropical gyres found in all basins except the northern Indian, the eastward-flowing Antarctic Circumpolar Current (ACC) south of approximately 40°S, and the cyclonic subpolar gyres in the North Pacific and North Atlantic. Eastward flow is observed in the northern tropical Pacific and Atlantic, associated with the North Equatorial Countercurrent in each basin. While the long-term mean flow in the northern Indian Ocean is weak, an analysis of the seasonal cycle (not shown) reveals a strong signal likely associated with monsoonal forcing. At this *p* level, the geostrophic flow shown here makes up only a part of the total velocity field, with Ekman flow constituting a significant portion, especially at lower latitudes.

At 200 db (Fig. 3), the subtropical gyres have all shifted poleward compared to their locations at 5 db. In the South Pacific, the subtropical gyre now has two distinct lobes, one to the east of New Zealand and one at about 25°S. Strong velocities are still found in the western boundary current regions and the ACC, while the circulation in the tropical regions is substantially weaker compared to the near-surface level.

In the mean flow at 1000 db (Fig. 4), the North Pacific subtropical gyre has become smaller and weaker relative to the upper levels. A generally southward flow is observed in the North Atlantic subtropics. The subpolar gyres, while somewhat reduced in strength, continue to occupy the entire North Atlantic and North Pacific north of approximately 45°N. This observation agrees well with the fact that the flow at higher latitudes is known to be less baroclinic than the flow in the subtropics (Talley et al. 2011). In the Southern Hemisphere, the subtropical gyres are shifted even farther poleward, with a strong recirculation east of southern Africa associated with the Agulhas retroflection. Anticyclonic flows in the western South Pacific, south of Australia, and in the southern Indian Ocean appear connected, indicative of a supergyre structure (Ridgway and Dunn 2007). The ACC continues to be quite vigorous at this pressure level, but elsewhere mean flow is relatively weak.

Last, the mean velocity field at 1500 db (Fig. 5) is similar in structure to the flow at 1000 db but exhibits more small-scale features, probably influenced by the smaller signal-to-noise ratio here. The increased southward flow along the western edge of the North Atlantic is suggestive of a deep western boundary current. A relatively strong southward flow along the east coast of South America is also present. In the Southern Hemisphere the ACC dominates the flow at this level, with significant anticyclonic flows persisting east of southern Africa and New Zealand and in the Tasman Sea.

The results presented here are broadly consistent with modern ideas of the large-scale circulation of the upper ocean. The most comprehensive observation-based description of the subsurface circulation up to this point was given by Reid (1989, 1994, 1997, 2003), which relied on decades of hydrographic and tracer data. Contrasting our results with these syntheses shows a high degree of similarity at the upper pressure levels in many areas of the global ocean, especially the North Pacific and North Atlantic. In the Southern Hemisphere, where the historical dataset was sparser, more differences between our results and Reid’s are apparent. For example, several small-scale undulations in the Agulhas retroflection and in the Pacific sector of the ACC (between 140° and 165°E and also between 160° and 120°W) stand out in our maps, where Reid’s analyses show purely zonal flows. Deeper in the water column, more substantial disagreement is found. For instance, at 1000 and 1500 db the Southern Hemisphere subtropical gyres are oriented predominantly zonally in our results, whereas Reid gives a stronger meridional component at these pressures. In general, the results from Reid show a smoother and more coherent flow at the deeper levels than is seen in the maps computed from Argo observations.

One might expect there to be dissimilarities between Reid’s results and the circulation estimates in this study due to the difference in the time periods over which the data were collected. While our estimates span from December 2004 to November 2010, the data compiled by Reid for his studies were collected from 1932 to 1991. Another major difference is the amount of data used in this study (over 480 000 Argo profiles) compared to Reid (5932 ship-based stations). In addition, our results employ a direct estimate of the absolute velocities at the reference level, whereas Reid relied on a subjective examination of tracer data to estimate absolute flow.

### b. Sverdrup balance

The wind-derived transport [*V*_{Sv} − *V*_{E} in (5)], averaged over the entire QuikSCAT period, exhibits the expected large-scale pattern of equatorward flow in the subtropics and tropics and poleward flow in the high latitudes (Fig. 6). The line of zero wind-derived transport is approximately aligned with the axes of the Kuroshio and Gulf Stream Extensions in the Northern Hemisphere and with the ACC in the Southern Hemisphere. The largest mean wind-derived transports are located in the higher latitudes, but there are local maxima in the eastern parts of most gyres. All of the following results will use wind transport averaged from August 1999 to October 2009 and Argo velocities averaged from December 2004 to November 2010. The longer averaging periods were adopted because they offered a slight reduction in the uncertainties and using the shorter contemporaneous periods produced only negligible differences.

The results of our first test of SB, using a priori estimates of the bottom boundary of the wind-driven circulation, are shown in Fig. 7. The maximum wintertime mixed layer *σ*_{θ} values along the line of zero Sverdrup transport, computed from the mapped mixed layer *θ* and *S*, are 26.24, 27.24, and 27.25 kg m^{−3} for the North Pacific, Southern Hemisphere and northern Indian, and North Atlantic, respectively. Because the zero Sverdrup transport line in the Southern Hemisphere was located in the circumpolar ACC, the same value was used for all basins there. The normalized difference between the two sides of (5) was determined as

with *V*_{g} computed from the mean Argo-derived velocity and density fields using the depth of the given isopycnals for *h* and *V*_{Sv} − *V*_{E} computed from the QuikSCAT wind stress fields. The difference Δ defined in this way is zero where SB is exactly correct and is bounded by ±1. Figure 7 gives the minimum Δ after factoring in the uncertainty on *V*_{g}. Yellow indicates perfect agreement, while those areas where the designated isopycnal was not found in the mean are shown in dark gray. Overlaid on the normalized difference are contours of *ψ* at 5 db.

Using these three isopycnals, we infer that SB provides a quantitative description of meridional transport for a large fraction of the global ocean, primarily in the subtropics and tropics. However, SB fails to accurately represent meridional transports in almost all of the ACC and subpolar North Pacific and North Atlantic for this choice of *h*. In the North Pacific subpolar gyre, the wind-derived transport predicts more poleward flow than is found in the geostrophic transport, which is perhaps not too surprising given that *σ*_{θ} 26.24 was chosen based on a ventilated thermocline model of the subtropical gyre. Using *σ*_{θ} 27.24, the northward geostrophic transport in the North Atlantic Current region is larger than the wind-derived transport, but north of this area the discrepancy reverses (at least in the part of the subpolar gyre where this isopycnal is present). For the most part, we find more poleward geostrophic transport than wind-derived transport in the northern part of the ACC and less in the southern part. Furthermore, meridional transports in the western boundary areas and the boundary current extensions do not adhere to SB. The poleward wind-derived transport is generally much smaller than the observed geostrophic transport in these regions. Last, the eastern boundary regions in the Pacific and Atlantic also exhibit poor agreement with SB.

Examining SB in a different way, we find the minimum *p* that produces agreement between the two sides of (5), within one standard deviation of *V*_{g} (Fig. 8). Regions where no *p* was found that satisfied SB, within the depth range of these velocity fields (2000 db), are shown in dark gray. While the existence of a *p* that satisfies (5) does not prove per se that SB is valid, it does indicate areas that are suitable candidates for SB. Although the results are somewhat noisy, the same picture emerges as in the first test: potentially good agreement with SB in the subtropics and tropics and poorer agreement elsewhere. Where agreement is possible, the wind-driven circulation appears to be confined to shallower pressure levels in lower latitudes than in higher latitudes, although small-scale features are prevalent, as expected for an average over 6 years (Wunsch 2010).

Even if there exists a value for *h* that balances the two sides of (5), is that enough to claim that SB is accurate on anything more than a purely mathematical level? Sverdrup’s original derivation assumes a flat bottom ocean, entirely wind-driven, above a quiescent abyss. While in reality the deep velocity certainly is not zero due to local thermohaline and topographic forcing, we expect that for SB to be valid, the wind-driven circulation must be sufficiently isolated from the deep circulation. To ensure this, we impose as an additional condition needed to satisfy SB that *υ*_{g}, within ±*e*_{υ}, is zero for some depth range below *h*. Because *υ _{g}* is related to ∂

*w*/∂

*z*through (1) and

*w*(

*z*=

*h*) = 0 when SB is satisfied, this criterion guarantees that

*w*≈ 0 in this depth range and thus that the wind-driven layer is well separated from the deeper circulation. The range chosen here corresponded to 200 db, which is arbitrary, although the results changed insignificantly for any choice from about 150 to 500 db. Stipulating this additional condition substantially decreases the area where SB provides an accurate description of the meridional geostrophic transport (Fig. 9). For the most part, the regions where SB fails under this stricter definition are in higher latitudes and western boundary regions, as well as areas in the central North Atlantic and Indian subtropical gyres.

These findings align reasonably well with the results of Wunsch (2011), who used a 16-yr average of data-assimilated model output to test the form of SB given by (2), with *w*_{E} given by *w*(*z* = −117 m). In that analysis, roughly 40% of the global ocean was found to agree with SB. The agreement was largely restricted to the interior subtropical gyres and the tropics, similar to what is presented here, although differences between those results and ours are apparent in the small-scale features. Considering the different datasets and methods of testing SB, the overall agreement suggests that these results are robust.

Given the evaluation of SB offered by Figs. 7–9, the question that arises is why does SB succeed in some regions and not others? One obvious reason for an apparent mismatch between the observed geostrophic and wind-derived transports in this study is that *w* may reach zero at a depth deeper than 2000 db, the limit of the geostrophic velocity dataset used here. Similarly, if *w* never actually vanishes, (2) is typically integrated from *z* = 0 to the seafloor so that a kinematic constraint on the flow at the bottom, *w*(*z* = *H*) = **u** · **∇***H*, can be used. This is likely a contributing factor to the lack of agreement in the high latitudes and the ACC, where the flow is more barotropic. At the points in these regions where SB is found to be valid, *h* is typically closer to 2000 db than is generally the case in the subtropics and tropics.

One potential reason that *w* might not become sufficiently small is that there are strong vertical motions associated with flow over topography. A simple test comparing the mean depth and bottom slope between the grid points where SB succeeds, as defined by the test shown in Fig. 9, and those where no agreement is found results in statistically significant differences (*p* value < 0.001) between the two groups. Places where SB is valid are on average deeper and flatter than those where SB is inaccurate. While this result does not demonstrate causality, it does validate our notion that SB is most accurate for a deep, flat bottom ocean, much as Sverdrup originally envisioned.

Another possible reason for disagreement in the high latitudes is that this study uses velocities averaged over 6 years, whereas the transit time for a first baroclinic Rossby wave is much longer in the high latitudes due to the reduction in the Rossby deformation radius (Chelton et al. 1998). While it is unknown whether a considerably longer time average would improve SB in the high latitudes, given the variability of the wind stress and other factors, the time needed for linear baroclinic adjustment there would certainly preclude the 6-yr time period used here. Furthermore, many high-latitude currents, such as the North Atlantic Current, are forced by both wind and buoyancy gradients. In these cases, we would expect that a separate wind-driven layer cannot reasonably be isolated from other elements of the circulation. Last, in the western boundary currents and their extensions, nonlinear dynamics likely play an important role. Because SB is derived from the linear vorticity balance given in (1), which neglects such effects, the lack of agreement in these regions is hardly surprising.

## 6. Conclusions

In this study, we have exploited the fact that the Argo array of profiling floats provides an exceptional tool for the estimation of large-scale absolute geostrophic velocities in the upper 2000 db of the ocean with unparalleled spatial and temporal resolution. The maps of global mean geostrophic velocities presented here correspond well with what has been previously reported based on traditional hydrographic data in the upper depth range, while at deeper pressures a more variable and complicated flow is observed. Limitations of the gridded velocity fields stem from the sparse coverage of the Argo array in some regions, specifically the inability to resolve boundary currents. As a purely data-based analysis, these velocity fields make an important contribution to our knowledge of the detailed circulation of the global upper ocean and should be useful for initializing and validating models and state estimations. These fields can also provide the foundation for observational analyses of meridional heat and freshwater transport and help further our understanding of the role of ocean circulation in influencing these important climatic quantities. As long as the Argo array is maintained, these fields can be continually augmented with more data and expanded to cover a longer time period.

The absolute geostrophic velocities presented here have been used to assess the extent to which Sverdrup balance, a simple but ubiquitous theory of the relationship between wind forcing and ocean circulation, accurately predicts observed geostrophic transports on a point-by-point basis. We find that within the uncertainties there is good agreement over large regions of the world’s oceans, namely, the subtropics and tropics away from the boundaries. Thus, for a substantial part of the global ocean, Sverdrup balance provides a reasonable quantitative picture of the observed circulation. However, the depth of the wind-driven circulation, defined as the location where the wind-derived transport matches the meridional geostrophic transport, has considerable geographic variability. While the simple model described by the Sverdrup balance is a useful starting point for a theoretical treatment of the large-scale circulation, the results presented here imply that it should not be unreservedly accepted as truth.

## Acknowledgments

The data used in this work were collected by the International Argo Program and the many national programs that have been involved. We thank each of the many agencies and individuals that have contributed to Argo over the past decade. In the United States, this work was supported by the National Oceanographic and Atmospheric Administration (NOAA) through Grant NA17RJ1232 Task 2 to the Joint Institute for the Study of the Ocean and Atmosphere at the University of Washington. We are especially indebted to Dr. Steve Piotrowicz of NOAA for his able and enthusiastic management of the U.S. Argo program since its inception. In addition, we thank the two reviewers for their insightful and helpful comments.

### APPENDIX

#### Objective Mapping of Argo Data

The technique used in this study to map the quantities of interest [*D*, **u**(*p*_{0}), *θ*, *S*, *p*, and mixed layer properties] closely follows that developed by LeTraon (1990). This method estimates both large- and small-scale signals from scattered data. The large-scale component, for which statistics are not known and perhaps not even clearly definable, is computed using the objective function fitting procedure of Davis (1985, 2005). A set of functions with spatial scales greater than a given limit, chosen from a complete set of functions, is fit to the data, subject to a zero-bias constraint. The large-scale signal is then removed from the data to produce small-scale anomalies that are objectively mapped following the conventional techniques of Bretherton et al. (1976).

In this study, the large- and small-scale signals were defined by different spatial and temporal scales. Spatially, the basis functions chosen to represent the large-scale field were spherical harmonics up to order 15. Thus, the large-scale component was the part of the signal with scales ≳24° latitude/longitude. This order was chosen for the spherical harmonic basis functions to accurately depict the large-scale signal and avoid overlap with the mesoscale field. Temporally, the large-scale signal was computed as a mean monthly seasonal cycle over the 6-yr study period. All data from each month, regardless of year, were used to compute a large-scale field according to the protocol outlined by LeTraon (1990) for nonsimultaneous data, whereby the covariance between data from different years was set to zero. The small-scale signal, in contrast, was mapped separately for each individual month. It is important to note that separation into large- and small-scale signals is strictly statistical and thus does not necessarily correspond to distinct physical processes.

The statistics of the small-scale signal, that is, the autocovariance of the anomalies, must be specified a priori and are used in estimating both the small-scale field and the large-scale field (for which the least squares fit to the given basis functions takes into account the correlation of the anomalies). A Gaussian form for the spatial autocovariance was used here, although other forms were tried and did not substantially affect the results. The autocovariance was assumed to depend only on separation distance *s* according to

for variance *C*_{0} and decorrelation length scale *L*.

While typical analyses chose a value for *L* somewhat arbitrarily, we employ a novel method that determines *C*_{0} and *L* in an iterative fashion directly from the data for each *p* or *σ*_{θ} level. An initial guess for *C*_{0} and *L* was used to first estimate the large-scale signal via objective function fitting and remove that signal from the data. Here, we used *C*_{0} = 1 and *L* = 250 km, but testing showed that the final parameters were insensitive to the choice of the initial values. The resulting anomalies were then used to calculate the observed autocovariance using all pairs of data collected in the same month and separated by spatial lag *s*:

where the angle brackets denote an ensemble mean taken over all valid pairs of data, and the prime denotes anomalies from the large-scale signal. The expression (A2) is written for *D*, but analogous equations were used for all scalar properties being mapped. All data pairs were grouped according to separation distance into 25-km bins and averaged to give a mean autocovariance curve. Next, a least squares fit of the observed mean autocovariance curve to (A1) gave both *C*_{0} and *L*. The data at the smallest separation distances (≤25 km) were excluded from the fit due to probable space–time aliasing. Using the fitted autocovariance parameters, the mapping of the large-scale signal and calculation of the autocovariance were then repeated. This procedure was iterated until the parameter values converged. This novel iterative approach has been tested and validated using both analytically constructed fields and output from a high-resolution model, the results of which can be found in A. R. Gray and S. C. Riser (2014, unpublished manuscript).

To maximize the number of data pairs included in the autocovariance estimation, and hence the degrees of freedom, the autocovariance is assumed to be isotropic. While this is not strictly true, Davis (1998) found that in the South Pacific the variability is generally isotropic except within about 5° of the equator. Thus, the assumption of isotropy most likely introduces significant errors only in the equatorial region, where the validity of the geostrophic relationship is already questionable.

Once the covariance parameters were determined for each *p* or *σ*_{θ} level, the large- and small-scale signals were both mapped onto a 1° × 1° horizontal grid at the above-mentioned temporal resolution. Observational error and subgrid-scale noise were specified to account for an additional 50% of the measured variance in the data; this value was chosen to be reasonably conservative, following Gille (2003).

To map **u**_{rel} and **u**(*p*_{0}), slight adjustments were made to the basic procedure. LeTraon (1990) and Bretherton et al. (1976) outline methods for mapping multivariate data and for determining fields related to the data through linear operators. As (7) and (8) show, the relative and reference components of velocity are related linearly to *D* and *ψ*_{0}, respectively. Using partial derivatives of both the autocovariance function (A1) and the spherical harmonic basis functions allowed us to map **u**_{rel} from the *D* data and *ψ*_{0} from the **u**(*p*_{0}) data. In addition, because the **u**(*p*_{0}) data were vectors, a simple application of (A2) was not possible. These data were first multiplied by *f* to account for the latitudinal dependence of (8). Then for each pair of velocities the vectors were decomposed into longitudinal and transverse components, which are aligned parallel and perpendicular to the axis between the two data points, respectively. Longitudinal and transverse autocovariance curves were computed following Batchelor (1970), who showed that these can be derived from the covariance function describing the velocity streamfunction given that the statistics of the field are isotropic, homogeneous, and stationary. Last, a least squares method was used to simultaneously fit the mean observed longitudinal and transverse curves to the analytical forms of those functions.

## REFERENCES

*J. Mar. Res.,*

**40,**559–596.

## Footnotes

A comment/reply has been published regarding this article and can be found at http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-14-0127.1 and http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-14-0215.1