Energy flux is a fundamental quantity for understanding internal wave generation, propagation, and dissipation. In this paper, the estimation of internal wave energy fluxes 〈u′p′〉 from ocean observations that may be sparse in either time or depth are considered. Sampling must be sufficient in depth to allow for the estimation of the internal wave–induced pressure anomaly p′ using the hydrostatic balance, and sufficient in time to allow for phase averaging. Data limitations that are considered include profile time series with coarse temporal or vertical sampling, profiles missing near-surface or near-bottom information, moorings with sparse vertical sampling, and horizontal surveys with no coherent resampling in time. Methodologies, interpretation, and errors are described. For the specific case of the semidiurnal energy flux radiating from the Hawaiian ridge, errors of ∼10% are typical for estimates from six full-depth profiles spanning 15 h.
Energy flux FE = 〈u′p′〉 = cgE is a fundamental quantity in internal wave energetics to identify energy sources, wave propagation, and energy sinks. Internal wave radiation transports energy from the boundaries into the stratified ocean interior for turbulence and mixing (Munk and Wunsch 1998). Arguably, it is the piece that is missing from 1D boundary layer parameterizations (e.g., Mellor and Yamada 1982; Price et al. 1986; Large et al. 1994; Baumert and Peters 2004; Johnson et al. 1994a, b), representing a potential sink of boundary energy in local budgets.
Until recently, internal wave energy fluxes in ocean observations were estimated by measuring the energy and using the dispersion relation to quantify the group velocity (e.g., Kunze and Sanford 1984; Mied et al. 1986; Kunze et al. 1995; D’Asaro et al. 1995). This requires measuring the vertical, zonal, and meridional wavenumbers, which is only possible in wave fields that are dominated by a single wave and sampled with profile surveys that are conducted rapidly enough to avoid temporal aliasing.
A more powerful and flexible means of estimating the net energy flux is with the velocity–pressure correlation 〈u′p′〉. Here, the principal complications are estimating the internal wave–induced baroclinic pressure perturbation p′ and determining how to average over the wave phase ϕ. For hydrostatic internal waves (ω ≪ N), the pressure anomaly p′ can be obtained by vertically integrating full-depth profiles of the density perturbation ρ′ using the hydrostatic balance. There remains the problem of the integration constant (which can be thought of as the internal wave–induced surface pressure perturbation).
Early work apparently did not recognize the contribution from surface pressure, and so the vertical distribution was incorrect (Garcia Lafuente et al. 1999). However, Ray and Mitchum (1997) and Cummins and Oey (1997) recognized that the vertically integrated baroclinic energy flux is independent of the integration constant because the depth integral of the baroclinic velocity vanishes, ∫0−H u′ dz = 0.
Kunze et al. (2002) took advantage of the fact that, like the baroclinic velocity, the baroclinic pressure perturbation also has a zero depth average ∫0−H p′ dz = 0. This allows one to estimate the integration constant from full-depth profiles of density perturbations. They used their measurements to show that convergence of up-canyon energy flux in the deeper Monterey submarine canyon balanced turbulence dissipation rates. Carter and Gregg (2002) did the same in the shallow end of the canyon. Since then, similar energy-flux computations have been used to quantify internal tide generation at ridges in observations (Althaus et al. 2003; Rudnick et al. 2003; Lee et al. 2005, manuscript submitted to J. Phys. Oceanogr., hereafter LEE05; Nash et al. 2004a, manuscript submitted to J. Phys. Oceanogr., hereafter NA04A; Rainville and Pinkel 2005, manuscript submitted to J. Phys. Oceanogr., hereafter RP) and models (Merrifield et al. 2001; Merrifield and Holloway 2002; Simmons et al. 2004) to examine the loss of a low-mode internal tide to high modes and turbulence over a near-critical continental slope (Nash et al. 2004b), to examine internal wave energy fluxes on a broad continental shelf (MacKinnon and Gregg 2003), and to quantify semidiurnal and inertial energy fluxes globally using historical mooring data (Alford 2003). Rigorous error estimates of FE have not been made in the above studies.
In this paper, methodologies for calculating the internal wave energy flux from ocean observations will be described (sections 2, 3). Uncertainties and biases that are introduced by limited temporal sampling in profile time series are described in section 4 and those with sparse vertical sampling in current meter moorings are described in section 5. Conclusions are presented in section 6. The effect of wave advection of horizontal density and velocity gradients is discussed in appendix A. Our method of generating synthetic space–time series consistent with a Garrett and Munk (1975, henceforth GM) internal wave spectrum is described in appendix B.
2. The internal wave energy equation
The role of the energy flux FE = 〈u′p′〉 = cgE can best be assessed from the internal wave conservation of energy equation
where E = KE + APE is the energy density, KE = (〈u2〉 + 〈υ2〉 + 〈w2〉)/2 is the kinetic energy density, APE = N2〈ξ〉/2 = 〈b′2〉/(2N)2 is the available potential energy density, ∂E/∂t is the accumulation or depletion of energy density, Q+ represents sources such as wind forcing of near-inertial waves (D’Asaro et al. 1995; Alford 2003) or internal tide generation by tide/topography interactions (Bell 1975; Baines 1982; Ray and Mitchum 1997; Althaus et al. 2003), and Q− sinks are associated with loss to turbulent dissipation and the background geostrophic flow.
In a steady-state balance, the energy flux represents transport of energy from sources Q+ to sinks Q−. In the context of the global internal wave field, recent evidence indicates that sources can lie thousands of kilometers away from sinks (Dushaw et al. 1995; Cummins et al. 2001; Nash et al. 2004b; Alford 2003), consistent with altimetric observations of the mode-1 internal tide radiating far northward from the Hawaiian ridge (Ray and Mitchum 1997), and inferences of a near-inertial wave propagating >300 km from its source to its detection point (Alford and Gregg 2001).
a. Flux calculations
To compute the baroclinic energy flux, the internal wave–induced perturbations in pressure p′ and velocity u′ must be inferred from density ρ and velocity u profiles. First, the density anomaly is estimated as
where ρ(z, t) is the instantaneous measured density and (z) is the time mean vertical density profile. Alternatively, ρ′(z, t) may be defined in terms of the vertical displacement of an isopycnal ξ(z, t) relative to its time mean position, so that ρ′(z, t) = (/g)N2ξ(z, t).
The pressure anomaly p′ (z, t) is calculated from the density anomaly using the hydrostatic equation
Although the surface pressure psurf(t) is not measured, it can be inferred from the baroclinicity condition that the depth-averaged pressure perturbation must vanish:
The perturbation velocity is defined as
where u(z, t) is the instantaneous velocity, u(z) is the time mean of that velocity, and uo(t) is determined by requiring baroclinicity
The baroclinicity conditions (4) and (6) require full-depth profiles to determine the barotropic contribution. If data are obtained from a small number of discrete moored instruments, identifying the barotropic signal may, thus, be difficult.
Ocean variability spans a broad range of temporal and spatial scales. In the frequency domain, u and ξ are usually dominated by a few narrow peaks (often in tidal and inertial bands) overlying a broad red continuum (Fig. 1). If long and continuous time series are available, the computation of FE in any frequency band is straightforward. Unfortunately, data are usually sparse, so it is important to determine what sampling and averaging (i.e., 〈. . .〉) is necessary for FE = 〈u′p′〉 to represent a meaningful energy flux for the spectral peak of interest, and to place bounds on its accuracy. Moreover, a number of the quantities in (2), (3), and (5) are difficult to measure from sparse data. First, (z), p(z), and u(z) represent vertical profiles of an undisturbed water column in the absence of internal waves [i.e., u(z) = T −1∫T0 u(z, t) dt, where T is the averaging interval, with T ≫ wave period]. For coarsely sampled data, it may be difficult to differentiate the background state (which may be slowly varying as a result of mesoscale activity) from a wave-perturbed state.
The background state slowly evolves because of mesoscale processes. Typical geostrophic near-surface velocity and pressure fluctuations that are associated with a 50 km × 200 m deep feature are umeso ∼0.2 m s−1 and pmeso ∼600 Pa. While these have the potential to produce instantaneous ∫umesopmeso dz ∼ 1 kW m−1, the mesoscale contributions are at a low frequency and are spectrally distinct from the internal wave band (Fig. 1). Over short time scales, these act to define the mean, but do not alter internal wave perturbation fields. Furthermore, umeso scales with dpmeso/dy, not pmeso, so that extrema in pmesooccur at near-zero umeso and vice versa. As a result, spatially integrating over a smooth geostrophic front produces ∫∫umesopmeso dx dz = 0. Mesoscale processes may also Doppler-shift narrowband peaks into broadband wave fields. In extreme cases, it may be most appropriate to use full spectral methods, such as those of RP, to compute energy fluxes.
In the following analysis, we treat the wave field as a number of narrowband peaks that are superimposed onto a broadband continuum. To estimate FE in a specific frequency band of interest, harmonic analysis or time-domain filtering must be performed to extract the frequency band of interest. Standard time series techniques determine one’s ability to separate different frequency constituents, because the frequency resolution is Δω = 2π/T, where T is the record length.
b. Example time series
A synthetic time–space series of the internal wave field, based on a site of strong semidiurnal internal tide propagation, is used to illustrate a number of aspects of the energy-flux calculation (Fig. 2). This wave field is consistent with the density and cross-ridge velocity observations that are obtained at absolute velocity profiler (AVP) station 14 (in Kauai channel) on 28 October 2000 during the Hawaii Ocean Mixing Experiment (HOME; NA04A; LEE05). The station is located approximately 20 km from a major generation site of semidiurnal internal tides, and is characterized by a strong depth-integrated energy flux (−19 kW m−1), a vertical structure that exhibits both vertically propagating and vertically standing semidiurnal internal waves, and a weaker near-inertial signal that contains a significant high-vertical-wavenumber shear, but little energy flux (<1 kW m−1). Both beams and low vertical modes contribute to the semidiurnal flux. Its characteristics are summarized in Table 1.
Figure 2 displays the vertical and temporal structure of each wave component, and the interplay between these in forming the total energy flux. First, we examine the vertical structure. The cross-ridge velocity u′ and vertical displacement ξ (upper two rows) have a similar vertical wavenumber content, but the amplitudes of u′ scale with N, while ξ scales with 1/N. Hence, u′ is surface intensified, while ξ is bottom intensified. Because p′ represents a vertical integral of ξ, its spectrum is redder by a factor of k2z than that of isopycnal displacement ξ (or velocity u′). It is, thus, dominated by the lowest few modes. As a result, p′ acts as a low-pass filter so that the u′p′ correlation is also dominated by low modes. Hence, only the low-mode component of u′ contributes to u′p′ for typical u′ and ξ (with spectra ∝ k−2z).
A striking feature is that the total instantaneous energy flux 〈u′p′〉 (lower-right panel, shading) is neither equal to the semidiurnal signal (lower-left panel) nor the sum of semidiurnal, near-inertial, and GM fluxes (lower-right panel, solid line), despite the fact that the energy fluxes in the near-inertial and GM fields are very weak. This is because the total instantaneous energy flux is not a linear combination of the constituent fluxes, but is a quadratic quantity containing cross terms, as illustrated below.
The observed oceanic velocity and pressure perturbation may be written as
where M2, f, and GM represent the semidiurnal, near-inertial, and Garrett–Munk variability, respectively, and adv represents the contribution to u′ and p′ from the horizontal advection of low-frequency features by wave motions (see appendix A). The total energy flux is then
While u′M2p′M2 is the leading term, the instantaneous off-diagonal contributions are not negligible, as illustrated by comparing the right and left columns of Fig. 2. For a long, well-sampled time series, these contributions average to zero, so that perfect sampling produces unbiased results. It is the purpose of the following sections to determine the error and bias for the case where sampling is sparse in either the spatial or temporal domains.
Although the following analysis focuses on the semidiurnal wave band (because it is often the most energetic; see Fig. 1), these techniques apply equally well to spectral peaks at other frequencies. And, while the chosen vertical structure and phasing of u′ and ξ is consistent with a site dominated by M2 variance, we assess the error for a wide range of relative amplitudes of semidiurnal, near-inertial, and GM signals, making this analysis applicable to sites where semidiurnal velocity variability is not dominant.
For a narrowband wave field, which is perhaps dominated by near-inertial or semidiurnal waves, as is often the case in the ocean, averaging 〈. . .〉 is optimal over an integral number of wave periods. For a single propagating wave, u′ = u′o cos(kx + mz − ωt − ϕu), phase averaging can be accomplished by averaging in the horizontal x, depth z, or time t. For a vertically standing wave, u′ = u′o cos(kx − ωt − ϕu) cos(mz), phase averaging must be over distance x (or time t) and depth z. For a horizontally standing wave, u′ = u′o cos(mz − ωt − ϕu) cos(kx), phase averaging must be over depth z (or time t) and distance x. For a broadband wave field with no single spectral peak, time series measurements of a sufficient resolution and duration to cover all frequencies are needed, and the energy flux is best determined spectrally.
Because oceanic measurements cannot span an integral number of periods for all waves of interest, the sample mean will contain a contribution from the wave field itself. If the wave field is dominated by a few spectral peaks (Fig. 1), harmonic analysis can efficiently extract the mean and perturbation quantities from relatively coarse-resolution observations, that is,
where the amplitude u′o(z) and phase ϕu and u(z) are determined through a least squares minimization of the residual Ru(z, t). By performing a similar harmonic analysis to p′(z, t), a narrowband estimate of the semidiurnal energy flux may be computed as
where p′o and ϕp represent the semidiurnal amplitudes and phases of the pressure perturbation.
4. Data with coarse temporal sampling
In this section, we illustrate how coarse temporal sampling affects the error in 〈u′p′〉 estimates. A framework is presented to quantify the standard error given a dataset that is completely resolved in depth but coarsely resolved in time, such as that obtained from free-fall profilers [e.g., AVP, expendable current profiler (XCP), high-resolution profiler (HRP)] or the shipboard lowered acoustic Doppler current profiler (LADCP)/CTD. We then apply this framework to a synthetically generated profile time series to illustrate the effects of temporal sampling for several different baroclinic tide regimes.
We wish to compute the time-averaged energy flux in a single-frequency band. In the following, we consider the semidiurnal wave band, denoted as 〈u′M2p′M2〉, although the analysis pertains equally well to any frequency waveband. If the wave field is dominated by semidiurnal variance, then a harmonic analysis to as few as four profiles over a 12.4-h period may extract the signal if phasing of the sampling is fortuitous. Error will be decreased by increasing the number of samples and/or the time series duration. It is the purpose of this section to establish error bounds on such a calculation, and to help guide sampling strategies.
Without rapidly acquired time series or a priori knowledge of the frequency content, it is not possible to determine how much signal has been aliased into the data from nonsemidiurnal constituents. One means of estimating the contamination is through statistics from data-inspired Monte Carlo simulations. We apply the following procedure to our observational datasets.
At each depth, perform a harmonic analysis using as many frequency constituents as the data allow [i.e., (n–1)/2 frequencies or less, where n is the number of samples]. Resolved time series or physical intuition should guide the choice of frequencies. Equivalently, harmonic analysis may be performed on time series of vertical mode amplitudes. In the following, we consider the semidiurnal and near-inertial bands, shown in mooring records to dominate the variability near the Hawaiian ridge (i.e., Levine and Boyd 2004, manuscript submitted to J. Phys. Oceanogr.).
For the purpose of error analysis, assume that the above harmonic analyses have captured all of the narrowband variance. Many realizations of synthetic time series may then be generated by adding random barotropic phases to each frequency component; a randomly phased broadband GM wave field is also added (see appendix B).
Compute error statistics by sampling time–space realizations that are consistent with the observations.
In the following, we use data from the Hawaii Ocean Mixing Experiment to illustrate the above procedure and to determine the error that is associated with typical energy-flux calculations. Time series of full-depth density and velocity were acquired with the AVP (Sanford et al. 1978, 1985; LEE05) at the 3000-m isobath in Kauai Channel during October 2000. This site (station 14; NA04A; LEE05) exhibits strong internal tides (−19 kW m−1 depth integrated) and significant near-inertial activity.
At HOME AVP station 14, six full-depth profiles over a ∼15 h period were acquired. Data were processed as outlined in the previous section and u′ and ξ were computed for each profile. Harmonic analyses were performed to determine the amplitude and phase at semidiurnal and inertial frequencies so that the measured velocity is u′ = u′M2 + u′f + Ru and vertical displacements are ξ′ = ξ′M2 + ξ′f + Rξ, where Ru and Rξ represent least squares–minimized residuals. For simplicity, we consider only the cross-ridge component of velocity and energy flux in the following analysis, with ^ defined as 37°E of north; this is roughly aligned with the depth-averaged semidiurnal energy flux at this site.
We employed Monte Carlo methods to test the sensitivity of the energy-flux calculation to various sources of contamination. Synthetic time series were generated with different strengths of each semidiurnal, GM, near-inertial, and instrument noise component. Time series of velocity and vertical displacement were generated as
where the an are nondimensional amplitudes that represent the strengths of the internal wave field components relative to the values in Table 1. Each an represents a wave amplitude, so that the energy or energy flux that is associated with each contribution scales with a2n and not with an (i.e., FE ∝ a2M2). A random instrument noise of u′ins = 0.02 m s−1 and ξ′ins = 0.5 m was arbitrarily prescribed and added to the wave fields with insignificant effect on the flux estimates.
First, benchmark synthetic time series of utest and ξtest were generated using aM2 = af = aGM = ains = 1, and estimates of 〈u′p′〉 that were computed using the following three sampling strategies: minimal (4 profiles in 12 h), efficient (6 profiles over 15 h), and well sampled (12 profiles over 25 h). In each case, 200 energy-flux estimates were computed from harmonic analyses that extract the semidiurnal signal (12); the near-inertial signal was also extracted for n > 5. Vertical profiles of the energy flux, standard error, and bias, computed through M2 harmonic analyses for these three cases, are shown in Fig. 3. A number of observations are immediately evident. The first is that even the minimal strategy qualitatively captures the vertical structure of the energy flux. Second, the standard deviation scales with the magnitude of the energy flux; that is, the error is smallest at middepths where the energy fluxes are the weakest, with 20% error being typical for a single depth estimate. Next, the sampling does not systematically bias the estimates; any bias is a factor of 10 less than the scatter. Finally, while increasing the number of samples and duration does reduce error, the gain in increasing the sampling from six in 15 h to twelve in 25 h is only marginal. This highlights the effectiveness of the harmonic analysis to extract the desired signal (semidiurnal, in this case).
To determine the sources of contamination, a sensitivity study was performed in which the strengths of the semidiurnal, GM, and near-inertial components were independently varied, while holding the strength of each remaining internal wave field constant. That is, 1) aM2 was varied from 0 to 2, while af = aGM = 1; 2) aGM was varied from 0 to 2, while af = aM2 = 1; and 3) af was varied from 0 to 2, while aM2 = aGM = 1. For each combination of aM2, af, and aGM, 200 realizations of u′test and ξtest were generated with (13) and (14) using random phasings of each wave component. These were then sampled according to each of our sampling strategies that are outlined in the following subsections and statistics of the depth-integrated energy flux ∫FE dz computed to diagnose the sources and magnitude of error and bias. Estimates were computed using both simple averages (10) and harmonic analyses (12).
The cross-ridge energy flux was computed from either the simple average of the perturbation quantities as FE (10) or from a harmonic analysis as FM2E (12). The standard deviation σ and bias μ of ∫FE dz or ∫FM2E dz was then computed from 200 independent realizations of the wave field, and divided by the prescribed ∫FM2E dz to yield the fractional error and bias. Using this definition, energy fluxes are considered to be biased if a near-inertial signal aliases into our estimates. Results are shown in Figs. 4–6 and Tables 2–4. The lines in each figure represent the sensitivity of the fractional error or bias to the amplitude of semidiurnal signal (Figs. 4a, b; 5a, b; and 6a,b) or contamination [GM (Figs. 4c, d; 5c, d; and 6c,d) or near inertial (Figs. 4e, f; 5e, f; and 6e,f)] for a particular sampling strategy.
In the following, we investigate errors and biases that are associated with two types of oceanic observations. The first is sampling over a short duration (12–24 h) using a uniform ∼3 h time step (section 4a). This is common of intensive process studies (i.e., Kunze et al. 2002; Althaus et al. 2003; Nash et al. 2004b; NA04A) that are designed to capture a snapshot of the spatial structure of an internal wave field. The second is data collected at irregular sampling intervals (section 4b). This represents data of opportunity, such as that obtained during the World Ocean Circulation Experiment (WOCE) transects where sampling is irregular in time or space with few repeated stations.
a. Regularly sampled time series
We first consider the case of a time series that is sampled n times at equal intervals over the period T. Synthetic time series were sampled using one of the three sampling strategies (minimal, efficient, or well sampled), from which statistics of the energy flux were calculated as simple averages [Eq. (10), Fig. 4, and Table 2] or from harmonic fits [Eq. (12), Fig. 5, and Table 3].
For the benchmark case of aM2 = 1, the error in ∫FE dz is 9%–16% and bias is less than 5%, with lower percentages corresponding to increased n. Using simple averages, the GM and near-inertial wave fields contribute to the error with a similar magnitude. Harmonic analysis is able to reject near-inertial contamination almost completely for n ≥ 6, because such sampling permits both M2 and near-inertial harmonic analyses to be performed at this latitude. In contrast, a harmonic analysis can neither extract near-inertial variance for n = 4, nor reject GM contamination, which is broadband and aliases into the M2 band from many frequencies. The wave field, hence, contributes most to the error.
Bias is substantially reduced through harmonic analysis for two reasons. The first is because energy-flux estimates from simple averages contain contributions from all frequencies, including those that are near inertial. In the example presented here, the near-inertial contribution has the same sign as the M2 energy flux, and, hence, the estimates are biased high. The near-inertial energy flux is 1.9% of the M2 for the benchmark case (Table 1); whereas our bias estimates for the well-resolved n = 12 case is 1.3%.
Bias also results from error in the sample mean u and p, which is subtracted from the observations to determine perturbations. Any such error reduces the variance of the perturbations.1 This produces a bias toward lower energy fluxes, which we find for all poorly sampled estimates (n < 12). Only for the case of aM2 = 0.4 are the energy fluxes biased consistently high. This bias results from the near-inertial contribution, which carries 12% of the energy flux for aM2 = 0.4 (i.e., 〈u′M2p′M2〉 = −3 kW m−1, whereas 〈u′f p′f〉 = −0.36 kW m−1). These trends are evident in Fig. 4b and are reduced through harmonic analysis (Fig. 5b).
We conclude that harmonic analyses are effective in extracting the M2 signal from one contaminated with GM and near-inertial waves. For large energy fluxes, such as those in HOME, a harmonic analysis to as few as n = 4 samples provides a reliable estimate of ∫FE dz. For weaker fluxes (i.e., aM2 = 0.4, ∫FE dz = 3 kW m−1), n ≥ 6 is required to reduce the error to 25% or less. We emphasize that these error bounds are based on depth integrals over a 3100-m water column, with GM and near-inertial wave fields that are specific to the Hawaiian ridge. We also note that harmonic analyses can only be performed on profiles that are collected at the same location, because temporal phase is distorted if profile locations are displaced as an appreciable fraction (i.e., >10%, ∼15 km) of the dominant horizontal wavelength.
Aspects of this analysis may be specific to the HOME site. First, the separation of near-inertial and semidiurnal frequencies depends on latitude. Second, the spectral bandwidth of these peaks changes with location (i.e., distance from generation site, mesoscale intensity, degree of Doppler shifting, etc.). Both factors alter the effectiveness of harmonic analyses to extract a signal of interest. At high latitudes, for example, it may be impossible to distinguish semidiurnal from near-inertial variability if the time series duration T < 2π/(ωM2 − f ).
b. Incoherently sampled wave field
We now consider the case of a space or time series that is sampled n times at irregular intervals over some long spatial or temporal period. While not an ideal sampling scheme, this section is included to evaluate the usefulness of data of opportunity to compute internal wave fluxes. For example, spatial transects with periodic LADCP/CTD data may be considered as a random collection of independent profiles. Such sampling does not capture all phase of the M2 signal equally, and alters the sample estimates of both the mean and perturbation fields. Here, we consider n = 4, 6, 8, 12, or 24 profiles acquired at random times over a T = 30 day duration (with a minimum ΔT = 3 h between the profiles). These may represent profiles that are acquired at a single location over a long duration, or at a number of different locations over which the internal wave field may be considered homogeneous. The length of the total sampling period is unimportant, provided that each wave field is stationary and randomly phased (as we assume). Statistics of the energy flux are shown in Fig. 6 and Table 4 for simple average estimates.
A comparison of Fig. 6 (Table 4) with Fig. 4 (Table 2) indicates the error and bias for n profiles that are acquired randomly are 2–3 times larger than those obtained at regular intervals. Also note that these estimates assume stationarity of the wave field over the sampling period. Given the intermittency of the internal tide (Wunsch 1975) on O(5 day) time scales, T should be chosen as short as possible, while still permitting an appropriate n.
5. Data with imperfect vertical sampling
Vertical deficiencies of typical ocean data fall into two general categories: sparse, discrete measurements at a number of depths (e.g., moorings), and highly resolved data over a portion of the water column (e.g., profiles). We will consider these two cases separately, assuming that the data are perfectly well sampled in time. (For discrete-depth moorings, this is generally an excellent assumption.) In this case, contamination by motions of other frequencies (e.g., GM background and near-inertial signals) is less of an issue, because filtering or harmonic analysis can remove these signals. Otherwise, the principles in the previous section apply.
As discussed earlier, the determination of baroclinic pressure anomaly requires full-depth, continuous data. Because neither type of moored dataset satisfies this requirement, we rely on normal modes to generate the full-depth profiles. Success will depend on the density of the sampling and the redness of the ocean vertical wavenumber spectrum at the frequency of interest.
The modal amplitudes are determined at each measurement time by solving a weighted least squares problem,
where ûj and ξ̂j are the jth modal amplitudes of velocity and displacement at each time, and Zuj(z), Zξj(z) are the normal mode structure functions. For velocity, the zeroth mode is rigorously the barotropic mode. Each set of i equations (one for each measurement depth) can be put in matrix form, and solved using standard overdetermined inverse methods (Dushaw et al. 1995). Formally, Mmeas instruments can resolve Mmeas modes. In practice, the specific mooring geometry and the shape of the spectrum can further limit this value.
In analogy with the temporal case, contamination arises when unresolved variance is projected onto resolved modes. The severity of the problem depends on the redness of the ocean spectrum, the geometry of the vertical sampling, and, in some cases, the number of modes being solved for. The aim of this section is to provide guidance in both designing vertical arrays and interpreting flux measurements of opportunity with imperfect sampling. Therefore, in parallel with the previous section, we conduct Monte Carlo simulations of energy-flux calculations. For each combination of spectral redness and mooring geometry, a synthetic signal was generated by superposing 30 modes with random phases for velocity:
The parameter q governs the partition of energy among the modes. It is not known a priori, but will vary from place to place and in time. Guided by HOME data, we choose values of q spanning −1, the observed HOME value.
Several schemes for randomizing phase were employed, none of which materially affected our results. Completely randomizing the phase of all the modes would cause the “true” flux in each realization to vary wildly. To avoid this, we chose to fix the phase of the second mode relative to the first, so that the integrated flux in each realization was near −19 kW m−1, as in the previous section.
For each pairing of vertical geometry and spectral redness, 100 realizations of synthetic data were generated, and the flux was computed over one M2 cycle for each using the discrete and a perfectly resolved array. The fractional bias and standard error were computed for each as the mean and standard deviation of the flux difference between the discretely sampled and perfect arrays.
a. Sparse vertical sampling
We first consider fluxes determined from current and density measurements at discrete depths. To collapse the infinite number of possible mooring geometries into a spannable parameter space, we assume that the instruments are “ideally situated,” which we define to be evenly spaced in a WKB-stretched depth coordinate z′, defined as
In this coordinate, the modes are nearly sinusoidal and the instruments, at z′ins = (H/Mmeas)[(1, 2, 3, . . . , Mmeas) − 1/2], are located to best resolve them. Because it is generally desirable to resolve the dynamical modes, most moorings are designed with this in mind.
The fractional error and bias resulting from various numbers of ideally spaced instruments are shown in Fig. 7. For each number of instruments Mmeas, Mmeas−2 modes were solved for (as will be seen in the next section, solving for fewer modes results in a moderate decrease in precision, but can result in greater stability). Each curve represents a different spectral redness, where the HOME case is in the middle. With 40 instruments, error and bias are zero for all input spectra because more modes are resolved than are present in the ocean spectrum (30). As the number decreases, bias remains near zero, but the standard error increases, reaching 0.75 for four instruments in a spectral slope of −1. Thus, the number of realizations that are required for stable averages grows for sparse moorings. The lack of bias indicates that long-term averages will eventually converge on the correct value.
These errors and biases are for an estimate of the flux in a single M2 cycle. Because the fractional errors for a small instrument number approach unity, it desirable to know how much averaging is required to yield more precise estimates. Obviously, the error decreases as the square root of the number of independent observations. Unfortunately, the decorrelation time scale of the high-wavenumber contaminants is unknown. For example, if a phase-locked internal tide beam lying outside or on the edge of the instruments’ range is the source of the contamination, no amount of averaging will help. At the opposite extreme, if the phase of the high-wavenumber motions is random from one cycle to the next, then the error in each single-cycle flux estimate will be uncorrelated. In this case, the fractional error in an average over, say 4 days, will be ≈3 times lower than the values cited here. In practice, this is likely to be about as good as can be done, because flux itself is variable on these time scales (Wunsch 1975; Alford 2003).
b. Profiles with gaps at the top or bottom
Because energy flux WKB-scales as buoyancy frequency N, it is typically extremely surface intensified. Consequently, the calculation can become unstable without measurements close to the surface. Here we consider the effects of a gap at the top of an otherwise well-resolved profile. The stability of the calculation depends, in this case, on the number of modes to be solved for. Therefore, we present results for two-, four- and eight-mode solutions (Fig. 8). Considering the former first (top panels), errors grow as the top instrument is deepened, but, again, bias remains nearly zero except for the “bluest” case. By 250 m, the fractional error is unity for the HOME case.
Increasing the number of modes to be solved for (middle, lower panels) decreases the error for small top gaps, as it does for well-instrumented discrete arrays (not shown). This is because the higher modes are better resolved and cannot project onto the flux-carrying low modes, introducing errors. However, solving for more modes carries the penalty of decreasing stability, as is evident in the much larger errors and strong biases in the eight-mode case, when solving for more modes. That is, higher-mode solutions are more precise but more sensitive to gaps.
Much larger gaps at the bottom are tolerable (Fig. 9). This is to be anticipated, given that the upper water column is weighted more strongly in WKB-stretched coordinates. The results are generally similar to the top-gap case, but much larger distances from the bottom are spanned. In the two-mode case, fractional error is nearly constant at 0.3–0.4, even when only 4300 − 3000 = 1300 m of the water column is spanned. This is a result of the surface intensification of the flux profiles and the dominance of the lowest modes in carrying the flux. Higher-mode solutions, as before, yield greater precision for the full-column coverage but have earlier instability and greater errors as the gap is increased.
A framework for assessing the error and bias of baroclinic energy-flux estimates has been presented. Our method employs data-based Monte Carlo simulations to assess the magnitude and parameter dependence of flux estimates made from (a) temporally or (b) vertically imperfect data.
In the former case, we conducted simulations by varying M2, GM, and near-inertial energy about realistic values. We find that a 10% error is typical for estimates based on n = 6 profiles spanning 15 h, such as those collected with AVP during the Hawaii Ocean Mixing Experiment (Rudnick et al. 2003; NA04A; LEE05).
In addition, we conclude that
unbiased semidiurnal energy-flux estimates can be computed from n = 4 profiles over 12 h;
the vertical structure of the energy-flux profile is qualitatively captured by as few as n = 4 profiles;
if spectral peaks are distinct, harmonic analyses applied to regularly sampled time series are effective at rejecting contamination from other narrowband frequencies, but only slightly reduce the error associated with broadband GM contamination;
wave advection of strong meso- and submesoscale fronts only weakly contaminates the depth integrated energy flux (appendix A); and
harmonic analyses should not be used to analyze sparsely sampled time series with nonuniform sample spacing, because GM and near-inertial contamination can alias unpredictably into the wave of interest.
We also considered the case of temporally well-sampled but vertically gappy data, as in the case of discretely instrumented moorings or partial water column–moored profilers. Here, temporal contamination can be eliminated via filtering or harmonic analysis, but normal-mode fits must be employed via least squares methods to determine the necessary depth integrals. We find the error to be highly sensitive to the slope of the internal tide spectrum. Bluer ocean spectra and sparser measurements yield poorer fits and larger fractional errors in the energy flux.
Our specific results for the vertically deficient case include the following.
For a “typical” mooring with six ideally placed instruments and the HOME spectrum, the fractional error is 40%.
Estimates are unbiased, except in the case of very large gaps at the top or bottom.
Estimates are sensitive to data near the surface. Errors increase rapidly as the depth of the top measurement increases.
As a result of the WKB weighting of energy flux, larger gaps can be tolerated at the bottom than at the top.
In all cases, solving for more modes reduces the error, but can affect the stability of the solutions if modes with a wavelength comparable to the WKB-stretched gap are solved for.
We conclude, with a final cautionary note, that even perfectly resolved flux estimates have subtleties that can require care in their interpretation. For example, NA04A and Nash et al. (2004b) considered a horizontal standing pattern resulting from two internal tides propagating in opposite directions. In such cases, the flux in the transverse direction displays a spatial periodicity, requiring spatial as well as temporal averaging. Naive consideration of its unaveraged divergence would lead to spurious conclusions.
The authors thank Tom Sanford, Craig Lee, and the technical support of the Applied Physics Lab (John Dunlap, Art Bartlett, and Bob Drever) for making these measurements possible. J. N. and E. K. were supported by NSF Grant OCE-98-19537 and ONR Grant N00014-97-10087. M. A.’s contribution to this work was supported by the Office of Naval Research Young Investigator award, Grant N00014-02-10526.
Wave Advection of Horizontal Gradients
Horizontal density and velocity gradients that are associated with meso- and submesoscale variability are advected by wave motions and contribute to u′ and p′ as measured at fixed (xo, yo). This appendix assesses the importance of u′adυ and p′adυ in contaminating FE estimates [i.e, through Eq. (10)].
Consider a wave field with velocity u′wave ≡ u′M2 + u′f + u′GM. At time t′, the fluid at position xo was located at xadv = xo − ∫t′to u′wave dt at time t = to. The advective contribution to ρ at t′ is, therefore, ρadv = ρo(xadv) − ρo(xo), and will produce a pressure perturbation padv through the hydrostatic balance. Similarly, the advective contributions to u′ and υ′ at t′ are uadv = uo(xadv) − uo(xo) and υadv = υo(xadv) − υo(xo).
For the purpose of this analysis, we assume that uo, υo, and ρo correspond to a surface-intensified front in geostrophic balance with the level of no motion at z = −500 m. Our benchmark case (afront = 1) corresponds to a series of 10-km-wide fronts with 20-m peak-to-peak vertical displacements that induce ±0.9 m s−1 surface currents and a ±0.43 m s−1 maximum average velocity in the upper 500 m. Frontal vertical displacements vary periodically in x as ξo = afrontξmax cos(2πx/λ)Z(z), where ξmax = 10 m, λ = 20 km, and Z(z) is piecewise linear in the vertical (increasing from 0 to 1 between 0 and −200 m, returning to 0 at −500 m). For afront = 1, the mean surface velocity gradient 〈|duo/dy|〉 is 0.18 m s−1 km−1. These gradients are extremely sharp and their vertical structure projects strongly into mode 1. Hence, we consider this to represent a worst-case scenario for advective contamination.
A front varying only with x is advected by the î component of wave velocity u′wave. This produces no uadv because uo = 0 in the front (only υo ≠ 0). Conversely, a front varying only with y is advected by υ′wave and alters u′ through uadv; because υadv = 0, υ′ is unaltered. It is, thus, necessary to consider both velocity components of the wave field in addition to their orientation with respect to the front.
We employ the Monte Carlo techniques of section 4 to assess the cumulative effect of all terms in (10) involving u′adv and p′adυ (as well as those in the υ′p′ equation involving υadv). Simulations were performed in which fronts were oriented in either î and ĵ directions. For each frontal orientation, the following two sensitivity studies were conducted: 1) frontal amplitudes were varied from 0 < afront < 2, while holding awave = 1 fixed; and 2) semidiurnal wave amplitudes were varied from 0 < awave < 2, while holding afront = 1 fixed. Here awave represents the total wave amplitude, so that aM2 = af = aGM = awave were covaried together. For each case, the error and bias was computed from 200 independent realizations with randomly phased wave fields that were superimposed onto randomly located fronts. In all cases, we assume perfect temporal and spatial sampling, and compute estimates using harmonic analysis. The results of these simulations are summarized in Fig. A1.
Despite the strength of the prescribed front, the advective contamination into FE is weak compared to the sampling error presented in section 4. Estimates are unbiased and contamination is restricted to the component of energy flux parallel to the frontal velocity. That is, a front with uo = 0 does not alter the î component of the energy flux. This is because the advection of ρo does not contaminate FE, because both ρadv and the induced p′adυ are in quadrature with u′M2. Hence, only the advection of frontal velocities alters the energy flux so that the error increases linearly with afront. For fixed afront, absolute error increases roughly with the square of wave strength, because both velocity contamination u′adυ and the pressure perturbation p′wave scale linearly with awave, and the induced error is dominated by the correlation 〈u′adυp′wave〉. The fractional error, however, decreases slightly with wave strength because both FE and the advected contamination scale approximately with a2wave.
In contrast to the depth-integrated estimates, the energy flux at a given depth within the upper 500 m has a standard error σ ∼10 W m−2, which is comparable to that associated with sampling errors (Fig. 3). Hence, advection can strongly contaminate the flux profile, while only weakly altering its depth integral. This is the direct result of the weak role of p′adυ in advective contamination, and the much higher wavenumber content of u′.
In summary, the advection of a surface-intensified front with 1.8 m s−1 changes in u over 10 km contaminates the depth-integrated energy flux by less than 20%. We conclude that frontal advection is not likely to strongly contaminate ∫FE dz. While the profile FE(z) may be contaminated by this extremely strong front, more typical frontal strengths will induce errors much smaller error in FE(z) than that associated with discrete sampling.
Generating GM Depth/Time Series
We wished to “contaminate” the tidal signals in the Monte Carlo simulations with time–depth series of velocity and displacement, uGM(z, t), ξGM(z, t), that have the Garrett–Munk spectrum ΦuGM(kz, ω), ΦξGM(kz, ω). Because the GM spectrum is separable by assumption, ΦGM(kz, ω) = ΦωGM(ω)ΦkzGM(kz) for both u and ξ. Essentially, we create Fourier amplitudes with this spectrum, randomize their phases, and inverse transform. Here we outline the detailed procedure for generating uGM(z, t); the procedure for ξGM(z, t) is the same.
For each realization, we generated a matrix of amplitudes ũGM(kz, ω) that was consistent with the internal wave wavenumber and frequency spectra of Garrett and Munk (1975), as modified by Cairns and Williams (1976). That is, ũũ* = ΦuGM(kz, ω). This was done for a set of discrete frequencies and wavenumbers, 0 < ω < ωN, 0 < kz < kzN, where the subscript N indicates the Nyquist frequency or wavenumber. We then randomized the phase of each amplitude by multiplying each by eiΨ, where 0 ≤ Ψ ≤ 2 π is a uniformly distributed random variable.
The amplitudes were then mirrored across zero frequency such that ũ(−ω) = ũ*(ω), to ensure real series. (There is no enforced phase relation between kz and −kz, allowing both vertically standing and propagating contributions.) At this point, the dimension of the amplitude matrix was the same as that of the desired time/depth series.
The amplitudes were then inverse 2D Fourier transformed, producing a depth–time series uWKB(z′, t) with the correct spectrum. However, the wavelength and amplitude are stationary in depth. That is, they mimic “true” ocean field after WKB stretching and scaling to remove refractive effects. The final “ocean”-contaminant depth–time series were produced by WKB “unstretching” using (5a), and “unscaling” using uGM(z) = uWKBGM[N(z)/No]1/2.
Corresponding author address: Jonathan D. Nash, 104 COAS Admin. Bldg., Oregon State University, Corvallis, OR 97331. Email: firstname.lastname@example.org