## Abstract

This study investigates the observed spectral character of eddy heat fluxes near the ocean surface, focusing on the distribution in wavenumber and phase speed space. Eddy heat fluxes in the eastern Pacific are calculated from concurrent satellite sea surface height and sea surface temperature data. A high-resolution coupled climate model is also analyzed in order to verify the physical mechanisms involved and to validate the model against observations. Wavenumber, frequency, and phase speed power spectra and cross spectra are constructed and presented as a function of latitude. These spectra reveal the dominance of coherent mesoscale eddies in both the length scale and phase speed of eddy heat fluxes. The breadths of the spectra are characterized via spectral moments; these moments show that the eddy fluxes are relatively concentrated around the dominant wavenumber and phase speed. Good agreement is found between the model and the observed spectra. The integrated heat transport and corresponding eddy diffusivity are shown to compare well with previous studies, but the results give a deeper insight into what determines the heat flux. Implications for eddy parameterization are discussed.

## 1. Introduction

Transient motions (also known as “eddies”) in the ocean and atmosphere drive significant transport of mass and tracers. Of particular importance is the meridional eddy heat transport, which contributes to the maintenance of Earth’s pole-to-equator temperature gradient (Trenberth and Caron 2001; Wunsch 2005). Although eddy heat fluxes in the ocean are relatively less significant than in the atmosphere, they are still an important part of the ocean heat budget, particularly at regional scales and in the Southern Ocean (Jayne and Marotzke 2002; Volkov et al. 2008; Hausmann and Czaja 2012; Abernathey and Cessi 2014). Because of their relatively small spatial scales, ocean eddy fluxes are more difficult to observe than those in the atmosphere, and their statistical properties are less well characterized. Satellites provide a uniquely powerful tool for observing eddies at the ocean surface. Climate models are also beginning to resolve the ocean mesoscale (McClean et al. 2011), and comparison between remotely sensed observations and high-resolution models is a crucial form of model validation.

A fundamental question is what determines the strength of the eddy flux and how this flux is related to more readily observable eddy properties such as eddy size, kinetic energy, and so on. Inspired by the classical “mixing length” arguments of Taylor (1915) and Prandtl (1925) regarding turbulent fluxes, many studies have assumed the eddy flux in the ocean to be proportional to the background tracer gradient (i.e., that it is diffusive) and to the product of a characteristic eddy size and eddy velocity (e.g., Holloway 1986; Keffer and Holloway 1988; Held and Larichev 1996; Visbeck et al. 1997; Stammer 1998). More recent studies have added a new ingredient to the equation: the eddy propagation relative to the background mean flow (Marshall et al. 2006; Smith and Marshall 2009; Abernathey et al. 2010; Ferrari and Nikurashin 2010; Klocker et al. 2012a,b; Abernathey and Marshall 2013). In particular, the simple stochastic model of Ferrari and Nikurashin (2010) demonstrates how zonal phase propagation suppresses meridional eddy diffusion and puts forth a quantitative theory for the magnitude of this effect.

The framework of Ferrari and Nikurashin (2010) was recently tested by Klocker and Abernathey (2014, hereinafter KA14) in a comprehensive way using kinematic passive tracer simulations in the east Pacific (the same sector studied here). Those results indicate that the extratropical, meridional eddy flux of a passive tracer due to mesoscale eddy stirring could be parameterized quite well in terms of a single wavenumber and phase speed at each latitude. The appropriate phase speed was found to be the long baroclinic Rossby wave phase speed, while the appropriate wavenumber was found to be proportional to the average diameter of tracked nonlinear coherent eddies from the eddy census of Chelton et al. (2011). The fact that the eddy flux can be parameterized in terms of an essentially monochromatic model seems at odds with the fact that the ocean contains a broad spectrum of variability in space and time (Richman et al. 1977; Stammer 1997; Hughes and Williams 2010; Wortham and Wunsch 2014). Therefore, a key motivation for our study is to attempt to reconcile the success of the monochromatic Ferrari and Nikurashin (2010) model with the broadband nature of the variability. We wish to assess how narrowly concentrated the eddy flux is around a single length scale and phase speed.

To answer this question, we employ satellite observations to investigate the spectral character of surface eddy meridional heat fluxes over a wide range of latitudes. This is achieved by calculating wavenumber–frequency cross spectra for sea surface temperature (SST) and the geostrophic velocity derived from sea surface height (SSH). There is an extensive literature on the analysis of spatiotemporal variance and covariance in different remotely sensed ocean surface datasets such as SSH, SST, and color-derived chlorophyll [see review by O’Brien et al. (2013)]. Many of these past studies focus on characterizing the propagation behavior of Rossby waves (Chelton and Schlax 1996; Polito and Cornillon 1997; Cipollini et al. 1997; Hill et al. 2000; Cipollini et al. 2001; Polito and Liu 2003; Killworth et al. 2004) and tropical instability waves (Polito et al. 2001; Contreras 2002; Chelton et al. 2000; Lee et al. 2012). The paper by Killworth et al. (2004) is particularly comprehensive and makes a convincing case that a large fraction of the variance in SST and surface chlorophyll arises by advective stirring by the surface geostrophic flow by motions propagating close to the long Rossby wave speed.

Being interested in the wave dynamics themselves, the studies cited above generally employed filters to isolate the spectral bands of interest. Here, the approach is slightly different: we consider the total, unfiltered eddy flux and examine its spectral density in wavenumber, frequency, and phase speed space. This perspective is inspired by the atmospheric study of Randel and Held (1991, hereinafter RH91), who made such a diagnosis for the eddy fluxes of heat and momentum in the troposphere. In particular, presenting the results of the cross-spectral analysis through 2D contour plots as a function of latitude and phase speed (or latitude and wavenumber) provides a novel view of the oceanographic data, revealing the strong latitudinal dependence in the spectra. Our results indicate that the extratropical meridional eddy heat flux in phase speed space is indeed concentrated around the long Rossby wave phase speed.^{1} Furthermore, the dominant length scale associated with the eddy heat flux is everywhere very close to the observed mesoscale eddy diameter. These conclusions help to explain the success of KA14.

In addition to analyzing the satellite data, we perform the same analysis on a state-of-the-art, global, eddy-resolving/eddy-permitting ocean model. This comparison serves two purposes. On one hand, it allows us to probe finer space and time scales than the observations can resolve. On the other, it provides a form of model validation in the spectral domain. The broad agreement between the model and the observations is encouraging on both fronts, suggesting that the observations are sufficient to resolve the dominant scales of eddy heat transport and that the model shares the spectral characteristics of the observations.

Our paper is organized as follows: Section 2 describes the satellite data and numerical model used to obtain SST and SSH. The basic results of the cross-spectral analysis are presented in section 3 as a function of latitude and wavenumber, latitude and frequency, and latitude and phase speed. In section 4, we calculate the first two moments of the spectra, diagnosing the dominant scales and also the breadth in wavenumber and phase speed space, and discuss the features observed. Section 5 examines the net meridional heat transport in the surface layer and compares the satellite results with the model and with previous estimates. Conclusions are given in section 6.

## 2. Data and models

To compute the meridional eddy heat flux, we need concurrent observations of meridional velocity and temperature. We focus our study on a sector in the east Pacific spanning 60°S to 50°N and ranging from 180° to 130°W. Beyond the intrinsic importance of the eastern Pacific for global climate variability (e.g., related to the El Niño–Southern Oscillation phenomenon), this sector was chosen specifically to facilitate comparison between our results and those of KA14 and Abernathey and Marshall (2013). Those earlier works picked this sector because it is relatively statistically homogeneous in longitude and contains very little land. These attributes are also well suited to our purposes here, which is to examine the dependence of the spectra on latitude. Understanding the spectral character of eddy fluxes in more inhomogeneous regions, such as western boundary currents, is an important problem but beyond the scope of this work.

### a. Sea surface height

Sea surface height data are used to estimate surface geostrophic velocities. The altimeter products were produced by SSALTO/DUACS and distributed by AVISO, with support from CNES (www.aviso.oceanobs.com/duacs/). For this study, we use the precomputed geostrophic velocities derived from the delayed time, two-satellite, “reference,” merged sea level anomaly fields. In these precomputed velocities, the method of Lagerloef et al. (1999) is applied in the equatorial band (±5°). This method, based on the “equatorial geostrophic” vorticity balance, has been validated with in situ current meters and allows us to obtain velocity estimates in this region. Nevertheless, we must maintain some skepticism of the results in the equatorial band; we focus primarily on the extratropics. A snapshot of the AVISO SSH field from this sector is shown in Fig. 1.

The horizontal spacing of the AVISO gridded data is ¼°. The effective resolution of the product is such that it “sees” eddies of approximately 50-km diameter and larger (Chelton et al. 2011); however, smoothing applied during the gridding procedure acts as a low-pass filter, attenuating the signal weakly at wavelengths below 200 km and strongly below 100 km (Ducet et al. 2000). Consequently, the SSH signal displays very little power at short wavelengths (see Fig. 3). This filtering would make it difficult to estimate, for example, spectral slopes characterizing turbulent inertial ranges from the gridded data. However, the focus here is not on the inertial range but of the prominent peak wherein most of the kinetic energy resides (Stammer 1997). This peak, which is everywhere at wavelengths larger than 200 km, is well resolved by the gridded data; mixing length arguments suggest that these large-scale, highly energetic motions should also dominate the heat flux (Larichev and Held 1995). Directly assessing the contribution of the filtered smaller scales to the heat flux could potentially be explored using along-track satellite or in situ data, but we do not take that route here. However, one motivation for examining the numerical model (with ° grid spacing, described below) is to attempt to probe smaller scales. In both model and data, we find the heat flux is dominated by wavelengths larger than 200 km, although the model shows a greater contribution from smaller scales.

AVISO produces a map every 7 days that represents a best estimate of the SSH field on that day. The data record begins in 1992, but we only consider the 9.3-yr period concurrent with the SST observations, as described below.

### b. Sea surface temperature

The SST data are the Group for High-Resolution Sea Surface Temperature (GHRSST) global level-4 sea surface temperature analysis produced by the NOAA National Climatic Data Center (Reynolds et al. 2007). An SST map is produced daily on a ¼° grid. We selected the version of the product that blends data from the 4-km Advanced Very High Resolution Radiometer (AVHRR), the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E), and in situ ship and buoy observations using optimal interpolation. The SST value represents the temperature at approximately 0.3-m depth. The data coverage for this product begins in June 2002 and ends in October 2011, the period of operation of the AMSR-E instrument.

### c. 1/10° POP model

The model analyzed here, a version of the Community Earth System Model (CESM), is described in Small et al. (2014). (The CESM code name for this run is hybrid_v5_rel04_BC5_ne120_t12_pop62.) This is a global climate simulation that includes ocean, atmosphere, sea ice, and land models, similar to that described by McClean et al. (2011). The ocean component, our focus here, uses the Parallel Ocean Program (POP) code and is hereinafter referred to as the “POP model.” The POP model ocean has a nominal grid spacing of 0.1° and an atmospheric grid spacing of 0.25°. The two are coupled every 6 h. This combination of high-resolution atmosphere, ocean, and coupling results in one of the most realistic global simulations currently available. We consider the use of a coupled model (rather than an ocean-only model) important for simulating the complex air–sea interactions that arise over mesoscale eddies (Small et al. 2008; Bryan et al. 2010), which may influence the eddy heat flux.

While the grid spacing of the POP ocean model is 0.1°, the effective resolution is not this fine. To characterize the effective resolution, we examined the wavenumber power spectra (described in the subsequent section) and found a modest high-wavenumber spectral rolloff at wavelengths shorter than 40 km. From this we conclude that the effective spatial resolution of the model is around 40 km, not significantly different from AVISO.

We extract daily surface velocity and SST fields from the same sector described above for a 5-yr period (model years 46–50). This is a significantly higher temporal resolution than the AVISO data. While higher-frequency motions are clearly present in the model spectra, our analysis below indicates that these high-frequency (superweekly) components make a negligible contribution to the meridional heat flux.

Note that we do not attempt to isolate the geostrophic component of the flow in the model; while this would allow for a more direct comparison with the AVISO results, we prefer instead to examine the full eddy flux produced by the model. As we will see, the similarity between the model and satellite results suggests that the flux in the model is indeed dominated by geostrophic motions.

### d. Preprocessing

Relatively little additional processing is applied to the data since the observational products we have chosen are already highly processed. For the analysis of the satellite data, we applied temporal smoothing to the daily SST data through a 7-day boxcar and then subsampled the smoothed SST data on the same days as the weekly AVISO output. This procedure eliminates aliasing in the frequency domain. For all datasets, we subtract the time mean at each point in space (this is already done in the case of AVISO fields, which are the anomaly relative to the 1993–99 mean) and then the zonal mean at each time step, effectively removing the basin-scale variability. Removing the zonal mean also filters out most, but not all, of the seasonal cycle. A small amount of smoothing in the frequency domain is also applied before interpolating to phase speed space (described further below). Everything remaining is included in our definition of eddy variability.

### e. Coherent eddy statistics

A large amount of the variance in midlatitude SSH has been attributed to coherent, nonlinear mesoscale eddies (Chelton et al. 2011). Throughout this study, we compare the length scales and phase speeds that arise from our spectral analysis with the coherent eddy characteristics from the eddy census of Chelton et al. (2011), whose results were made publicly available (http://cioss.coas.oregonstate.edu/eddies/). The observed eddy length scale *L*_{s} is defined by the average radius of all the eddies at each latitude in the sector. This radius itself is determined for each eddy from the area enclosed by the SSH contour corresponding to the maximum geostrophic flow speed, that is, where the eddy velocity is greatest. [See Chelton et al. (2011), their section 4.2, for further discussion of the eddy length computation.] To convert this length scale to a wavenumber, we follow the recommendation of Chelton et al. (2011) and assume the eddy streamfunction to be described by a Gaussian function with an *e*-folding scale of . We then define the corresponding eddy wavenumber as .

Two additional length scales are relevant to our study. The first baroclinic Rossby radius of deformation is a fundamental length scale for the dynamics of large-scale ocean circulation; of particular relevance here is the fact that the most unstable mode of baroclinic instability occurs near the deformation radius (Stammer 1997; Chelton et al. 1998; Scott and Wang 2005; Smith 2007). It also enters into the Rossby wave dispersion relation and exerts a strong control on the observed phase speeds. The deformation radii arise as eigenvalues of the vertical stretching operator in the quasigeostrophic potential vorticity for a resting ocean (Pedlosky 1987). Specifically, the first deformation wavenumber is the largest eigenvalue *K*_{d} given by the Sturm–Liouville equation:

where *f* is the Coriolis parameter, and *N* is the Brunt–Väisälä frequency. We obtained the deformation radius data from Tulloch et al. (2009). Our *K*_{d} values represent a zonal average over the sector. Given *K*_{d}, the long baroclinic Rossby wave phase speed in the zonal direction is then given by

where is the time- and depth-averaged zonal flow, and *β* is the meridional gradient of *f*. This expression represents the long-wave limit of the classical linear Rossby wave dispersion relation, with the addition of the Doppler shift by the depth-averaged flow. As argued by Klocker and Marshall (2014), the depth-averaged flow seems to be the relevant Doppler shifting for nonlinear mesoscale eddies. And as shown by KA14, poleward of 20° latitude, *c*_{R} ≃ *c*_{eddy} in this sector.

We also invoke the Rhines scale, defined as the wavenumber , where *u*_{rms} is the root-mean-square eddy velocity (calculated from the AVISO data). This represents the scale at which turbulent motions become effective at transferring energy into zonally elongated flow such as jets (Rhines 1975; Maltrud and Vallis 1991).

For a given wavenumber *K*, the wavelength is defined as *L* = 2*π*/*K*. While elementary, this conversion can be a source of great confusion. For example, under this terminology the Rossby deformation wavelength near 45°S is approximately 120 km (as in Tulloch et al. 2009); this is greater than the deformation radius of Chelton et al. (1998) by a factor of 2*π*. It is common to plot wavenumber spectra in terms of “cycles per meter,” which implies the division of wavenumber *K* by a factor of 2*π*. This quantity is best described as an “inverse wavelength,” not a wavenumber. Correct treatment of this issue is crucial, for example, in assessing the strength of the inverse cascade.

## 3. Cross-spectral analysis

Here, we describe the technical details of the spectral analysis and present the basic results. Extensive discussion of the results is deferred until section 4.

### a. Univariate power spectra for SST and surface meridional velocity anomalies

Here, we describe the calculation of wavenumber–frequency spectra for *θ*, the SST anomaly. An identical procedure applies to *υ*, the meridional velocity anomaly. (Note that because of the preprocessing described above, all variables are anomalies from the time and zonal mean.) In principle, *θ* is a continuous function of the zonal coordinate *x* and time *t*; *θ* = *θ*(*x*, *t*) at each latitude *φ* in the sector. However, our observations are discrete, with *N* spatial points in longitude (spaced by Δ*x*) and *M* points in time (spaced by Δ*t*) such that the total zonal length of the sector is *L* = *N*Δ*x* and the total temporal length of the record is *T* = *M*Δ*t*. The discrete space and time coordinates are denoted as *x*_{n} = *n*Δ*x* and *t*_{m} = *m*Δ*t*. For the satellite data used here, *L*(*φ*) = 2*πa* cos(*φ*)(50/360), where *φ* is latitude (giving the width of our 50°-wide sector), *N* = 200 (° spacing), *T* = 3402 days, and *M* = 486 (data every 7 days).

We write the discreet SST anomaly as

We can express *θ*_{mn} using a discrete inverse Fourier transform as

where are the complex Fourier components, is the wavenumber, and is the angular frequency. Equation (4) summarizes the normalization and unit conventions in our Fourier transform definitions. We adopt the convention of RH91 in which all wavenumbers are positive while frequencies take both positive and negative values. The values of are computed numerically from using the NumPy implementation of the fast Fourier transform (FFT) algorithm.

Parseval’s theorem states that the total power of the signal is the same in either basis. The normalization condition chosen in Eq. (4) means that each Fourier component represents a fraction of the variance, such that

where the asterisk denotes the complex conjugate, and the overbar is a sum over all wavenumbers/frequencies (here equivalent to a time and zonal mean).

We define the power density as a function of wavenumber as the sum over all frequencies:

The normalization by , the spacing of the discrete wavenumbers, means that represents a continuous power density function, giving results that are independent of *N*. Similarly, we define the power density as a function of frequency as the sum over wavenumbers:

where is the spacing of the discrete frequencies.

As in RH91, we construct phase speed spectra by interpolating the spectral density from space to space, where is the phase speed. Before interpolating, the density must be multiplied by *k*; this transformation ensures that the total power is the same whether integrating over *ω* or *c*. We also smooth the signal in frequency space before interpolating, using a Gaussian filter with an *e*-folding scale of two frequency bands to avoid aliasing. We interpolate to 1000 points in *c*, evenly spaced from −1 to 1 m s^{−1} in order to capture the wide range of phase speeds with high precision. This band contains 95.5% of the total SST variance and 99.3% of the velocity variance. After summing over wavenumbers, we obtain the power density as a function of phase speed .

Raw wavenumber–frequency power spectra at different locations in the ocean are shown in numerous other publications (e.g., Killworth et al. 1997; Wunsch 2010; Wortham and Wunsch 2014) and are not plotted here. Here, we are interested instead in the integrals of the spectra. In Fig. 2, we plot , , and , using a logarithmic color scale. This figure reveals the distribution of SST variance by wavenumber, frequency, and phase speed as a function of latitude.

From the surface meridional velocity data, we define (the velocity anomaly space/time coordinates) and (the Fourier transform) in the same way described above. Figure 3 shows , , and .

### b. Eddy heat flux cross spectra

Parseval’s theorem also applies to the product of *θ* and *υ*; the eddy heat flux is the same whether expressed as an average of space/time components or a sum of Fourier components. We express this mathematically as

Just as described above for the univariate spectra, we can sum the components of selectively to define , , and . These functions are plotted in Fig. 4. Unlike the power spectra described above, can take both positive and negative values, corresponding to northward and southward heat transport. As seen in Fig. 4, the eddy heat flux is poleward in both hemispheres, except near the equator, where it reverses. This is consistent with the mean SST gradient, which also reverses near the equator; the eddy flux is always downgradient. Further discussion and comparison of the spectra is deferred until section 4.

### c. Cross-spectral analysis of 1/10° POP model

The analysis of the model is identical. Only the spatiotemporal sampling and resolution are different. For the model *L*, the sector width is the same, *N* = 500, *T* = 1825 days, and *M* = 1825. The different spectra from the POP model are shown in Figs. 5, 6, and 7. To facilitate comparison in phase speed space, which is the primary focus of our study, in Fig. 8 we plot from both satellite data and the POP model in the same figure. Showing only the extratropics (poleward of 10°) in Fig. 8, where the magnitude of the phase speed is small compared to the equatorial region, allows us to examine the spectra in closer detail.

## 4. Spectral moments and discussion

While a visual comparison of the spectra in Figs. 1–7 is informative, a more quantitative comparison is desirable. In particular, we wish to assess whether the spectra have peaks in the same locations and how those peaks are related to the underlying dynamics. Furthermore, we wish to assess to what extent the spectra are narrowly concentrated around these peaks versus broadly distributed. In this section, we characterize the properties of the power spectra and cross spectra via their moments in wavenumber, frequency, and phase speed space. The first moment tells us about the dominant scale. The second moment tells us how concentrated the distribution is about that scale. Spectral moments play an important role in the theory of geostrophic turbulence and have frequently been used to characterize the length scales of ocean eddies (e.g., Rhines 1975; Scott and Wang 2005; Tulloch et al. 2011).

In our discrete notation, we define the first moment of a wavenumber spectrum at a given latitude [itself defined in Eq. (6)] as

The second moment is then defined as

We make analogous definitions for the *ω* and *c* moments and also for the moments of the other spectra and .

A Gaussian distribution is completely described by the mean *M*_{1} and the variance *M*_{2}. The interpretation of *M*_{1} and *M*_{2} is therefore clearest when the spectra have a clearly defined, dominant peak. In the case of , which is not positive definite, it is possible for the normalization factor in the denominator of Eqs. (9) and (10) to approach zero, leading to a meaningless result. To avoid this situation, we mask the moments at latitudes where the ratio . This is only the case in regions where the mean SST gradient is weak (primarily near 15°S) and the heat flux is vanishingly small and noisy. At most latitudes, the heat flux does display a clear spectral peak, as evident in Figs. 4 and 7.

The first moments , , and are plotted in the upper panel of Fig. 9, with the second moments , , and in the middle panel, for both the satellite data and the POP model. The *c* moments are shown in Fig. 10. Equipped with this more quantitative description, we are now prepared to discuss the features and characteristics of the spectra.

At all latitudes, the first moment of the eddy heat flux is at wavelengths of 250 km or greater (Fig. 9). Generally speaking, these are the same wavelengths containing most of the kinetic energy [as indicated by ]. The dominance of these wavelengths in the kinetic energy has been noted by other authors (e.g., Wunsch 2010; Wortham and Wunsch 2014) and has been attributed to a weak geostrophic, turbulent inverse cascade of energy from some source scale to the deformation radius (Stammer 1997; Arbic and Flierl 2004; Scott and Wang 2005; Tulloch et al. 2011). It is interesting to see that these wavelengths also dominate the heat flux. This implies that the temperature anomalies generated by eddy stirring are not strongly damped by interaction with the atmosphere.

To further analyze the dominant scales of the heat transport, the bottom panel of Fig. 9 plots the ratio of the three important wavenumbers identified in section 2 (the coherent eddy wavenumber *K*_{eddy}, the deformation wavenumber *K*_{d}, and the Rhines wavenumber *K*_{β}) and . This ratio is very close to one for *K*_{eddy} poleward of 20°. The consistent proportionality between the dominant length scale of the heat flux and the diameter of coherent mesoscale eddies strongly suggests that the mesoscale eddies are responsible for most of the heat transport. The deformation wavenumber is consistently larger than , especially at high latitudes, where the ratio approaches 4. The Rhines wavenumber is smaller.

Of course, when interpreting the wavenumber spectra, we must keep in mind the low-pass filtering effect of the AVISO processing, which artificially attenuates small scales. This filtering attenuates the inertial range at wavelengths below 200 km (Ducet et al. 2000), in which power-law dependence is observed in the along-track altimetry (Xu and Fu 2012). One probable consequence of this filtering is to bias toward larger length scales. This effect can be seen by comparing with the POP results. It is clear from Figs. 1–7 that there is significantly more small-scale variance in the POP fields. This is certainly reflected in larger values of (Fig. 9) for the POP model by a factor of 2–3 for all three spectra. The values of are also significantly different between the POP model and the data; a dominance of shorter wavelengths is seen in the model SST. Since the SST spectra are red throughout (in contrast with the energy spectra, which contain a clear peak), such a shift is expected. There is also a shift to shorter wavelengths in the POP and ; however, the magnitude of the shift is much smaller (generally less than 20%), and at many latitudes, the model and data coincide well. This indicates that the first moment of is primarily determined by the mesoscale peak in and only weakly affected by higher wavenumbers, where the energy typically falls off as a power law in *k*.

We now turn to the phase speed moments in Fig. 10. In general, the first moment of the eddy heat transport tracks in both the model and data, especially poleward of about 10° latitude. Moreover, this moment tracks quite closely the observed eddy propagation speed *c*_{eddy}, which itself is well described by the long Rossby wave phase speed *c*_{R} poleward of 20° latitude (Tulloch et al. 2009; KA14). This further suggests that the extratropical heat flux is generated through mixing by coherently propagating eddies. Closer to the equator the moments are noisier because of the regions of vanishing temperature gradient and heat flux. Nevertheless, the heat flux and energy spectral moments are consistent with tropical instability waves, which propagate at speeds of 0.3–0.4 m s^{−1} (see also Contreras 2002; Polito et al. 2001). In general, the first moments of the temperature spectra tend toward slower phase speeds in both model and data, reflecting a wider range of mechanisms that cause temperature variance.

The second moments show that the temperature variance (in both model and data) is spread much more widely across phase speeds than either the velocity variance or the heat flux. For the satellite data, both and have values of <0.05 m s^{−1} in the extratropics, indicating that the corresponding phase speed spectra are highly peaked about their first moment. This result provides an answer to one of our key motivating questions. Given this narrow peak, the success of the monochromatic model of KA14 in the extratropics seems unsurprising; there is clearly a single dominant phase speed in the data. However, the model moments indicate somewhat greater spread in phase speed, a fact not immediately evident in Fig. 8 on visual inspection. As with the wavenumber moments, is uniformly higher in the model than in the data. In contrast to the wavenumber spectra, the satellite temperature data actually contain more phase speed spread than the model.

## 5. Net meridional heat transport and eddy diffusivity

In this section, we consider the total meridional eddy heat transport at each latitude. This section has two goals. First, we demonstrate that our “direct” method of estimating near-surface heat transport from satellite data provides results that are consistent, in magnitude and spatial structure, with previous model-based estimates. Second, we demonstrate that the differences in heat transport between the POP model and satellite data can be explained through the lens of the Ferrari and Nikurashin (2010) diffusivity closure.

We calculate the net meridional eddy heat transport from both the satellite data and POP model. This heat transport is defined as

where *ρ*_{0} = 1027 kg m^{−3} is a reference density and *c*_{p} = 4186 J K^{−1} kg^{−1} is the specific heat of seawater at constant pressure. To translate the surface eddy temperature flux into a heat flux, it must be multiplied by a finite depth *h*. This choice of depth strongly influences the estimate of . We opt for a lower bound estimate.

We expect that the remotely sensed surface velocities and temperatures are representative of the mixed layer. Observational studies have shown the vertical amplitude of eddy heat transport in the subtropical North Pacific is at its maximum near the surface and decays to zero over a depth scale of several hundred meters (Roemmich and Gilson 2001; Qiu and Chen 2005). Since our surface observations cannot assess the subsurface component without further assumptions, we simply present them as a lower bound on the vertically integrated heat transport. Based on the Argo-derived mixed layer estimates of this sector from Holte and Talley (2009), we use a spatially constant value of *h* = 50 m. This is itself a lower bound; in the Southern Ocean, the annual-mean mixed layer is considerably deeper. Furthermore, mixed layer depth exhibits considerable temporal variability. However, a spatially constant value facilitates direct comparison between the model and the data, which is more important here than the precise magnitude of .

The value of is shown in Fig. 11 for both satellite data and POP model, normalized to give units of TW (10^{12} W) per degree of longitude. The order of magnitude in the extratropics is roughly 0.1–0.3 TW per degree or 10–30 TW across the roughly 100° width of the Pacific. For comparison, Jayne and Marotzke (2002) estimate 50–100 TW in the upper 25 m of the global ocean (their Fig. 3). Of course, the heat transport is not zonally homogeneous across the Pacific, with transport typically concentrated in the western boundary current regions. The point is that the present eddy heat transport estimates are of the same order of magnitude of other estimates (Jayne and Marotzke 2002; Volkov et al. 2008; Dong et al. 2014) and therefore represent a significant contribution to the total global meridional heat transport. The meridional structure is also similar to the studies cited.

The model and data generally agree quite well, both in spatial structure and the magnitude of . Disagreement in magnitude, by factors of 2–3, arises in two principal locations: near the equator and in the ACC. (This disagreement is also visible, though not quite as obvious, in comparing Figs. 4 and 7.) In these areas, the model produces a larger-eddy heat transport. What is the physical origin of this disagreement? The closure model of Ferrari and Nikurashin (2010) posits that the eddy flux depends on four principal factors: the background meridional tracer gradient, the eddy kinetic energy (EKE), the eddy size, and the eddy propagation speed relative to the mean flow. Using kinematic simulations of passive tracer advection, KA14 verified that such a model can accurately describe eddy fluxes in this sector. While we do not attempt to quantitatively fit that model in this study, it is instructive to consider these different factors when comparing the model with the data. Figures 9 and 10 indicate that the dominant length scales and phase speeds in the model and data are very similar. Therefore, we can expect the differences in to be because of the differences in background meridional temperature gradient and eddy kinetic energy.

In Fig. 11b, we plot the time- and zonal-mean meridional SST gradient from the model and the data. We see that the gradients are nearly identical, except in the equatorial region, where the model gradients can be 50% larger. This partly explains the fact that the model is higher in this region. The EKE, defined as , is plotted in Fig. 11c. (Technically only enters directly in the meridional flux, but the eddy velocity statistics are relatively isotropic in the extratropics.) The model and data EKE are quite similar in the extratropics but differ significantly at low latitudes and in the ACC. The model is uniformly more energetic by a factor of 2 in the ACC and at the equator; the regions of EKE mismatch are the same as those of differing eddy heat flux. This mismatch is largely due to the presence of significant ageostrophic currents in the model (F. Bryan 2014, personal communication). Therefore, we can also attribute some of the discrepancy in heat transport to the discrepancy in EKE.

Finally, we plot the surface eddy diffusivity for heat, defined as

Where the gradient vanishes, *D* blows up. To avoid this, we mask *D* wherever the absolute value of the gradient is less than one-tenth of its global-mean value. This occurs in different locations for model and data, making a direct comparison in the equatorial region difficult. Nevertheless, the differences observed in EKE are clearly reflected in *D* as well. Moreover, the general pattern in *D*, low at high latitudes and high at low latitudes, is very consistent with KA14. That study examined the eddy diffusivity of a synthetic passive tracer advected by the AVISO velocity fields. The similarity in the magnitude and spatial structure of *D* between this study and that one suggests that the eddy heat flux is generated through advection of the mean temperature gradient by the geostrophic eddy field.

## 6. Conclusions

This study was largely motivated by KA14, who found that the extratropical eddy flux of a simulated passive tracer in this sector could be parameterized in terms of a single wavenumber and phase speed at each latitude. How can such a monochromatic model be compatible with the broadband variability of the ocean? To resolve this apparent contradiction, we examined the surface eddy heat flux directly using satellite products. We calculated the wavenumber frequency power spectrum at each latitude for surface geostrophic velocity and SST as well as the cross spectrum of these two fields. We integrated these spectra at constant wavenumber, frequency, and phase speed in order to form Figs. 1–7. To quantify the dominant scales and the breadth of their distribution in spectral space, we computed first and second spectral moments. We also performed the same analysis on a high-resolution coupled climate model.

We found that, poleward of 20°, the dominant length scale for both velocity and heat flux is very close to the scale of the nonlinear coherent mesoscale eddies identified by Chelton et al. (2011). For the data, the spectral breadth in wavenumber space, as characterized by the second moment, is between 1 and 2 cycles per 1000 km. For the model, the spectral breadth was greater by about a factor of 2. This probably reflects the fact that the AVISO data are attenuated at scales below 200 km, resulting in narrower spectra. The close correspondence between the length scales of the heat flux and the eddy scales themselves implies a limited role for the atmosphere (e.g., damping of temperature anomalies) in modulating the heat flux. Nevertheless, the way air–sea interaction influences the cross spectra is an intriguing topic for future study.

Similarly, poleward of 20°, the dominant phase speed for both velocity and heat flux is very close to the propagation speed of tracked nonlinear coherent mesoscale eddies. As shown by Tulloch et al. (2009), KA14, and Klocker and Marshall (2014), this speed is approximately equal to *c*_{R} [Eq. (2)], the long Rossby wave phase speed. The spectral breadth in phase speed space is approximately 5 cm s^{−1} for the data. The model spectral breadth is again greater by a factor of 2 or more. Taken together with the wavenumber spectra, these results indicate that the eddy heat flux is generated predominantly by the coherent eddies observed by Chelton et al. (2011).

KA14 reached the same conclusion through a very different approach. They performed kinematic advection–diffusion experiments on a passive tracer driven by AVISO velocity fields and attempted to fit the model of Ferrari and Nikurashin (2010) to the resulting eddy diffusivities. The best agreement was found using the coherent eddy length scale together with the coherent eddy propagation speed—the exact same scales that emerged from the present spectral analysis.

Both these studies are primarily diagnostic, emphasizing the importance of coherent eddy kinetic energy, size, and phase speed in determining eddy fluxes. They do not explain why the eddies have the energy, size, and phase speed that they do. Each of these topics composes an active field of research. The main energy source for mesoscale eddies is the baroclinic instability of the large-scale density field (Gill et al. 1974; von Storch et al. 2012). However, in order to explain the magnitude of the equilibrated eddy kinetic energy, one must also understand how the eddy energy is dissipated, a far less obvious question (Cessi 2008; Arbic et al. 2009; Ferrari and Wunsch 2010; Abernathey et al. 2011). There is growing evidence for the importance of bottom processes in dissipating eddies (Scott et al. 2011; Wright et al. 2012, 2013). Similarly, for the length scales, there is a general belief in the existence of a weak inverse energy cascade in the ocean that transfers energy from a source near the deformation radius to the moderately larger observed scales (Scott and Wang 2005; Tulloch et al. 2011). However, it is far from clear what halts this cascade; proposed mechanisms invoke the importance of *β* (via Rhines scale arguments) and also the role of friction (Held and Larichev 1996; Smith et al. 2002; Lapeyre and Held 2003; Thompson and Young 2006, 2007). The theoretical explanation for the observed eddy propagation speed also remains an active topic of debate (Killworth et al. 1997; Chelton et al. 2007, 2011; Wortham and Wunsch 2014; Klocker and Marshall 2014). Progress on any of these topics would lead directly to progress on the eddy parameterization problem.

One clear shortcoming of the approach used here, which is constrained by the satellite data, is that it does not address the subsurface. Although eddy energy and flux peak near the surface, significant energy and transport exist at depth, particularly in the Southern Ocean. Many methods have recently been proposed to extrapolate satellite data into the interior (Lapeyre and Klein 2006; Isern-Fontanet et al. 2008; Scott and Furnival 2012; Wang et al. 2013; Smith and Vanneste 2013). Future studies could attempt to take advantage of such methods to calculate the cross spectra of subsurface eddy heat fluxes. Novel methods have been applied to estimate subsurface eddy heat transport from Argo data (Dong et al. 2014), but such methods are unlikely to yield results in the spectral domain.

Finally, the inconclusive nature of our findings in the equatorial region is a clear call for further study. The eddy heat flux here is dominated not by mesoscale eddies but by much larger tropical instability waves (Jochum and Murtugudde 2006). The model and data disagree significantly near the equator. This is likely due in part to the difficulty of reconstructing equatorial currents from SSH data alone, but model biases cannot be ruled out either. Furthermore, the vanishingly weak SST gradients result in weak, noisy heat fluxes. KA14 also noted that the Ferrari and Nikurashin (2010) theory did not fare well in the equatorial region. Therefore, a deeper understanding of how tropical instability waves mix and transport heat in this climatically important region should be a top priority.

## Acknowledgments

R. Abernathey gratefully acknowledges support by NASA ROSES Grant NNX14AI46G. Comments by two anonymous reviewers greatly improved the manuscript. We thank Justin Small and Frank Bryan for making the POP model output available through the Earth System Grid. The CESM simulation was conducted with resources provided by the Computational and Information Systems Laboratory of NCAR under the Advanced Scientific Discovery Program.

## REFERENCES

*J. Geophys. Res.,*

**114,**C02024, doi:.

**23,**6277–6291, doi:.

*Geophys. Res. Lett.,*

**34,**L15606, doi:.

*Nat. Commun.,*

**5,**3294, doi:.

*J. Geophys. Res.,*

**105,**19 477–19 498, doi:.

*J. Geophys. Res.,*

**115,**C10048, doi:.

*J. Geophys. Res.,*

**113,**C09005, doi:.

*J. Geophys. Res.,*

**109,**C07002, doi:.

*Geophys. Res. Lett.,*

**39,**L12610, doi:.

*Geophysical Fluid Dynamics.*2nd ed. Springer-Verlag, 710 pp.

*J. Geophys. Res.,*

**116,**C09029, doi:.

*J. Adv. Model. Earth Syst.,*

**6,**1065–1094, doi:.

*J. Geophys. Res.,*

**114,**C02005, doi:.

*Geophys. Res. Lett.,*

**35,**L20601, doi:.

*J. Geophys. Res.,*

**117,**C03049, doi:.

## Footnotes

^{1}

One slightly confusing yet well documented fact to keep in mind is that coherent nonlinear eddies propagate at this phase speed, but larger (apparently linear) Rossby waves propagate somewhat faster (Chelton et al. 2007, 2011).