The effects of upward buoyancy on the accuracy with which Lagrangian floats can measure the Eulerian mean variance 〈ww〉E and skewness SwE of vertical fluid velocity w in the wind-driven upper-ocean boundary layer is investigated using both simulated floats in large-eddy simulation (LES) models and two float datasets. Nearly neutrally buoyant floats are repeatedly advected by the turbulent velocities across the boundary layer. Their vertical position Z is determined from pressure measurements; their W variance 〈WW〉F and skewness SWE are determined from the time series of float W = dZ/dt. If the float buoyancy is small, then the simulated floats measure the Eulerian velocity accurately; that is, δW2 = 〈WW〉F − 〈ww〉E and δSW = SWF − SwE are small compared to 〈ww〉E and SwE respectively. If the floats are buoyant, and thus have an upward vertical velocity Wbias relative to the water, then δW2 and δSW can become significant. Buoyancy causes the floats to oversample both shallow depths and strong vertical velocities, leading to positive values of δW2 and negative values of δSW. The skewness SZ′F of depth measures the degree to which shallow depths are oversampled; it is shown to be a good predictor of Wbias/〈WW〉F1/2, δW2/〈WWF〉, and δSW/SWF over a wide range of float buoyancies and boundary layer forcings. Float data collected during two deployments confirm these results, but averaging times of several float days are typically required to obtain stable statistics. Significant differences in the magnitude of the effect may occur between datasets from different ocean forcing regimes and float designs. Other measures of float buoyancy are also useful predictors. These results can be used to correct existing float measurements of 〈ww〉E for the effects of buoyancy and also can be used as a means to operationally assess and control float buoyancy.
The upper-ocean boundary layer is a key conduit between the ocean and atmosphere. Unlike the ocean interior, the boundary layer is commonly turbulent with its associated eddy fluxes having a first-order role in upper-ocean physics, chemistry, and biology. A detailed understanding of this region has been inhibited by the lack of comprehensive measurements of this turbulence—its intensity, fluxes, scales, structures, and dynamics. Some measurements have been limited to the microscale dissipation rates of energy and scalar variance (e.g., Terray et al. 1996). Measurements of larger-scale energy and fluxes are difficult, but not impossible (Gerbi et al. 2008), because the turbulent velocities are typically much smaller than the velocities resulting from surface waves.
One promising technique is the use of high-drag, neutrally buoyant Lagrangian floats (D’Asaro et al. 1996; D’Asaro 2003a). These meter-sized devices are designed to follow the average three-dimensional motion of the water immediately surrounding them. This is achieved by the float matching the ambient density and having a large vertical drag, so that any residual buoyancy will result only in small deviations from the vertical velocity of the fluid. These floats measure their depth from the pressure at a central location on the float, and determine vertical velocity from its rate of change. Because there are no variations in pressure along Lagrangian trajectories because of linear surface waves, and oceanic surface waves are well approximated as linear, the velocities and accelerations measured by Lagrangian floats are nearly uncontaminated by surface wave velocities. Float measurements of vertical velocity have the magnitude expected of turbulent velocities (Harcourt et al. 2002; D’Asaro 2001; Steffen and D’Asaro 2004), correlate well with atmospheric forcing, and, when combined with scalar measurements on the float, yield useful measurements of scalar covariance fluxes (D’Asaro 2003b, 2004; D’Asaro and McNeil 2006).
Ideally, the float’s velocity should equal the local water velocity at the pressure sensor. A real float, however, is significantly larger than the smallest scales of turbulence; its motion is therefore more representative of a local volume average over fluid parcel trajectories. Lien et al. (1998) investigate how finite float size affects the ability of a float to measure turbulence statistics. In turbulent flows there is a strong theoretical basis (inertial subrange theory) for understanding the form of the Lagrangian frequency spectra of velocity at high frequency. Lien et al. (1998) derive the spectral forms for finite size floats. Measurements made in a variety of energetic turbulent flows (Lien et al. 1998; D’Asaro and Lien 2000) show remarkably similar spectral shapes that accurately fit the Lien et al. (1998) form. Spectra are isotropic within the inertial subrange, which is also in agreement with theory (D’Asaro and Lien 2000). Lien and D’Asaro (2006) show that the values of the turbulent kinetic energy dissipation rate estimated from these spectra agree well with those estimated from smaller-scale inertial subrange wavenumber spectra. These studies provide considerable confidence that these floats are indeed approximately Lagrangian, the effects of the finite float size can be modeled, and the estimates of vertical velocity statistics from the float data are meaningful.
Two additional mechanisms may cause differences between the turbulence statistics computed from Lagrangian floats and from an Eulerian average. First, exactly Lagrangian floats may produce biased statistics if they are not distributed uniformly. For example, floats deployed uniformly within an entraining ocean mixed layer, but not below, will undersample recently entrained fluid. This results in biased statistics, particularly near the mixed layer base. This bias is apparent in Lagrangian floats programmed to stay within the mixed layer. Second, and the main focus here—floats whose buoyancy does not match that of the ambient fluid—will rise relative to their surroundings, resulting in biased statistics, not simply because of the direct effect of this motion, but primarily because it causes the distribution of floats to be inhomogeneous and correlated with local fluid properties (D’Asaro et al. 2002; Harcourt et al. 2002; Noh et al. 2006).
We address the ability of the Lagrangian floats to measure quantitatively the variance and skewness of vertical velocity W in the upper-ocean boundary layer in the case where the float’s buoyancy does not match that of the surrounding water exactly. D’Asaro (2003a) and Tseng and D’Asaro (2004) report remarkably high correlations between the variance of W, or the vertical kinetic energy (VKE) measured by Lagrangian floats, and the wind speed or wind stress in wind-forced boundary layers. This conflicts with the recent analysis of numerical results (Harcourt and D’Asaro 2008, hereafter HD08) demonstrating an additional strong dependence of W on wave properties. Similarly, Sullivan et al. (2004, 2007) suggest an important role for wave breaking in the boundary layer dynamics associated with large values of vertical velocity skewness. We therefore focus on the mixed layer turbulence statistics of W to identify and profit from their sensitivities to the physical mechanisms responsible for upper-ocean mixing. Using such statistics to discriminate between competing theories of mixed-layer dynamics requires empirical data where both random errors and systematic biases are well understood. This paper addresses the systematic biases in Lagrangian float measurements of the turbulent moment’s vertical velocity in mixed layers resulting from imperfect ballasting.
Our approach is to embed simulated imperfect floats in a turbulence-resolving large-eddy simulation (LES) of the wind- and wave-forced upper-ocean boundary layer, as described in section 2. The characteristics of the wind, waves, and floats are varied over realistic ranges and the vertical velocity statistics measured by the floats are compared to the accurate Eulerian averages available from the LES model in section 3. The scaling of the biases in float-measured “bulk,” or layer-averaged w statistics with the float buoyancy is described in section 4, resulting in consistency tests, statistical estimates of the float buoyancy and biases, and algorithms to correct the biases. The model fidelity is tested in section 5 by comparison with float data from two different deployments. The results and implications are summarized and discussed in section 6.
2. LES model and simulation cases
a. LES model
HD08 report steady-state LES solutions for vertical turbulent kinetic energy (TKE) in mixed layers forced by wind and waves. As in previous LES studies of Langmuir turbulence (Skyllingstad and Denbo 1995; McWilliams et al. 1997; Min and Noh 2004; Noh et al. 2004), the interaction of waves and turbulence is modeled by the Craik–Leibovich vortex force (Craik and Leibovich 1976). Horizontally uniform surface stress τ0 and Stokes drift profiles us(z) are specified from the 10-m wind speed U10 and the surface wave age Cp/U10, where Cp is the spectral peak phase speed, using empirical surface wave spectra and an associated wave age–dependent neutral drag coefficient CD. Wave breaking effects are not otherwise included. Mixed layer depths HML vary from 30 to 120 m, with 0.6 ≤ Cp/U10 ≤ 1.2 and 8 m s−1 < U10 < 70 m s−1, thereby addressing most possible oceanic conditions where TKE production is dominated by wind and wave forcing. HD08 obtain a robust scaling for the average mixed layer VKE in terms of us(z), HML, and friction velocity u* = . Further details of the LES model and the scaling of Langmuir turbulence can be found in HD08. Here, we focus on the response of Lagrangian floats to the turbulence resolved in a subset of the HD08 LES cases.
b. LES model cases
Four ocean boundary layer simulations from HD08 were used to span a range of wind speeds U10 from 10 to 54.4 m s−1 and varying mixed layer depth and wave age Cp/U10. Each of these four cases had vertical resolution Δz = 1.42 m and horizontal resolution Δx = Δy = 1.5 m. In HD08, LES cases are labeled “Σi(U10, Cp/U10),” indicating the wind speed U10 and wave age Cp/U10, and their membership in one of five model case subsets (Σi; where i = 1, 2, 3a, 3b, 4). Their forcing and dimensional scales are detailed in HD08 by five correspondingly numbered tables of U10 versus Cp/U10.
1) Case LW: “Low wind”
2) Case MW: “Moderate wind”
3) Case H1: “High wind 1”
The first high-wind case (case H1) is referred to as case “Σ3b(U10 = 54.4 m s−1)” in HD08. It has dt = ⅛ s, a deeper mean mixed layer of 87.8 m, U10 = 54.4 m s−1, and very young seas at Cp/U10 = 0.6. For further forcing information, see Table 3b in HD08. This simulation uses a “saturated” high-wind surface drag coefficient CD = 2.30 × 10−3.
4) Case H2: “High wind 2”
The second high-wind case (H2) is referred to as case “Σ3a(U10 = 54.4 m s−1)” in HD08. It has dt = ⅛ s, a deeper mean mixed layer of 90.9 m, U10 = 54.4 m s−1, and very young seas at Cp/U10 = 0.6. For further forcing information, see Table 3a in HD08. This simulation uses a higher drag coefficient CD = 3.98 × 10−3, based on a formulation wherein CD continues to increase with U10 at high winds.
c. Simulated float types
LES fields represent grid volume averages of the underlying velocities or thermodynamic scalars. Although each model float is considered to be located at a point, the local properties that are spline interpolated to this point correspond to mean properties in that neighborhood. Ideally, real floats are exactly ballasted at all times, with buoyancy changing continuously to match the density of water at its ambient temperature, salinity, and pressure. Each LES simulation included model floats with behaviors ranging from exactly ballasted “ideal” floats to Lagrangian floats that are highly buoyant. The float’s horizontal velocity equals that of the water. The difference between the float’s vertical acceleration and the acceleration of the surrounding water is given by the sum of buoyancy and drag forces on the float divided by the sum of the float’s mass and the virtual mass of the surrounding water [see Harcourt et al. 2002, their Eq. (1)]. The vertical buoyancy force is computed using the assumed float buoyancy and a volume of Vfloat = 0.05 m3. The vertical drag force is computed assuming a horizontal area of 0.5 m2, a drag coefficient Cd = 1, and an added mass from entrained fluid of 9 Vfloat. While accurate drag and inertial characteristics of Lagrangian floats in turbulent boundary layers are not well known, the float model here is designed to have similar dimensional scales and generally exhibits behavior consistent with the real ones. The float model, based on field experience and empirical trends, assumes that sampling biases arise from these dominant vertical drag dynamics. It does not include any of the several other effects of less certain significance, such as departures from lateral fluid velocities or lift resulting from the tilt of the drag plate in sheared flows. The types of simulated floats are as follows:
Ideal float model: These floats are always exactly ballasted with a density that matches that of the ambient water. However, their trajectories are not perfectly Lagrangian because the float model still includes virtual mass inertial terms.
Correct float model: These floats include realistic design features. The float density is set to that of the water (i.e., they are “ideal”) whenever the float is between 15 and 25 m, and otherwise it remains constant at its last value.
Uniform float model: These floats are the same as the “correct” model floats, but with additional upward buoyancy from mass deficits of 1, 2, 4, 8, 16, 32, or 64 g.
Bubble float model: These floats are the same as the “uniform” model floats at the surface, but with the upward buoyancy decaying with float depth −Z as (1 − Z/10 m)−1 corresponding to the gas law for attached air bubbles that do not dissolve. This mimics observed float buoyancy profiles.
Statistics are calculated on 100 model floats trajectories of each type, subsampled to 25 LES time steps to reduce storage and computational requirements. This preserves the high resolution of layer dynamics (at 2.5–20-s intervals) at ratios of sampling to mixing time scales that remain similar between the LES cases. These subsampled rates are still several times larger than in real Lagrangian float data. Floats are initially distributed uniformly to a depth greater than the mixed layer in each initial LES model domain. These statistics from each float model are identified using a short phrase. For example, “20 g uniform” implies the uniform model with 20 g of excess buoyancy, while “20 g bubble” implies the bubble model with 20 g of buoyancy at the surface. All four LES cases are sampled by ideal and correct model floats, as well as by the seven differently weighted uniform and bubble buoyancy model types, making 1600 model floats in each LES case. The high-wind cases H1 and H2 are additionally sampled by 20-, 40-, and 80-g floats of both uniform and surface-intensified bubble buoyancy type. Eulerian and float turbulence statistics were averaged either over the mixed layer or into bins of width in Z/HML, over a period of 300 000 time steps, except that case LW was averaged over only 250 000 time steps.
Special rules governing the float behavior near the surface are necessary to prevent the surface boundary condition w = 0 from trapping floats at the surface. Floats are excluded from the upper half of the grid layer by resetting the float depth to Z = −Δz/2 when this level is exceeded. This near-surface exclusion zone is quasi realistic because the large drag plates are centrally located about the O(1)-m-tall Lagrangian floats, and the mean float pressure is around this depth when the instrument begins to protrude from the surface. Buoyancy forces on the floats are extinguished gradually over the remaining upper LES grid layer: for −Δz/2 > Z > −Δz, the buoyancy force is scaled by 4(−Z/Δz − ½)2. The designation of the surface buoyancy of the bubble floats is therefore a slight misnomer, because the maximum buoyancy at Z = −Δz = −1.424 m is only 10 m/(10 m + Δz) = 0.8753 times the bubbles’ nominal surface buoyancy.
3. Simulated effects of excess float buoyancy
a. Float trajectories and the construction of mean profiles
Simulated depth–time float trajectories for three types of floats in case MW show the vertical extent of the float depth excursions, which qualitatively outline the depth of the actively mixing layer (Fig. 1). A more quantitative examination of the trajectories and their distribution, however, shows clear effects of the float buoyancy.
The floats’ motion relative to the water resulting from buoyancy forces affects the probability distribution function (PDF) Π(z) of their vertical position. For an ensemble of Ni floats at Zi(t), this mean PDF over an interval of time tI < t ≤ tF may be written in terms of the Dirac delta function δ,
where the (Ni dependent) unit normalization makes N the total number of float data points on the interval tF − tI, and Nz = NΠNi(z) the number at level z. Some ambiguity may be removed from this definition in the large ensemble limit,
and more generally, mean profiles constructed from any data Ξi(t) collected along float trajectories may be similarly defined as
Here, the superscript F distinguishes mean float profiles (ΞF) from Eulerian averages (ξ) taken uniformly over space or time without Π weighting. Upper- and lower-case symbols are used to distinguish float measurements (e.g., Z, W, Ξ, and Γ) from their Eulerian counterparts (e.g., z, w, ξ, and γ).
The definitions in Eqs. (2) and (3) still depend ambiguously on the initial distribution of the parcel-sized floats. The sum over the indexed float ensemble in the limit Ni → ∞ might also be specified equivalently as an integral over initial positions x0 in the initial volume Φ0 corresponding to tagged trajectories Z(t, x0). Here we study the statistics of mixed layer floats, by which it is implied that Φ0 is approximately the well-mixed upper-ocean layer at the beginning tI of an averaging period, and that when tF − tI is much larger than the integral mixing time scale of this boundary layer, vagaries in the precise definition of Φ0 become insignificant. We therefore omit explicit references to the initial ensemble distribution in our discussion of the mean distributions and float data statistics from mixed layer floats.
In practice, depth PDFs Π(z) are estimated from both synthetic and real float data by averaging contributions into finite depth bins Δz using N = NiNj data points at Zi(tj) from a finite ensemble of Ni floats sampled at Nj times on intervals tI < tj ≤ tF Mean profiles are correspondingly estimated from any simulated or real Lagrangian data Ξi(tj) observable by a finite ensemble of floats with limited spatiotemporal resolution,
where H is the Heaviside step function. Note that averages based instead on equally weighted contributions from the interpolation of trajectories onto gridded depths z may give different results that do not converge to Eqs. (1)–(3).
Because Lagrangian floats sample a volume that is expanding due to entrainment, their profiles are computed by averaging data from along the trajectories into bins of nondimensional depth z/HML, where HML is determined from the well-mixed depth in corresponding Eulerian profiles. In general, such profile averages will be indicated by an overbar, and are constructed here using bins of width Δz = 0.04HML.
In a large ensemble limit, mean turbulence statistics Ξi(tj) constructed from Lagrangian data will converge to the Eulerian average ΞF(z) = ξ(z) for “perfect” Lagrangian sampling trajectories that have both a uniform distribution Π(z) and zero bias velocity. These averages will diverge from Eulerian means when buoyancy effects cause the Lagrangian float measurement Ξ to differ from (hypothetically undisturbed) local ξ, and, more critically, create sampling that is both nonuniform and coherent with variations in ξ. These measurement biases may also arise when the bias velocity is zero but the distribution is not uniform, as occurs when a float ensemble within a mixed layer expands during entrainment with WF = −wF = −(z/HML)dHML/dt. However, such contributions from layer deepening are generally much smaller than those resulting from the float buoyancy effects studied here.
b. Profiles of vertical velocity bias
The mean profiles of bias velocity are constructed from simulations of case MW using Eq. (4). The LES-resolved fluid velocity along trajectories for float models with mass deficits less than 8 g (Fig. 2) show that the buoyancy of the floats corresponds closely, but not exactly, to their upward vertical motion relative to the water. When the mass deficit is small and the drag force is large, a float’s fall rate approaches its terminal velocity Wterm, related by
to the float buoyancy mass mb, the water density ρ, the horizontal area A = 0.5 m2, and the drag coefficient CD = 1. The ideal floats have no relative buoyancy, the correct floats have very little within the mixed layer, the uniform floats have a nearly uniform buoyancy offset, and the bubble floats have an upward buoyancy that varies along trajectories resulting from a depth-dependent mb that increases toward the surface (Fig. 2).
Figure 2 shows the relationship between and the corresponding profile of mean terminal velocity from Eq. (5) shown in gray. In the top bin, both and are reduced because model float buoyancy is extinguished adjacent to the excluded zone above z = −Δz/2, but because of contributions at z = −Δz/2 in updrafts and weak downdrafts. The difference between and disappears with time as floats are pulled down into the mixed layer. New departures from the terminal velocity can arise at depth along simulated float trajectories resulting from model float buoyancy changes to variations in pressure, temperature, and salinity along simulated trajectories, and the LES subgrid diffusion of momentum. Correct and bubble model floats that sample strong density variance in the lower layer show such departures from temporarily mixing into and rising out of denser entrained fluid.
The overall extent of simulated departures in and can depend slightly on forcing and depth scales in the LES models, but their layer averages are within 5% of each other across all simulations here for mass deficits larger than 1 g. Real Lagrangian mass deficits and drag properties are often uncertain, but recent design generations are estimated to have 3 times larger ACD. Their bias velocities are therefore expected to differ even less significantly from . Model results will show how systematic errors in higher turbulence moments can be predicted and corrected as functions of However, uncertainties and variability in any individual Lagrangian float’s drag and buoyancy make using as a predictor of these effects less useful than several other equivalent predictors introduced in section 4.
c. Depth distribution
The motion of the floats relative to the water affects their depth distribution. Figure 3 shows typical depth distribution for different classes of floats. Ideal floats fill the entire mixed layer so that their distribution therein is nearly uniform. It drops abruptly at the bottom of the mixed layer, here defined as the depth above the pycnocline at which the mean stratification profile falls to 10% of its maximum value within the pycnocline. Some ideal floats penetrate into the pycnocline within the deeper mixing layer, accurately reflecting the continual exchange of fluid between the mixed layer and the pycnocline. The relatively rapid return of the ideal floats to the more homogeneous mixed layer is enhanced by steady entrainment deepening in these simulations. For example, near day 1.4 an ideal float leaves the mixed layer, lingers in the pycnocline for about 0.1 day, and is then entrained back into the mixed layer (Fig. 1). Examples of this behavior in real floats are shown in D’Asaro (2003b, their Fig. 3). For a sufficiently long simulation with some mixing throughout, ideal floats will eventually be mixed uniformly throughout the entire simulation volume. The nonuniform distributions in these simulations are artifacts of starting floats only within the mixed layer.
Correct floats match their buoyancy to that of the mixed layer only in the 15–25-m depth interval. They are distributed uniformly within the mixed layer, but spend much less time below it. Recall that the density of the boundary layer is continually increasing because of entrainment. Correct floats that remain within the mixed layer, and thus frequently revisit the 15–25-m layer, will have a density that matches that of the mixed layer closely and thus will behave much like ideal floats. However, those that remain in the more quiescent region near the bottom of the mixed layer or in the upper pycnocline will gain buoyancy and thus tend to move upward out of this region. This creates a low bias in the number of correct floats at the bottom 20% of the layer. Nevertheless, the ballasting strategy of the correct floats appears to be effective in keeping floats within the mixed layer while sampling it nearly uniformly.
Uniform floats are buoyant at all times and depths and thus continually move upward relative to the water. Their depth distribution increases toward the surface with a greater surface bias for more buoyant floats. However, even the most buoyant floats reach the bottom of the mixed layer occasionally. D’Asaro et al. (2002) and Harcourt et al. (2002) examine the behavior of such floats. Buoyant floats concentrate near the surface because they move upward relative to the water. Classically, the upward drift Wbias of the floats is balanced by an effective vertical eddy diffusivity KF of the floats in the boundary layer turbulence. Modeling the PDF of floats as a scalar concentration with mean profile Π(z) and assuming horizontal homogeneity leads to
If Wbias and KF are constant, then this advection–diffusion equation yields an exponentially decaying steady-state solution
where λ = KF/Wbias = Π−1dΠ/dz and Π0 = λ−1 gives the integral of Π(z) over z < 0 unit normalization. This is a simplified formulation, but it captures the essence of the physics determining Π(z). The decay in concentration Π(z) with depth increases with (Wbias > 0) bias velocity and decreases with the strength (KF) of mixing. It is used here to guide our analysis.
Bubble floats have extra buoyancy relative to correct floats near the surface. This concentrates the PDF for bubble-type floats near the surface (Fig. 3b) where they are buoyant, but retains a nearly uniform distribution at depth, where they are more neutral. These results are qualitatively consistent with a local interpretation of Eq. (7) in which Wbias and dΠ/dz are large near the surface and small at depth where the PDF is more uniform.
d. Notation for bulk averages
Bulk averages representing mean statistics over the mixed layer are indicated by the use of brackets 〈 … 〉 to distinguish them from the profiles of horizontal means denoted by overbars. A subscript further distinguishes the bulk Eulerian average of a quantity Ξ
where ζ = z/HML from its bulk Lagrangian average
over the time series of data Ξi from i = 1, … , N mixed layer float trajectories of duration tF − tI.
Even buoyant Lagrangian float trajectories provide a visually compelling image of the mixing depth (Fig. 1), but quantitative estimates can be sensitive to small biases. One simple measure of the mixed layer depth is twice the average depth of floats initially deployed within it. Such a bulk mixed layer average is Habs = −2〈Z〉F. If the depth distribution of Z were perfectly uniform extending from the surface to depth HML, and zero below, then 〈Z〉F = 〈Z〉E = Habs and Habs = HML. The value of Habs is plotted for each case in Fig. 3. For the ideal floats Habs is slightly larger than HML, because the distribution extends below HML; for correct floats, the correspondence is good. For increasingly buoyant floats, Habs is increasing shallower than HML.
e. Buoyancy effects on VKE and skewness
Float buoyancy affects turbulence statistics computed from the floats. Profiles of vertical velocity variance WWF (Fig. 4) and vertical velocity skewness
(Fig. 5) computed from ideal and correct floats match closely the corresponding Eulerian statistics [ww and SWE(z) = W3/W23/2, respectively] computed from the resolved LES model fields that advect the floats. Only small differences arise because of float ensemble sizes, interpolation errors, and the small buoyancy effects that keep the correct floats within the mixed layer. Vertical velocity variance increases to a maximum in the upper mixed layer and then decays toward the bottom of the mixed layer with peak values well above u*2, the wave-free surface friction velocity scaling. These high values are due to surface wave forcing in addition to stress forcing, as parameterized in HD08. Skewness is negative, about −0.7, and relatively uniform, indicating a dominance of larger, less frequent downward velocities and weaker, more frequent upward velocities at all depths.
Buoyant floats, with either uniform or bubble buoyancy, show an enhanced level of velocity variance and decreased skewness in the region of buoyancy at all depths for the uniformly buoyant floats and near the surface only for the bubble floats. These effects come about because buoyant floats oversample the near-surface convergence zones into which they are swept and the downdrafts beneath them. Here, these are the convergent surface jets above sheet-like downdrafts between counter-rotating turbulent Langmuir cells. Buoyant floats therefore do not encounter more quiescent fluid as often at depth, and when they do it is sampled for shorter durations because they tend to rise out. These behaviors cause float velocity variance to be high for significantly buoyant floats.
While buoyant floats oversample the structures responsible for strongly negative fluid velocity skewness, the corresponding float measurement is reduced in magnitude. This is because while rising against these downdrafts, they record reduced downward velocities; and while passing through updrafts on their return to the surface, their velocities are higher. The magnitude of W skewness is reduced toward zero by both this shift in reported velocity and the inherently related shift in sampling between up- and downdrafts that is counter to their skewness in the Eulerian distribution.
f. Buoyancy effects on the Lagrangian budget of W
In the absence of directly measured float biases, the effects of float buoyancy on the Lagrangian budget of vertical velocity can also be identified in tests for self-consistency in the Lagrangian budget of vertical velocity (Harcourt et al. 2002; D’Asaro et al. 2002). For Eulerian averages or perfect Lagrangian sampling, a profile of the mean rate of change Dw/Dt ≡ ∂w/∂t + in vertical velocity along float trajectories will be related in a steady-state, horizontally homogeneous environment to the VKE profile ww (Fig. 4) by Dw/Dt = ∂ww/∂z. Departures γw = Dw/Dt − ∂ww/∂z from this equivalence may be due to unsteadiness, horizontal inhomogeneity, or statistically biased sampling. For nearly Lagrangian float trajectories this divergence of Lagrangian and Eulerian momentum budgets appears as ΓW ≠ 0, where the momentum budget difference
is generalized by using the mean profile DFW/DtF of float accelerations along imperfectly Lagrangian trajectories. Where the mixed layer is horizontally homogeneous and evolving slowly, a nondimensional measure of buoyancy error can be defined on the integral of this discrepancy in the vertical velocity budget
Computed in layer-relative depths, EW(z/HML) profiles from the LES case MW show the steady increase of this measure with float buoyancy, particularly at middepth (Fig. 6). The middepth value EW(z = −HML/2) was used in D’Asaro (2004) as one measure of float error resulting from excess buoyancy.
The relationship between the W budget discrepancy ΓW and the effect of float buoyancy on the depth PDF Π(z) becomes better defined when it is noted that ΓW is nearly equal to
in simulated data. Substitution for ΓW on the left-hand side for the Lagrangian budget in Eq. (12) gives an alternate nondimensional error estimate
The apparent equality of
That ΓW and ΓΠ converge to give Eq. (16) in the steady-state limit is derived in the appendix as a special case of a more general relationship for float budgets of conserved quantities. Note that Eq. (16) is a potentially useful relationship between buoyantly biased observations, but it does not relate them to unbiased Eulerian statistics and therefore cannot be used directly to remove biases from float data.
The effect of local excess float buoyancy and Wbias is evident in the slope of EW (Fig. 6a,b), which is related by Eq. (14) to the value of λ in Eq. (7). However, pursuit of a corrective prescription for profiles of float statistics showed that the effect of Wbias also appears to have nonlocal effects within the float VKE and W skewness profiles. This motivated a search for proxies for the layer-averaged relative bias Wbias/〈WW〉F1/2 to predict the effect of float buoyancy on layer-averaged W statistics. Of such additional predictors of statistical bias based on ΓΠ and related quantities, the best candidate encountered was the PDF-weighted 〈ΓΠ〉F = 〈ΠΓΠ〉E, scaled on 〈WW〉F/HML.
4. Predictors for errors in bulk W statistics for buoyant floats
Vertical kinetic energy errors are measured by δW2/〈W2〉E, where δW2 = 〈WW〉F − 〈ww〉E. Vertical velocity skewness errors are computed in the following two ways: δSW〈 〉 = SW〈F〉 − Sw〈E〉 is the error in estimating the skewness Sw〈E〉 = 〈w3〉E/〈w2〉E3/2 of bulk Eulerian w from the bulk float W skewness Sw〈F〉 ≡ 〈W3〉F/〈W2〉F3/2; δ〈SW〉 = 〈SWF〉E − 〈SwE〉E is the difference between the uniformly weighted layer average 〈SWF〉E of the float skewness profiles SWF [see Eq. (10) and Fig. 5] and the corresponding mean of the Eulerian profile SwE. Similar numerical results are given for fractional errors δSW〈 〉/Sw〈E〉 and δ〈SW〉/〈SwE〉E in both measures of W skewness. While the measurement of the profile (Fig. 5) of skewness is a more recognizable metric of mixed layer structures, the former is preferred for float measurements because it underweights contributions of large uncertainty from the deepest levels with the least data.
b. Wbias as a predictor
The most straightforward corrections for errors in bulk VKE and skewness are those based on the trajectory average 〈Wbias〉F of a directly measured float bias Wbias. Figure 7a plots Y = δW2/〈WW〉F against the predictor X = 〈Wbias〉F/〈W2〉F1/2. A linear relationship holds between VKE error and X for X < 0.75. At larger X, Y decreases as the buoyant floats become increasingly trapped at the surface and thus have small values of W.
A more linear relationship is found in Fig. 7b for the error in Sw〈F〉 for biases that are up to twice as large. The larger spread in VKE error dependence between LES model cases may depend on the different rates of entrainment between the lower and higher wind cases, and the corresponding differences in the extent to which floats penetrate the pycnocline.
c. Z skewness and other measures as predictors
Practically, Wbias is often difficult to estimate. We have therefore evaluated five less direct measures of float buoyancy. Of these, the skewness SZ′F = 〈Z′3〉F〈Z′2〉F−3/2 of Z′ = 〈Z〉F − Z was found to be the best predictor, as well as simple to evaluate. It is well correlated with 〈Wbias〉F/〈W2〉F1/2 (Fig. 8b) and δW2/〈WW〉F (Fig. 9a). As in Fig. 7a, δW2/〈WW〉F increases, but then saturates and decreases with increasing values of SZ′F. Note also that SZ′F is slightly positive, even for ideal and correctly ballasted floats; that is, the Y intercept is nonzero in Fig. 8b or in the last row of Table 1. This represents the ongoing effect on the depth PDF Π from the displacement of mixed layer floats by recently entrained fluid.
All five quantities and their relationship to nondimensional Wbias are detailed in Table 1. These additional predictors of statistical bias are EW(z = −HML/2), HML〈ΓΠ〉F/〈WW〉F, and three quantities based on moments of Π:3〈Z2〉F〈2Z〉F−2 − 1, 2 〈|Z′|〉F|〈Z〉|F−1 − 1 and SZ′F. Table 1 lists the linear fits between these and 〈Wbias〉F/〈W2〉F1/2 for 〈Wbias〉F/〈W2〉F1/2 < 0.5. Table 2 lists the linear fits between each quantity and the normalized errors δW2/〈W2〉E, −δSW〈 〉/Sw〈E〉, and −δ〈SW〉/〈SwE〉E. The range for fitting was determined subjectively, and in several fits there were outliers within this range from the most extreme buoyancy cases (e.g., 64 g buoyant in case LW) that had to be discarded manually. As in Figs. 7a,b and 9a,b, there is no significant difference between the scaling of bulk errors from the uniform and bubble float models within each LES case set. There is also little notable difference in scaling behavior at small biases between different LES cases, weighted equivalently in these fits. However, Lagrangian float data shown below in section 5 suggest such sensitivity, possibly for cases with different surface buoyancy fluxes.
The relative errors in the slopes of linear fits suggest that 〈Wbias〉F/〈W2〉F1/2 is the best predictor of variance and skewness errors (i.e., only 2% uncertainty for slope fit in Fig. 7b). The second most accurate (∼3.5%) is Sz′F (lowest row in Tables 1 and 2), and the least accurate (5%–6%) is EW(z = −HML/2). Predictive skill for W variance error is more uniform between the six different measures of float buoyancy, reflected by relative errors of 4%–5% for all slopes except for a 5.4% uncertainty for EW(z = −HML/2). Fitted lines all have intercepts near zero error, but given the small but fundamental differences between Eulerian and Lagrangian statistics in even the correctly ballasted float cases, it is not surprising that there is some difference at zero ballast error. Differences in the uncertainty of the slopes of the fits here appear to stem from differences between the LES cases, rather than from the size of the float sample or the duration of the LES simulation relative to turbulent integral time scales. In applying such fits to recover Eulerian bulk statistics from Lagrangian float data, an estimate of uncertainty introduced (in addition to sampling error) may be obtained by rescaling the Table 2 uncertainties bracketing the linear fit Y = [m ± σ(m)]X + [b ± σ(b)] (and only about half the data points) so that its standard deviation from the mean fit over the ordinate range matches that of the data σ(ΔY), which is generally larger.
5. Comparisons with float data
a. Does depth skewness predict Wbias in data?
A 50-day-long (April–June 2008) deployment of a second-generation mixed layer float (MLFII; D’Asaro 2003a) in the North Atlantic south of Iceland (nominally 61°N, 21°W) was part of the 2008 North Atlantic Bloom Experiment. On each day, the float cycled through periods of vertical profiling, satellite communication, autoballasting, and an approximately 22-h-long Lagrangian drift. During these drifts, the float was in the mixed layer, nearly neutrally buoyant, and was repeatedly carried across the boundary layer by the turbulence in a manner very similar to that shown in Fig. 1. The float depth was measured using pressure every 50 s and the vertical velocity was computed from the variations in depth.
During the drifts, the float continuously adjusted its volume using the correct float method so as to retain a constant buoyancy based on measurements from two onboard CTDs and analytical expression for its own equation of state and that of seawater. The float’s equation of state was modeled as
(D’Asaro 2003a), where V is the float volume, B is the volume of the ballasting piston, V0 is a reference volume, T is temperature, P is pressure above a reference atmospheric pressure Pat = 10 dbar, α is the float’s thermal expansion coefficient, γ is its compressibility, and A is the amount of trapped air. The value of α was set from the known materials of the float and B was varied to control the float. During its daily autoballasting period (not shown in Fig. 10) the float settled on an isopycnal within the pycnocline. The value of V0 was computed from each of these periods using Eq. (17), assuming that the density of the float matched that of the water as measured by the CTDs. Isopycnals at different depths were chosen on different days, so that over time the values of V0, γ, and A could be determined. These measurements showed that after yearday 105, the float had a stable equation of state, with Eq. (17) describing the volume to an accuracy of about 10−6 m3 (1 cc). Thus, the float’s buoyancy is known throughout this period to about 1 g. However, this stability was not apparent in real time and the overly aggressive daily adjustments of the ballasting parameters resulted in a float buoyancy that varied from nearly neutral to as much as 30 g. Fortuitously, this provides a large dynamic range in float buoyancy, which is used here.
A plot of the stratification, float trajectory, wind speed, and float buoyancy for this deployment shows that during stronger winds and more neutral buoyancy, for example, on yearday 122, the float oscillated across the mixed layer as expected (Fig. 10). Early in the mission, the float was often trapped at the surface during times of weak wind and positive float buoyancy, for example, on yearday 105. Later, the buoyancy control was adjusted so that it settled near 10-m depth for safety, for example, on yearday 137. Negative buoyancy, for example, on yearday 143, resulted in the float settling onto an isopycnal below the boundary layer.
The value of Wbias was computed from the buoyancy using Eq. (5). During drifts the drag is primarily due to a hexagonal cloth drogue extending from the middle of the float with an area of about 1 m2. Approximating the drogue as a rigid plate with flow perpendicular to the plate, CD ∼ 1.1 so that CDA = 1.1 m2. A databased value of CDA was estimated from plots of vertical velocity against buoyancy similar to those in D’Asaro (2003a, Fig. 11 therein) yielding CDA = 1.5 m2. The larger value is easily justified by the additional drag resulting from the float itself. The estimate is based primarily on periods when the float was moving downward. During periods of upward motion faster than 0.2 m s−1 the drag appears to be significantly less, probably because the drogue folds under force of the flow. However, this does not happen at the low relative speeds during drift mode. Accordingly, Wbias was computed using CDA = 1.5 m2.
The skewness of perturbation float depth SZ′F = 〈Z′3〉F/〈Z′2〉F3/2 was computed from the measured float pressure for each drift period. The perturbation depth Z′ = 〈Z〉F − Z and the resulting skewness are sensitive to how the mean 〈Z〉F is computed. As in Fig. 1, a float trajectory can linger near the surface or at depth for extended periods. These variations may be due to real variations in the boundary layer depth (e.g., shoaling in response to diurnal heating) or merely random. Using a short time interval to compute 〈Z〉F causes the mean to track these variations, resulting in a more Gaussian distribution around the mean and smaller skewnesses. Using a longer time interval reduces this effect but may lead to anomalously large skewnesses resulting from real variations in the boundary layer depth. Here, each drift is divided into blocks approximately 12 ks (1 ks = 103 s) long, and the mean is computed in each block. Skewness is computed using this stepwise-varying mean. Error in each skewness estimate is computed from Monte Carlo simulations of the standard deviation of skewness for a Gaussian distribution as function of degrees of freedom. The degrees of freedom for each drift are estimated assuming a correlation time of 2〈Z〉F/〈W2〉F1/2, noting that the mixed layer depth is given by about 2〈Z〉F.
Figure 11 tests the predictions of Fig. 8b using these data. A clear relationship between scaled Wbias and skewness is evident, either for all data (thin line) or only including data with Wbias > 0 (thick line). The results agree well with the predictions from the simulations (dashed line).
b. Does depth skewness predict δW2 in data?
Existing float data (e.g., D’Asaro 2001) show a strong empirical correlation between 〈WW〉F and the wind stress ρu*2 computed using bulk formulas. Accordingly, the deviations from this relationship can be used to assess any additional effect of float buoyancy on the measured value of 〈WW〉F.
Figure 12a shows 〈WW〉F u*−2 as a function of Z skewness for the same float data as in Fig. 11. The friction velocity u* is computed from an operational National Centers for Environmental Prediction (NCEP) model wind product obtained through the National Oceanic and Atmospheric Administration (NOAA)/National Weather Service (NWS), interpolated to the float location using the neutral 10-m Large and Pond (1981) drag coefficient. Each point is averaged over a single drift. A cubic fit (thick line) to heavy symbols shows an increase in estimated 〈W2〉 with Z skewness at low SZ′F, saturating near SZ′F ≅ 1, and then decreasing; this is –the same trend that is revealed by the simulations (Fig. 9a).
A second set of float data provides a more robust relationship. D’Asaro (2001) describes ship-based deployments of three Lagrangian floats for 2 weeks in deep water off of Vancouver Island, British Columbia, Canada, during January 1995. The acoustically tracked floats (D’Asaro et al. 1996) were each deployed for approximately 24 h, and then recovered and redeployed at the center of the tracking array. The floats were reballasted by hand before each deployment based on changes in the water density measured by the shipboard CTD and qualitative examination of their depth PDF. Typically, two or three floats were in the water simultaneously, separated by either a few kilometers for floats launched separately or by tens to hundreds of meters for float groups launched together. The presence of multiple nearby floats with the same meteorological forcing and mixed layer depth makes the computation of skewness more robust than for records using a single float.
Skewness was computed using all float data in a given time bin. The bin width is necessarily a compromise between reducing statistical error, requiring longer bins, and sampling over a statistically homogeneous region, requiring shorter bins. The cleanest result, found using half-overlapping 240-ks (2.8 day) intervals, is shown in Fig. 12b. The components are shown in Fig. 13. The vertical velocity variance is nearly proportional to u*2 (or wind speed squared) and is thus largest during the two windy periods near yeardays 19 and 28. The skewness decreases at the higher winds consistent with the scaling in Fig. 11. Note that the skewness is larger at the same wind speed during the second windy period, presumably because of systematic changes in the manual ballasting accuracy. Using the scaling in Fig. 11, estimates of Wbias vary from 5 to 30 mm s−1, similar to those estimated in D’Asaro (2001), generally increasing with time.
The two panels of Fig. 12 show the same pattern of increasing, saturating, and then decreasing 〈WW〉F u*−2 with increasing Z skewness SZ′F. Both panels show an extrapolated value of 〈WW〉F u*−2 of about 1 for zero skewness and a peak near unit skewness. These are sufficient to support the claim that skewness is a useful indicator of float buoyancy and the resulting measurement biases. However, the magnitude of the peak is significantly higher for the North Atlantic data than for the North Pacific data. This may be an artifact of the large errors in the North Atlantic data. If not, it may represent differences in the properties of the two different float types, differences in the manner with which u* was computed, differences in the averaging interval used, differences in the variability of the mixed layer, or differences in wind, wave, or other forcing at the two locations. These issues are the subjects for future studies.
A set of large-eddy simulations spanning a wide variety of realistic wind and wave cases were used to study the ability of imperfect Lagrangian floats to measure the turbulent properties of the upper-ocean mixed layer. The following three properties were studied: the depth of mixing computed from the depth of float mixing, the vertical velocity variance 〈WW〉F computed from the float velocity, and the skewness of vertical velocity computed from the float velocity.
Correctly ballasted Lagrangian floats that match buoyancy within the mixed layer yield excellent estimates of these turbulent statistics, with biases of less than 10%. Such floats are distributed nearly uniformly across the mixed layer.
Excess float buoyancy causes model floats to move upward relative to the water with average speed Wbias close to the terminal velocity Wterm and have a depth distribution that increases toward the surface. Various measures of this increase are useful predictors of Wbias〈WW〉F−1/2. Of these, the skewness of float depth is the most useful because of its computational simplicity and good predictive skill.
Buoyant floats bias the mixing depth low, the W variance high, and the W skewness high. The W variance increases linearly with Wbias〈WW〉F−1/2 and Z skewness until they reach about 1. Here, the variance estimated from the floats is about 160% of the true value. For higher values W variance saturates or decreases as the floats become increasingly trapped near the surface. The W skewness increases monotonically with Z skewness.
The model predictions were tested using data from two float datasets. The Z skewness was found to be a good measure of Wbias〈WW〉F−1/2., as predicted. The measured W variance scaled by friction velocity 〈WW〉F u*−2 increases, saturates, and decreases with increasing Z skewness, as predicted for both datasets. For typical (10 m s−1) wind speeds, accurate estimates of Z skewness require several float days of averaging. Natural variability of the boundary layer on these time scales can bias the Z skewness high, so that even neutrally buoyant floats can have a positive Z skewness.
These results show the importance of accurate float ballasting for making accurate measurements of turbulence statistics within the mixed layer. They provide a mechanism for correcting biases in imperfectly ballasted floats resulting from float buoyancy and open the possibility of autonomously ballasting floats in the boundary layer.
This work was supported by the National Science Foundation (OCE0549887 and OCE0628379), the Office of Naval Research (N00014-08-1-0446, N00014-09-1-0174, N00014-08-1-0575, N00014-08-1-0447, and N00014-08-1-0577), and by a grant of HPC resources from the Department of Defense High Performance Computing Modernization Program.
Relationships between Biased Measurements
The exact form of the steady-state relationship ΓW = ΓΠ [Eq. (16)], between the PDF of depth Π and the float-weighted mean profiles WWF and DFW/DtF, may be obtained as follows: dropping for brevity from Eq. (3) the notation of the large ensemble limit, as well as the tags on the ith trajectory Zi(t), the profile of covariance between W(t), and any float observable Ξ(t) is
Because Nz = NΠ,
The generalized scaling property of the Dirac delta function is
with summation over the Ni roots f (ti) = 0. Substituting
in the integrand of Eq. (A2) yields
Because , the right integral is ΠDFΞ/DtF, using the modified material derivative notation DF/Dt for the time derivative along the float trajectory. The result
can also be stated in terms of Eulerian horizontal averages weighted by float concentration Π,
Corresponding author address: R. R. Harcourt, 1013 NE 40th Street, Applied Physics Laboratory, University of Washington, Seattle, WA 98105. Email: email@example.com