## Abstract

Observations of surface waves, currents, and turbulence at the Columbia River mouth are used to investigate the source and vertical structure of turbulence in the surface boundary layer. Turbulent velocity data collected on board freely drifting Surface Wave Instrument Float with Tracking (SWIFT) buoys are corrected for platform motions to estimate turbulent kinetic energy (TKE) and TKE dissipation rates. Both of these quantities are correlated with wave steepness, which has previously been shown to determine wave breaking within the same dataset. Estimates of the turbulent length scale increase linearly with distance from the free surface, and roughness lengths estimated from velocity statistics scale with significant wave height. The vertical decay of turbulence is consistent with a balance between vertical diffusion and dissipation. Below a critical depth, a power-law scaling commonly applied in the literature works well to fit the data. Above this depth, an exponential scaling fits the data well. These results, which are in a surface-following reference frame, are reconciled with results from the literature in a fixed reference frame. A mapping between free-surface and mean-surface reference coordinates suggests 30% of the TKE is dissipated above the mean sea surface.

## 1. Introduction

Turbulence at the ocean surface is important to the exchange of gasses, heat, momentum, and kinetic energy between the atmosphere and ocean. Turbulence introduced through wave breaking (e.g., Craig and Banner 1994) as well as wave–turbulence interactions (e.g., Thorpe 2004) elevate turbulence levels beyond the predictions for classic rigid boundary layers (Agrawal et al. 1992). Extensive work over the past three decades has improved understanding of the oceanic surface boundary layer through field measurements (Agrawal et al. 1992; Terray et al. 1996; Drennan et al. 1996; Gemmrich and Farmer 2004; Jones and Monismith 2008; Gerbi et al. 2009; Sutherland and Melville 2015; Thomson et al. 2016), development of models (Craig and Banner 1994; Burchard 2001; Umlauf et al. 2003; Umlauf and Burchard 2003; Carniel et al. 2009), and laboratory measurements (Duncan 1981; Rapp and Melville 1990; Lamarre and Melville 1991; Drazen and Melville 2009). This prior work has focused on wave conditions in deep and intermediate water depth. A more limited literature has focused on measurements and models in the surf zone (Feddersen and Trowbridge 2005; Feddersen et al. 2007; Feddersen 2012a,b; Grasso et al. 2012) and at river inlets (Thomson et al. 2014; Zippel and Thomson 2015; Moghimi et al. 2016), where modifications to the wave field from currents and bathymetry alter surface boundary layer processes.

Very near the surface (within a few wave heights), turbulent kinetic energy (TKE) dissipation rates are balanced to a first order by turbulent transport (Scully et al. 2016), which can be modeled as a diffusive process (Craig and Banner 1994; Burchard 2001; Umlauf and Burchard 2003). Wave breaking provides a source of turbulence, which is modeled as a TKE flux input at the surface. In deep water, the equilibrium of short wind waves (Phillips 1985; Thomson et al. 2013) is often assumed, and the surface flux into the ocean is estimated from wind parameters (Gemmrich et al. 1994; Terray et al. 1996; Sutherland and Melville 2015; Thomson et al. 2016). In the surf zone, the breaking of long waves injects additional TKE to the surface, on the order of 10%–15% of the total wave energy flux gradient (Feddersen 2012b; Zippel and Thomson 2015).

The surface flux of turbulence is difficult to prescribe at river inlets, where wave breaking is different from purely wind-driven whitecapping or depth-limited surf. At river inlets, strong currents and gradients in currents can shoal and refract surface waves, often causing breaking in intermediate depth (Zippel and Thomson 2017). Indeed, even wave dissipation (distinct from the turbulent dissipation) in such environments is still an active area of research (e.g., Rapizo et al. 2017).

In addition to the magnitude of the TKE surface flux from wave breaking, the vertical fate of this turbulence remains an active research area. Many studies agree that the decay scale is set by the significant wave height *H*_{s} and that the vertical decay is a power law. However, measurements have yet to converge on a single decay exponent *λ* for TKE dissipation rate. Estimates are typically constrained to 1 < *λ* < 2, but this appears to be sensitive to the choice of reference frame. Many studies using fixed frame instruments, such as Terray et al. (1996), Drennan et al. (1996), Jones and Monismith (2008), and Gerbi et al. (2009), found decay scales of *λ* = 2. Studies measuring inside wave crests with wave-following platforms found different values, for example, Gemmrich (2010) found 1 < *λ* < 1.5, Sutherland and Melville (2015) found *λ* = 1 for *z* < 0.6*H*_{s}, and Thomson et al. (2016) found *λ* = 1.4.

There is a lack of consensus, then, on the appropriate surface flux of TKE and its vertical decay at river inlets. This has, in part, lead to difficulty in understanding how wave-breaking turbulence influences such regions and whether wave-breaking turbulence has distinct properties in these regions. Certainly, the bathymetry and currents at river inlets can enhance wave breaking, but once the waves have broken, the resulting turbulence may not be any different than it is in the open ocean. There is a small, growing body of work on how wave-breaking turbulence might interact with buoyant layers. For example, Gerbi et al. (2013) modeled a buoyant river plume during upwelling-favorable winds and found that the inclusion of wave-breaking turbulence increased plume thickness. Using field measurements, Thomson et al. (2014) showed large wave energy flux gradients across a plume front and observed wave-breaking turbulence levels at the surface that were as large as published turbulence values at the subsurface plume interface. Further studies have investigated surface boundary layer effects where buoyancy is relevant (Vagle et al. 2012; Gerbi et al. 2015).

### Turbulence scalings

Craig and Banner (1994) presented one of the first analytic results for wave-breaking turbulence. Assuming turbulent transport can be modeled diffusively, they presented a solution where vertical diffusive transport balances TKE dissipation near the ocean surface,

where *ν*_{k} is diffusivity of TKE, *q*^{2} is TKE, and *ε* is the TKE dissipation rate. In Craig and Banner (1994), diffusivity is a function of turbulent length scale , a constant *S*_{q}, and the turbulent velocity scale such that . We note that it has become more common since Craig and Banner (1994) to express the diffusivity of TKE as a function of the eddy diffusivity and the turbulent Schmidt number *σ*_{k}, such that (Umlauf and Burchard 2003), where *C*_{μ} is a shear-dependent stability function (Canuto et al. 2001). Craig and Banner (1994) also assumed a linearly increasing turbulent length scale with distance from the surface *z* (defined positive down):

where *z*_{0} is a roughness length, and *L* was assumed to equal Von Kàrmàn’s constant *κ*, such that *L* = *κ* = 0.41. The common closure assumption (sometimes called the cascading relation) between TKE, , and the TKE dissipation rate *ε* was also made:

where is a model constant. Last, using a surface boundary condition relating the flux of TKE through the surface *G* to wave energy dissipation, Eq. (1) was shown to yield power-law decay functions for TKE dissipation and TKE, respectively:

Umlauf and Burchard (2003) showed the power-law decay *λ* could be expressed in terms of model parameters,

where . In the absence of mean shear (as is expected in the surface boundary layer with wave breaking), the stability function *C*_{μ} is assumed to equal the model constant , such that *R* = 1 (Canuto et al. 2001; Umlauf and Burchard 2003). The constants used in Craig and Banner (1994) result in *λ* = 3.4, which would require *σ*_{k} ≈ 2 under this framework [Feddersen and Trowbridge (2005, their section 2c) offers a useful discussion on *σ*_{k}]. Burchard (2001) therefore suggested the turbulent Schmidt number be a function of the ratio of production and dissipation, *σ*_{k} = *σ*_{k}(*P*/*ε*).

Terray et al. (1996) proposed that the surface roughness length is proportional to the significant wave height, *z*_{0} ∝ *H*_{s}, and used field measurements from small waves on a lake to fit the scaling,

A power-law decay with *λ* = 2 was proposed to hold below a breaking layer with depth, *z*_{brk} = 0.6*H*_{s}, above which dissipation rate was assumed to be constant. This constant breaking layer concept has been refuted with recent measurements from surface-following platforms, where the decay slope 1 ≤ *λ* ≤ 2 is found very near the free surface (Gemmrich 2010; Sutherland and Melville 2015; Thomson et al. 2016). These recent measurements are in moving frame coordinates , where *η* is the ocean free surface.

Currently, there is no clear consensus on how to map measurements from the coordinate frame to coordinates referenced to the mean sea surface; that is to say, mapping from to . Both Gemmrich (2010) and Thomson et al. (2016) did it directly, using raw time series of *η*. Without a general coordinate transform, it is difficult to fully interpret comparisons of the various field measurements and model predictions. It is also important to note the change from (*z*/*z*_{0} + 1) in the analytic solution to (*z*/*H*_{s}) in the scaling, which can give similar functional values near the surface for different values of *λ* and *z*_{0}/*H*_{s}.

Choosing a constant turbulence length scale, results in an exponential decay, rather than a power-law decay, as discussed in Umlauf et al. (2003) for deep water, and independently in Feddersen (2012b) for shallow water. This solution is particularly interesting as it may apply in the region *z* ≤ *z*_{0}, where the assumption of a linearly increasing length scale near the surface may not hold. Following Feddersen (2012b), the TKE dissipation rate could then be expressed as

with and , such that specification of sets the decay scale.

Here, we present field measurements of turbulence and waves from the mouth of the Columbia River to examine the validity of these surface turbulence models under a wide range of wave conditions. The uniqueness of the river mouth, relative to the open ocean, remains an open question, but the practical effect is to provide a natural laboratory with ample wave breaking. We focus in particular on determining an appropriate model roughness length and length-scale decay constant for the surface turbulence. A description of the field site, the dataset, and wave and turbulence processing techniques are presented in section 2. Data processing includes a method for correcting buoy velocities for platform motion and compares two methods for estimating TKE dissipation rates. Field measurements are compared with existing open-ocean turbulence models in section 3, along with a limited exploration of the interaction of surface turbulence with the subsurface stratification. Section 4 discusses the choice of model constants, and the implications of the measurement reference frame on the results. The results are summarized in section 5.

## 2. Methods

Measurements of waves and turbulence were collected from freely drifting Surface Wave Instrument Float with Tracking (SWIFT) buoys (Thomson 2012) at the mouth of the Columbia River as part of the Riverine and Estuarine Transport (RIVET) Experiment II (RIVET-II) between April and September of 2013. Up to six buoys were deployed at a time on drifts lasting a few hours each. On ebbs, the research vessel would often wait for the tide to change in order to safely cross the Columbia River Bar, such that buoys were not tended during the drifts. On floods, the research vessel could operate throughout the domain, and buoys were tended during the drifts (including being reset if they approached shore, thus avoiding grounding). Therefore, ebb deployments lasted longer, and more data were collected on ebb tides. Buoys were deployed in pairs, and they typically stayed within a few hundred meters of each other throughout a drift. Figure 1 shows drifter tracks over 10-m bathymetry contours (bathymetry prepared by Akan et al. 2017). Measured wave heights ranged from 1 to 4 m, winds were typically 5–10 m s^{−1}, and drift speeds were up to 3.5 m s^{−1} on strong ebbs.

### a. Surface waves and wave breaking

Wave statistics were estimated with velocity data collected at 5 Hz with Qstarz BT-Q1000eX GPS loggers mounted in the buoy center at the water line. Following the Herbers et al. (2012) method, horizontal velocities are converted to sea surface elevation statistics using linear theory:

where *E*(*f*) is the wave spectrum, *S*_{uu}(*f*) and *S*_{υυ}(*f*) are the spectra of horizontal velocity components, *c*(*f*) is the wave celerity determined from dispersion, and *g* is acceleration due to gravity. Spectra were computed over short 5-min time series to better keep stationarity while the drifter moves rapidly through a heterogeneous environment. Significant wave heights were estimated as . The wave celerity *c* and wavenumber *k* are altered via wave current interaction and vertically sheared currents, and therefore the Kirby and Chen (1989) dispersion relation was used to estimate wavenumber via an iterative scheme. Wave breaking is observed on board each buoy using a GoPro Hero 2 camera. A full description of the wave-breaking statistics and wave spectral processing is available in Zippel and Thomson (2017).

### b. Raw turbulence data and motion correction

Velocities were measured using 2-MHz Nortek Aquadopp profilers. The Aquadopps were mounted inside the buoy spar, with the length of the Aquadopp body vertical. The Aquadopp heads were mounted in line with the body, such that the three acoustic beams were looking 20° off axis with vertical, and toward the surface. Samples were recorded at 4 Hz in pulse coherent burst mode, with a 5-min sampling interval (1024 samples per burst), in 16 profile bins spaced 4 cm apart with a 10-cm blanking distance. The profiler heads were mounted at 0.67 m depth such that the farthest bin from the Aquadopps was approximately at the ocean surface. Because the buoys were free drifting, the velocities measured in this reference frame were primarily turbulent fluctuations. The velocity range of the pulse coherent instruments was 1.15 m s^{−1} in the horizontal and 0.48 m s^{−1} in the vertical, allowing accurate measurement of turbulence in strong currents (drift speeds were measured over 3 m s^{−1}, but drifter slip relative to these currents was less than 10 cm s^{−1}). More details on the Aquadopp settings and sensitivity are in Thomson (2012).

The drifting platform primarily tracks with the free surface, such that velocity contamination by wave orbital motions is small. However, measured along-beam velocities **u**_{beam}(*z*, *t*) are contaminated by buoy motions, both translational (bobbing) and rotational (tilting) motions. We remove these motions from the time-domain-measured velocity as follows.

Pressure measurements on board the Aquadopp are converted to water depth using *z* = *P*/*ρg*, where *ρ* is the water density and *P* is the measured pressure. The vertical velocity relative to the free surface is then estimated using a centered difference of the measured water depth. The component of this velocity projected into along-beam coordinates is *u*_{bob,beam} = , where is the rotation matrix based on the measured heading, pitch, and roll, and is the position unit vector of the Aquadopp measurement bin relative to the beam transducer. Because the roll is measured in the yawed, pitched reference, and pitch is measured in the yawed reference, the transformation matrix is , with

where *h* is the heading, *p* is the pitch, and *r* is the roll.

The Aquadopps were mounted away from the center of motion of the buoy, looking off axis resulting in a nonorthogonal beam vector relative the rotational motions. Following Edson et al. (1998), the fixed-frame angular rate pseudovector **Ω** was estimated,

where represents a centered difference estimate of the derivative. The expected rotational velocity measured along a beam is then , where **m** is the center of motion of the buoy. There is a small amount of overlap, as the bobbing correction implicitly contains a component of the vertical rotational motion. We estimate the rotational effects to be relatively small when compared with the bobbing correction. Still, this may overcorrect motions in some cases and result in a bias in TKE.

Raw measured along-beam velocities were corrected in the time domain *u*(*t*)_{cor,beam} = *u*(*t*)_{meas,beam} − *u*(*t*)_{bob,beam} − *u*(*t*)_{rot,beam}. Because of the centered difference estimates for , , , and , the number of corrected samples in each burst became 1022, rather than the measured 1024. Before further processing, velocity data were quality controlled using the reported backscatter amplitude (*a* > 30 counts) and correlations (cor > 50; Elgar et al. 2001; Thomson 2012). Removed points were replaced with cubic interpolation (because a continuous series is needed for spectral analysis); however, if more than half of the 4-Hz samples were removed, the 5-min burst would not be used. Velocity spectra were estimated using Welch’s method, where each 1022-sample, motion-corrected time series was split into windows of 64 samples each with 50% overlap, and a Hamming taper was applied to each window. The FFTs of each window were then averaged, resulting in a power spectral density estimate with approximately 16 degrees of freedom.

Example velocity spectra are shown in Fig. 2. Typically, the vertical velocities due to bobbing (and estimated from the pressure measurement) accounted for most of the platform motion contamination, and the effects of rotational motion were relatively small. The translational bobbing motions were most apparent near 1 Hz, the estimated natural frequency of the buoy (Thomson 2012). A second peak associated with the tilting motions was also apparent at a slightly lower frequency. The bobbing motions contaminate the frequencies where the equilibrium range is typically observed. Once corrected for platform motions, however, a region with *f*^{−5/3} slope is evident in most spectra, consistent with an inertial subrange. Two of the Aquadopp beams faced away from the vane and were thus oriented away from any flow disturbance caused by platform (i.e., away from self-wake). The third beam, however, often experienced flow distortion and was therefore not used.

### c. TKE dissipation rates and TKE

Two methods were used to estimate TKE dissipation rates, the structure function method (Wiles et al. 2006) and a spectral method based on Tennekes (1975). We will briefly overview the former here, and a more complete description can be found in Thomson (2012). The second-order structure function is defined

where *b* is the spatial distance between measurements, denotes a time average, and *u*(*z*) is the demeaned along-beam velocity. The structure function can be related to the dissipation rate *ε* through

where *N* is an offset introduced through uniform noise across all measurements and is a constant (Wiles et al. 2006).

The spatial structure function has been the preferred TKE dissipation estimate in previous studies using SWIFT drifters because it is robust to platform motion. In this present study, motion correction is done directly, and thus the dissipation rate can also be estimated using the *f*^{−5/3} slope region of the spectrum. In the spectral method, velocity measurements in frequency are converted to turbulent wavenumber with an advective scale, *k*_{t} = 2*πf*/*U*_{adv}, and the assumption of a frozen turbulence field. This presents a problem in the free-drifting measurement frame as the relative velocity of the ambient water relative to the platform is near zero, *U*_{adv} ≈ 0. Here, we follow Tennekes (1975), where velocity spectra in the absence of an ambient current are expected to follow a self-advected form,

where *u*_{rms} is the root-mean-square of the demeaned velocity, *χ* = 8 is a constant (De Silva and Fernando 1994), and *ω* = 2*πf*. Equation (12) is inverted, and the mean slope in the inertial subrange over frequencies 0.68 < *f* < 1.5 Hz is used to solve for *ε*.

Example profiles of TKE dissipation rates are shown in Fig. 3. The methods agree favorably in magnitude across the two acoustic beams outside of the buoy wake. The spectral method shows more variation vertically. This may be due to increased spatial independence in estimating dissipation. That is to say, the structure function method uses distributed spatial information in estimating TKE dissipation rates, which may blur existing spatial gradients. The spectral method is more localized in space, with strict separation between estimates in depth, which may be the cause of increased vertical slope in Fig. 3. This is consistent with the work of Guerra and Thomson (2017), where structure function estimates also showed less spatial variation than the spectrally estimated dissipation rates (i.e., Guerra and Thomson (2017, their Figs. 6 and 12). A comparison of the spectral method and the structure function method across all measurements and depth bins is shown in Fig. 4.

TKE is estimated using the variance of the motion-corrected velocities along each beam. The variances from beams 2 and 3 (beam 1 is in the wake of the platform and thus avoided) are averaged, such that , assuming isotropy. In the case of nonisotropic turbulence, horizontal eddies larger than the drifting platform are not expected to be measured, as they would result in platform motion, not velocity fluctuations relative to the free-drifting platform. It is not clear if these eddies will retain importance when the measured turbulence analysis is in the free-drifting reference frame.

The velocity measurements in this study are referenced to the free-water surface ( coordinates) rather than the mean sea surface (*z* coordinates). While the balance of diffusion and dissipation has, to this point, been referenced to the *z*-coordinate reference, we will nonetheless test these equations from the surface-following reference. We note that the same set of equations can be derived from a diffusive–dissipative balance in substituting for *z*, if no extra terms arise in the TKE balance from the coordinate transform. In fact, the surface boundary condition of TKE flux is more physically justified at the free surface than the mean surface. The numeric implications of this reference frame mapping are discussed in section 4a.

## 3. Results

### a. Wave steepness

Zippel and Thomson (2017) showed that wave steepness is a strong indicator of wave breaking at river inlets, and that the relevant steepness is between the deep-water formula for whitecaps and the shallow-water formula for surf. The turbulence results suggest that this wave breaking is the dominant source of near-surface turbulence throughout the Columbia River mouth. Figure 5 shows a strong correlation of depth-averaged TKE and TKE dissipation rates with wave steepness. The strong correlation of turbulence values with wave steepness holds for estimates from both the spectral method and the structure function method. Appendix A evaluates non-wave-breaking sources of turbulence, including shear production, buoyancy, surface convergence, and bottom stress; the conclusion is that wave breaking is the dominant forcing for turbulence in the upper 0.5 m.

### b. Turbulent length scales

Figure 6 shows estimates from the measurement bin closest to the free surface (within 0.04 m of the instantaneous surface). The roughness lengths are related to both significant wave height and mean wavenumber. The correlation with wave height is strong, while the correlation with mean wavenumber is weaker (though still significant). Of course, wave height and wavenumber are not completely independent; waves are limited in steepness, such that *Hk*/tanh(*kd*) is roughly constant. Wave height estimates were more robust than mean wavenumbers (which assumed a dispersion relation for waves over a sheared current; see Zippel and Thomson 2017) and therefore may explain the difference in correlations.

For roughness length proportional to the wave height, , the data yield the parameter constraint . This is consistent with typical model values *L* = *κ* = 0.4, , and , but is underdetermined. Still, the results strongly support the wave height roughness length suggested in Terray et al. (1996) and moderately support the wavenumber roughness length suggested in Drennan et al. (1996) and in Jones and Monismith (2008). The results do not support a constant roughness length, or more complicated variations such as Gemmrich and Farmer (1999). The ratio *q*^{3}/*ε* has no correlation with drift speed *U* or drift speed normalized by the mean wave phase speed *U*/*C*_{p}. Thus, although the steepness leading to wave breaking may be controlled by the currents at the river mouth, the resulting surface turbulence appears to be unaffected by the currents.

Estimates of length scale just beneath the surface are consistent with the surface values, but do not further constrain model parameters. Figure 7 shows estimated length scales across all measurement depths plotted against the expected relationship with depth below the free surface. The more standard constant values and (Fig. 7a) and constant values found best fitting the dataset (see appendix B for more) and (Fig. 7b) fit the lognormal mean of the results well. Since measurements have a limited depth range, most of the variation in the data is explained by wave height, which ranges from approximately 0.25 < *H*_{s} < 3.5, where measured depths are limited to . This is highlighted in Fig. 7c, which shows that wave heights without depth variations are still correlated with turbulence results.

### c. Decay scales

Decay scales *λ* are estimated following the power-law model [Eq. (7)]. Applying the full model requires estimating the magnitude surface input TKE flux *G*, which is poorly constrained and not possible to directly estimate from these measurements. Instead of specifying *G*, the TKE and TKE dissipation rates are normalized by their measured near-surface values *q*^{2}(0) and *ε*(0) (rather than scaling them). The normalized values are expected to have the same decay scale, but may contain errors due to uncertainties in vertical placement *z*. For *ε*, this offset error is [1 + (Δ*z*/*z*_{0})]^{λ}, where Δ*z* is the error in position. For Δ*z* on the order of an Aquadopp bin size (0.04 m), we expect this vertical error to be small. Errors in normalized TKE and TKE dissipation rates are expected to be log distributed. The noise floor in such a normalization is therefore relative to the surface value. However, given surface values of *ε* ~ 10^{−3} m^{2} s^{−3}, and noise floors of *ε* ~ 10^{−5} m^{2} s^{−3}, we expect approximately two decades, with error increasing farther from the surface. Noise floor estimates of TKE are more complicated because of motion correction, but a similar range could be expected.

Normalized TKE and scaled depth (Fig. 8) show good agreement with open-ocean power-law decay models, albeit with very specific parameters. The more standard model parameters , (dotted line) do not fit the measurements well and tend to overpredict the amount of turbulence at depth. Decreasing the roughness length to and *L* = 0.3 (dashed line) gives better agreement, but is no longer fully consistent with the measurements shown in Figs. 6 and 7. The parameter values *L* = 0.21, , and (dashed dotted) do a reasonable job in matching the measurements and hold with the relation from Fig. 7b. The exponential solution also seems to fit the data well near the surface. However, the constants in the exponential solution are, again, relatively unconstrained. Taking from Fig. 7c, and constraining 0 < *L* < 0.4, we arbitrarily set as the approximate location where the power-law decay seems to dramatically lessen. A reasonable fit is achieved with , *L* = 0.135.

Normalized TKE dissipation rates are shown against scaled depths in Fig. 9. Here, parameters *L* = *κ*, , and predict more dissipation at depth than is measured (similar to TKE in Fig. 8). The other sets of constants all do reasonably well to predict the measured dissipation rate profiles. Constants *L* = 0.21, , and fit the TKE dissipation decay, TKE decay, and the length-scale arguments (Figs. 7b, 8). This set of constants gives *λ* = 3.6, which is larger than the decay reported in Terray et al. (1996). However, the scaling still represents the data well for , as shown by the red dashed–dotted line in Fig. 9. Therefore, the decay *λ* = 3.6 is still consistent with previous studies over the measured depth range. This highlights the importance of the roughness length in the argument of the power-law model, especially close to the surface. Estimates of normalized TKE dissipation from the structure function method are similar to those from the spectral method, but do show slightly reduced decays. Last, the exponential solution fits the top of the profile well, below which the power-law relation appears to hold.

### d. Fronts

Fronts are common in the vicinity of the river mouth, and these complicate the relation of wave steepness to turbulence shown in Fig. 5. The fronts are associated with strong horizontal gradients in currents −*dU*/*dx*. Such gradients are difficult to quantify with free-drifting buoys alone, because the drifters rapidly move toward the convergence zone and stay within it, providing no current or wave information on either side. Thus, the buoy estimates are highly localized within these gradients. Furthermore, these gradients may cause refractive focusing of the incident waves, which might cause increased wave breaking.

The dataset includes a few cases where spatial information is available in the form of surface current maps derived from airborne interferometric synthetic aperture radar (SAR) (Farquharson et al. 2014). An example of buoy 5-min drift positions are shown overlaid on a SAR composite velocity map in Fig. 10a, which is used as a case study of the interaction between the wave-breaking turbulence and the river plume.

The SAR velocity field is a composite of six aircraft passes over the estuary. Each pass took approximately 7 min to complete, so the surface velocity field shown in Fig. 10a evolved over a period of 42 min. This evolution, combined with calibration errors from pass to pass (<10 cm s^{−1}), accounts for some of the variation seen in data collected during each track in the composite field. Other pass-to-pass differences may be ascribed to a surface velocity measurement bias that depends on the SAR viewing geometry in areas of large subresolution (meter scale) waves (Thompson and Jensen 1993). This bias has not been characterized for this dataset because of the lack a comprehensive measurement of the surface wave field (including breaking) throughout the domain. Areas of noisy measurements are due to low backscattered signal from the surface.

Despite these sources of noise, a front can be seen by the rapid spatial change in velocities in the SAR-derived velocity field (approximate longitude −124.02; Fig. 10a). A large gradient in wave energy flux would be expected across such a current gradient −*dU*/*dx*, even on following currents, because of the rapid change in wave steepness required from conservation of wave action and dispersion (Chawla and Kirby 2002). The SWIFT buoys were visually confirmed to be caught in this front at approximate longitude −124.04, which matches with tower-based radar measurements of the front (not shown; see Honegger et al. 2017). The drifters stayed in the front until they were reset to avoid grounding, and therefore the drift track is discontinuous at longitude −124. SAR measurements lagged the timing of buoys in the convergence zone, and thus the front is not as visible at this leading edge in Fig. 10a. Measurements of wave-breaking rate (Fig. 10b) and TKE dissipation rates (Fig. 10c) increased in tandem where the horizontal gradient in currents was largest (although the shown SAR velocities lagged the buoy timing, and the front edge is offset). The increase in wave-breaking turbulence at the front is consistent with the results of Thomson et al. (2014), where a similar example is provided from an ebbing front offshore (using a different case from this same dataset). Unfortunately, a direct estimate the wave surface flux *G* used in Eq. (7) cannot be made for any of these cases, because the wave measurements are at the gradient (rather than across it). Still, the elevated turbulence associated with the waves can be compared to other sources of turbulence in the river mouth.

#### Relation to estuarine turbulence

Interfacial turbulence across plume fronts is typically attributed to shear production, which is opposed by stable buoyancy at a strong density interface. This shear-driven turbulence has been measured at similar levels to wave-breaking levels (*ε* in the range of 10^{−4} to 10^{−3} m^{2} s^{−3}; see Kilcher and Nash 2010; Horner-Devine et al. 2015). Therefore, it is worth investigating whether the surface turbulence *O*(10^{−3}) at the river mouth is not exclusively from wave breaking and is affected by shear and buoyant effects that are expected in the absence of waves. In a similar manner, there are questions on the effect of wave-breaking turbulence on these estuarine forced (shear production and buoyancy) terms.

For the front case already shown, CTD casts were made on either side of the velocity gradient using a YSI CastAway (cast locations shown in Fig. 10a). Vertical shear was estimated over the buoy’s measurement range near the surface. Profiles of density (Fig. 11a) are used to calculate stratification *N*^{2} and are shown along with the estimate of squared shear *S*^{2} in Fig. 11b. In addition, the second vertical derivative of TKE is shown as a proxy for diffusive transport. Squared shear and buoyancy frequency measurements are of a similar magnitude to values reported in Jurisa et al. (2016) and are much smaller than the wave-breaking proxy of diffusive transport near the surface. Figure 11c shows mean measured TKE dissipation rates and estimates for mean downward advective transport . The uncertainty in the advective transport is large; it may be a leading-order term very near the surface. Thus, downwelling at the front may be the most significant mechanism for the interaction of wave-breaking turbulence at the surface and estuarine turbulence at depth.

An estimate for the influence of buoyancy is made using the Ozmidov length scales in Fig. 11d. The Ozmidov length scale is defined and represents the smallest turbulent length scales affected by buoyancy. Profiles of the near-surface dissipation rate in the front (Fig. 11d) are extrapolated to the depths of measured density and used to estimate *L*_{O}, which is compared with the mixing length scale that best fits the turbulence measurements for the majority of the dataset [i.e., ]. These constants represent a solution for the cascading relation (Fig. 7) and are not limited to a diffusive–dissipative balance. Near the surface, mixing lengths are smaller than Ozmidov lengths, suggesting buoyancy is not affecting the turbulence at relevant scales.

The example presented here shows a single case where vertical transport (diffusive and/or advective) is large compared to buoyancy and shear near the surface, in a layer that is thin compared with plume thickness (typically ~10 m at the Columbia River; see Kilcher and Nash 2010). Therefore, it is likely that the near-surface turbulence examined herein is generally unaffected by buoyancy and shear production. Since mean downward velocities are only expected to be large when drifters are trapped in fronts (a small subset of the data), the diffusion–dissipation balance shown in the previous sections is expected to be valid. The data presented in this case study also suggest a region below the measurements ( cm) where wave-sourced turbulence (transport) and river effects (buoyancy and shear production) overlap and interact.

## 4. Discussion

### a. Model parameters

Specification of roughness length is clearly important in characterizing turbulence near the ocean surface, but it has been left relatively open as a model parameter in the literature. The direct physical implication of Eq. (2) is that is the length scale of turbulence introduced by the wave-breaking events. This length can be thought of either as that imposed by the front leading bubble plume (Longuet-Higgins and Turner 1974) or the spatial extent of the vortex tube created at the horizontal edges of the breaker (Clark et al. 2012). In the former case, the physical justification for can be seen through the work of Duncan (1981), where the length of the aerated fluid on a breaking wave is related to wave height. This physical justification, under the same measurements in Duncan (1981), could also be used to justify the scales and (where *C*_{p} is wave phase speed). In the latter case, Pizzo and Melville (2013) showed that the circulation under breaking waves can be scaled with both wave steepness and wave phase speed. Under this framework, the turbulent length scales in a Reynolds-averaged Navier–Stokes (RANS) model may be physically justified in a similar sense with wave parameters *H*_{s}, *k*^{−1}, and *C*_{p}. It is possible this study only finds better agreement between roughness length and *H*_{s} because SWIFT wave measurements more accurately quantify *H*_{s} than *k* and *C*_{p}. The physical justifications for scaling the vertical coordinate by *k*^{−1} or , therefore, may be equally valid to *H*_{s}.

As discussed in Burchard (2001) and Gerbi et al. (2009), setting *L* < *κ* suggests the turbulent length scale grows more slowly in the dynamic surface boundary layer than in rigid boundary layers. Gerbi et al. (2009) evaluated the ratio , but did not distinguish between the two parameters, finding the ratio was smaller than expected for a rigid boundary layer, . However, as shown in Figs. 7a and 7b, this ratio is also sensitive to choice of roughness length. Therefore, the much larger is only justified here through the smaller choice of . The measurements reported here suggest a different value for the stability function than typically used in turbulence models. However, many of these models assume turbulent equilibrium (that *P* + *B* = *ε*) in determining stability functions, which may not be reasonable near the ocean surface. Work to include the dynamic effects of waves on turbulence closure has been started recently (Harcourt 2013), but these have yet to be fully incorporated and tested within turbulence models.

The success of the Terray et al. (1996) scaling below and the exponential solution for the top of the water column (Fig. 9) leads naturally to a piecewise scaling for dissipation referenced to coordinates,

Appropriate choice of *α* could then rectify the range of decay scales (1 < *λ* < 2) found very near the surface (e.g., Gemmrich 2010; Sutherland and Melville 2015; Thomson et al. 2016). The constant *A* must be a function of decay scale such that the piecewise function is consistent at the interface .

### b. Reference frames

Some discrepancy in reported slope *λ* may be attributed to choice of reference frames. As reported in Thomson et al. (2016), turbulence lasting longer than one wave period is moved vertically with the free surface, and thus fixed frame measurements capture an effective average of the turbulence at depths varying from the free surface. In other words, the advection of a depth-varying turbulent field across fixed frame instruments creates a complicated mapping between *z*-referenced coordinates and -referenced coordinates. Lumley and Terray (1983) investigated a similar effect, where the advected turbulence is assumed uniform across wave orbitals. Given that turbulence decay is appreciable over a wave height (Figs. 7–9), this assumption of vertically homogenous turbulence across a wave particle excursion is likely not valid near the surface, where particle excursions (below the free surface) are roughly equal to sea surface deviations. Here, we attempt a simple numeric estimation for the orbital advection of a nonuniform field. The key simplification is the assumption that the TKE dissipation rate is a 1D field with a robust average in time. A more realistic approach would evaluate the effect on the velocity field [as was done in Lumley and Terray (1983) and Rosman and Gerbi (2017)], rather than directly on the higher statistical moments of the velocity field (i.e., the TKE dissipation rate). Such a description would require knowledge of how the velocity field statistics vary with depth a priori.

Following Andrews and McIntyre (1978), we take a Taylor series expansion of the -coordinate dissipation function about the coordinate *ζ*. Here, we approximate the wave orbital particle excursions by setting *ζ* = *η*, which is valid for long waves relative to the depth range *z*. In this way, we can relate the two reference frames as a function of coordinate and sea surface statistics. The first five terms in the expansion are

where is the normalized dissipation function; is a time average operator; is the surface referenced vertical coordinate, such that ; and ′ denotes differentiation with respect to *z*. Assuming *η* is Gaussian distributed, the first and third statistical moments are zero,

such that only the sea surface variance and kurtosis are needed. In this way, one can relate the -coordinate and *z*-coordinate mean dissipation rates [i.e., to ]. Assuming a power-law form, the nondimensional dissipation rate , the second and fourth derivatives are

This representation does not account for the partial drying of *z*-coordinate measurements, which are exposed to air in wave troughs. Given that represents a rate, the nonwetted time must be accounted for in the average. Here, a more general, probabilistic solution accounts for this issue,

where *P*(*η*) is the distribution of sea surface elevations. Unfortunately, this integral is not easily solved analytically for many probability distribution functions. A numerical solution Eq. (18) assuming Gaussian *P*(*η*) is shown in Fig. 12, plotted alongside the Taylor expansion approximate solution and the power law referenced to coordinates with constants *L* = 0.22, , and *z*_{b} = 0.32. Choice of *η* over the wave orbital particle excursions *ζ* simplifies much of the analysis, but likely adds more turbulence at depth in the *z*-coordinate reference since wave orbital excursions are necessarily smaller than the surface displacements, *ζ* ≤ *η*. The assumption of Gaussian distributed sea surface elevations does not describe nonlinearities common in wave fields, but is often a reasonable approximation (see, e.g., Schwendeman and Thomson 2017, their Fig. 6).

A few notable features appear in the *z*-coordinate numeric mapping shown in Fig. 12. First, on a logarithmic scale (Fig. 12a), the shape of the mapping is qualitatively similar to the Terray et al. (1996) model. In particular, there is a nearly constant layer of dissipation above a power-law decay region due to the fractional coverage of water in the region, and thus it must be accompanied by a net TKE dissipation above the mean sea surface. Integration of the fixed-frame dissipation profile shows 30% of the total dissipation exists above the mean sea surface and 50% exists above *z*/*H*_{s} = 0.6 (the assumed breaking layer depth in Terray et al. 1996). The conversion also suggests measurements of *λ* in the fixed reference frame could decay faster than those in the moving reference frame below the breaking-layer depth. Because of the number of simplifications used, the numeric mapping presented here is only intended as a rough estimate. It does, at least, provide results qualitatively similar to the direct coordinate mapping of Thomson et al. (2016), wherein the surface-following estimates are maximum, on average, at the mean sea level [i.e., similarities between Fig. 12 herein and Figs. 2f,h in Thomson et al. (2016)].

## 5. Summary

Measurements of waves and near-surface turbulence at the mouth of the Columbia River were compared with existing analytical models and turbulence scalings. The observed surface turbulence is consistent with wave-breaking parameterizations developed for the open ocean, despite the uniqueness of the wave-breaking mechanism at the river mouth. This may be related to the relatively thick river plume, which is generally at least twice as thick as the waves are tall. Thus, stratification and shear are not relatively strong at the surface where waves are breaking.

The vertical dependence of the surface turbulence is consistent with a classic analytic model balancing diffusion and dissipation. TKE dissipation rates also follow the canonical *λ* = 2 power-law decay with depth, but for a range of scaled depths different than those originally proposed. Further, turbulent length scales estimated from measurements are seen to increase linearly with depth, supporting the a priori assumption.

The model is sensitive to the choice of constants, primarily the roughness length. Measurements suggest this roughness length is proportional to the significant wave height, perhaps because of advection by wave orbital motions. We find that the method used in processing turbulence data moderately changes the result; the spatial structure function blurs the vertical gradient, relative to a frequency spectrum approach, and this yields slightly decreased decay scales.

A mapping of coordinates from a reference frame moving with the free surface to a reference frame fixed at the still-water level is discussed. This mapping is dependent on the sea surface statistics and may help explain some of the discrepancies in reported power-law decay exponents.

## Acknowledgments

Funding for this project was provided by the Office of Naval Research as part of the RIVET-II DRI, and for the DARLA group. SWIFT data are available at http://www.apl.washington.edu/SWIFT. Cigdem Akan prepared the bathymetry data. Hans Burchard, Johannes Gemmrich, Eric D’Asaro, and two anonymous reviewers provided valuable feedback. Thank you to the captains and crews of the R/V *Oceanus*, the R/V *Pt. Sur*, and the F/V *Westwind*, and the field engineers Alex de Klerk and Joe Talbert for invaluable help collecting data.

### APPENDIX A

#### Non-Wave-Breaking Turbulence Effects

In river plumes outside of the near-surface layer, strong vertical shear is the mechanism for creating turbulence in stably stratified buoyant layers through Kelvin–Helmholtz instability. This leads to a three-term balance between TKE dissipation rate, TKE production from shear, and stabilizing buoyancy, *ε* = *P* + *B*, where production and buoyancy are often modeled and , respectively. Here, *N*^{2} is the squared buoyancy frequency, *S*^{2} is the squared shear frequency, *ν*_{t} is the eddy viscosity, and is diffusivity of buoyancy. The size of shear production and buoyancy are estimated below.

If production and buoyant terms are included in the wave-breaking TKE equation, Eq. (1) becomes

Direct comparison of buoyancy, shear production, and transport divergence terms is made difficult by estimation of turbulent diffusivities , and and eddy viscosity . Here, we expand the diffusive transport using the product rule, such that . A component of the vertical transport term can be compared to the shear frequency and the buoyancy frequency without estimation of eddy viscosity and diffusivity and are only modified by the relevant Schmidt numbers. Since only represents a component of the total transport and is expected to be the same sign as , the ratios and are likely overestimating the importance of shear and buoyancy relative to vertical diffusive transport.

Shear number *S*^{2} can be estimated from velocity data collected on board SWIFT buoys. The along-beam, motion-corrected velocities are rotated into east–north–up coordinates using the heading pitch and roll data in the time domain and averaged over each 5-min burst period. The average east and north velocities are smoothed vertically using a moving-average filter, and then shear is estimated using a vertical centered difference scheme. These shear profiles are not corrected for Stokes drift effects, which may bias the shear in the lower water column measurements high. TKE profiles are smoothed in a similar manner, and is calculated using a second-order central difference scheme. A histogram of the ratio is shown in Fig. A1. Although the variance of the distribution is large, the majority of data show that *S*^{2} is less than 10% of , and the mean value of the ratio is . Given that this ratio overestimates the importance of shear due to Stokes effects and underestimates the importance of transport due to neglect of a term, Fig. A1 strongly suggests that shear production is small relative to the vertical diffusive transport very near the surface.

Estimates of buoyancy frequency *N*^{2} are unfortunately not available for the majority of the dataset. Jurisa et al. (2016) reports values of 0.009 < 〈*N*^{2}〉 < 0.02 from the River Influences on Shelf Ecosystems (RISE) project (Hickey et al. 2010). The Columbia River plume in particular has notable linear velocity shear and density profiles (Kilcher and Nash 2010), which result in constant vertical *N*^{2}. Extrapolating the larger 〈*N*^{2}〉 reported in Jurisa et al. (2016) to the surface gives a range very similar to the center of the squared shear distribution, *O*(10^{−2}) s^{−2}. Thus, we draw a similar conclusion, that buoyant effects are secondary to vertical diffusive transport in the near-surface layer.

Shear and buoyancy numbers have a secondary effect in modifying the stability function *C*_{μ}, which directly influences *ν*_{k}. For a roughly constant shear number, this would modify *R* in Eq. (6), adjusting the decay exponent. The ratio was calculated using the Canuto et al. (2001) stability function, but revealed no preferential sorting in decay scales of either TKE (Fig. 8) or TKE dissipation rate (Fig. 9). It is likely then that vertical shear effects are not biasing the decay scales estimates.

Bottom boundary layer turbulence is also estimated to be small relative to surface layer turbulence. The bottom boundary layer TKE dissipation rate is expected to scale as , where the subscript BBL refers to bottom boundary effects referenced from the seabed. For general values of and *z*_{BBL} = *d*, this scaling yields *ε*_{BBL} < 10^{−5} m^{2} s^{−3}, which is below the measurement noise floor of the SWIFT buoy. Shallow regions with fast flows are the most likely to exhibit bottom-generated turbulence near the surface; however, we also note that these regions are typically impeded by stabilizing buoyant plume effects. Data collected upriver were primarily on flood tides where velocities are smaller than on ebb tides.

Surface convergence at river plume fronts can create large downward velocities. These velocities can contribute to mean downgradient advection of turbulence, . Outside of convergence zones, measured vertical velocities were typically smaller than the measurement error, and vertical TKE gradients were on the order *O*(10^{−2}–10^{−1}) s^{−2} m s^{−2}, from which we estimate that the mean downgradient advection of TKE was typically ≤10^{−5} m^{2} s^{−3}, outside of the measurement range of the dissipation rates. However, convergence at plume fronts created mean downward velocities on the order of 10^{−2}–10^{−1} m s^{−1}, which resulted in downgradient mean advection being a leading-order term in the TKE equation.

### APPENDIX B

#### Estimates of Turbulence Constants

Values for constants , *L*, and were estimated by minimizing the sum of log residuals between measurements and models via Eqs. (2)–(5). Confidence intervals on these values were estimated using a smooth bootstrap method, where measurements of TKE and *ε* were resampled randomly, and best-fit values were recomputed. Resampled values of TKE and *ε* were adjusted by adding lognormal distributed noise with zero mean and standard deviation equal to 20% of the measured value. The measurements were resampled 1000 times, giving robust estimates for the mean and variance of , *L*, and . In addition, a second error estimate was made to roughly quantify the effects of TKE production with a magnitude on the order of 10% of dissipation rate. One-sided noise (additive only) was added in the decay-scale fits [Eq. (4)] during resampling. Mean best-fit values with confidence intervals of twice the standard deviation were , , and for the initial bootstrap and , , and for the bootstrap where ad hoc production bias was simulated.

## REFERENCES

*Proc. IEEE Int. Geoscience and Remote Sensing Symp.*, Quebec City, QC, Canada, IEEE, 2661–2664, https://doi.org/10.1109/IGARSS.2014.6947021.

## Footnotes

^{a}

Current affiliation: Woods Hole Oceanographic Institution, Woods Hole, Massachusetts.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).