An improved method for determining near-surface currents from wave dispersion measurements

TL;DR

This paper introduces a simple, parameter-free inversion method that significantly improves the accuracy of near-surface current profile estimation from wave spectra, validated by laboratory experiments with real-world relevance.

Contribution

The paper presents a novel, easy-to-implement inversion technique that enhances current profile estimates from wave data without requiring user input or complex parameters.

Findings

Achieved over 3x accuracy improvement in laboratory tests.

Validated method's effectiveness on realistic ocean-like flows.

Demonstrated method's simplicity and minimal additional effort.

Abstract

A new inversion method for determining near-surface shear currents from a measured wave spectrum is introduced. The method is straightforward to implement and starts from the existing state-of-the-art technique of assigning effective depths to measured wavenumber-dependent Doppler shift velocities. A polynomial fit is performed, with the coefficients scaled based on a simple derived relation to produce a current profile that is an improved estimate of the true profile. The method involves no user-input parameters, with the optimal parameters involved in the polynomial fit being chosen based on a simple criterion involving the measured Doppler shift data only. The method is tested on experimental data obtained from a laboratory where current profiles of variable depth dependence could be created and measured by particle image velocimetry, which served as "truth" measurements. Applying the new inversion method to experimentally measured Doppler shifts resulted in a $>3\times$ improvement in accuracy relative to the state-of-the-art for current profiles with significant near-surface curvature. The experiments are dynamically similar to typical oceanographic flows such as wind-drift profiles and our laboratory thus makes a suitable and eminently useful scale model of the real-life setting. Our results show that the new method can achieve improved accuracy in reconstructing near-surface shear profiles from wave measurements by a simple extension of methods which are currently in use, incurring little extra complexity and effort. A novel adaptation of the normalized scalar product method has been implemented, able to extract Doppler shift velocities as a function of wavenumber from the measured wave spectrum.

Tables

Table 1. Table 1: Summary of properties for the four laboratory current profiles.

Profile	Flow Type	Water Depth	Flow rate	$\Gamma$
		[mm]	$[m^3/$s$]$	$\times {10}^{-5}$ $\left[m^3{s}^{-2}\right]$
a	Stagnation region	95	0.021	$3.8\pm 0.05$
b	Stagnation region	95	0.017	$2.8\pm 0.1$
c	Stagnation region	95	0.014	$2.7\pm 0.1$
d	Curved mesh	80	0.014	$6.7\pm 0.1$

Table 2. Table 2: Summary of the optimal parameters and results of the PEDM applied to x 𝑥 x -components of experimentally measured Doppler shifts.

Profile	$n_{\max }$	$\delta z_T$	$\delta z_B$	$\Delta U_{\mathrm{RMS}}^{\mathrm{EDM}}$	$\Delta U_{\mathrm{RMS}}^{\mathrm{PEDM}}$	$\frac{\Delta U_{\mathrm{RMS}}^{\mathrm{EDM}}}{\Delta U_{\mathrm{RMS}}^{\mathrm{PEDM}}}$
		[mm]	[mm]	[mm/s]	[mm/s]
a	8	0.5	17.9	34.2	8.9	3.8
b	8	0.5	3.1	21.9	4.3	5.1
c	10	1.7	5.2	14.2	3.0	4.8
d	3	4.0	7.3	3.0	3.4	0.9

Equations

$${\omega }_{\mathrm{DR}}(\mathbf{k})={\omega }_0(k)+\mathbf{k}\cdot \tilde{\mathbf{c}}(k),$$

$$\tilde{\mathbf{c}}(k)=2k{\int }_{-\infty }^0\mathbf{U}(z){\mathrm{e}}^{2kz}dz,$$

$${\tilde{c}}_i=-\frac{S}{2k_i}+U_0=U\left(z=-(2k_i{)}^{-1}\right).$$

$$Z_{\mathrm{eff},\mathrm{lin}}(k)=-(2k{)}^{-1}.$$

$$\tilde{c}(k)=\sum _{n=0}^{\infty }n!u_n{\left(-\frac{1}{2k}\right)}^n$$

$$U_{\mathrm{EDM}}(z)=\sum _{n=0}^{\infty }n!u_n{z}^n.$$

$$U_{\mathrm{EDM}}(z)\approx \sum _{n=0}^{n_{\max }}{u}_{\text{EDM},n}{z}^n,$$

$$U_{\mathrm{PEDM}}(z)=\sum _{n=0}^{n_{\max }}\frac{1}{n!}{u}_{\text{EDM},n}{z}^n.$$

$${ϵ}_{\mathrm{RMS}}=\sqrt{\overline{({\tilde{c}}_i-{\tilde{c}}_{F,i}{)}^2}},$$

$$P^l(k_x,k_y,\omega )=|P_x^l(k_x,k_y,\omega ){|}^2+|P_y^l(k_x,k_y,\omega ){|}^2,$$

$$F(x,y,t)=\exp \left[-\frac{1}{2{\sigma }_m^2}\left(\frac{x^2}{L_x^2}+\frac{y^2}{L_y^2}+\frac{t^2}{T^2}\right)\right],$$

$${\omega }_0(k)=\sqrt{\left(gk+\frac{\sigma }{\rho }{k}^3\right)\tanh kh},$$

$$S(k)={\int }_0^{2\pi }{\int }_0^{\infty }\mathrm{d}\omega \mathrm{d}\theta \left(|P_x(\mathbf{k},\omega ){|}^2+|P_y(\mathbf{k},\omega ){|}^2\right),$$

$$G_i(\omega ,\theta ,{\tilde{c}}_{x,i},{\tilde{c}}_{y,i})=\exp \left[\frac{(\omega -{\omega }_{\mathrm{DR}}(k_i,\theta ){)}^2}{4a(\theta )}\right],$$

$$P_i^{\prime }(\theta ,\omega )=10\log P(k_i\cos \theta ,k_i\sin \theta ,\omega )+10\log P(2k_i\cos \theta ,2k_i\sin \theta ,2\omega ),$$

$$N({\tilde{c}}_{x,i},{\tilde{c}}_{y,i})=\frac{⟨G(\omega ,\theta ,{\tilde{c}}_{x,i},{\tilde{c}}_{y,i})P_i^{\prime }(\theta ,\omega )⟩}{⟨G⟩⟨P^{\prime }⟩},$$

$$k\Delta c_{\Delta k}(k)=\frac{\mathrm{d}{\omega }_0(k)}{\mathrm{d}k}\Delta k.$$

$$\Delta c_{\Delta \omega }(k)=\Delta \omega /k,$$

$${\mathrm{Fr}}_S(k)=\frac{1}{\sqrt{gk}}{\int }_{-\infty }^0\mathrm{d}z\mathbf{k}\cdot {\mathbf{U}}^{\prime }(z){\mathrm{e}}^{2kz}=\frac{⟨S{⟩}_k}{{\omega }_0(k)},$$

Full text

\draftfalse\journalname

JGR:Oceans

Department of Energy and Process Engineering, Norwegian University of Science and Technology, Trondheim, Norway

\correspondingauthor

Benjamin K. Smeltzerbenjamin.smeltzer@ntnu.no

{keypoints}

Simple and more accurate method for reconstructing near-surface current profiles from wave spectra

Laboratory setup where shear currents and waves can be well-controlled and measured

Novel adaptation of a method for extracting wavelength-dependent Doppler shifts from wave spectra

1 Introduction

Characterizing near-surface ocean currents is of importance to a vast range of applications. At a fundamental scientific level, near-surface currents influence the exchange of energy and momentum between the air and sea (Terray et al., 1996; Kudryavtsev et al., 2008), impacting climate models. At a more practical level, currents affect wave-body forces, and can be relevant for operational safety in coastal areas Dalrymple (1973); Zippel & Thomson (2017). Accurate measurements of the mean flow in the top meters of the water column are difficult to obtain, in large part due to the presence of waves which induce platform motions and additional sources of noise. Conventional methods such as acoustic Doppler current profiling (ADCP) typically discard data in the topmost few meters of the water column.

An attractive alternative to in situ techniques is to deduce currents from measurements of waves, whose dispersion is altered by the presence of a background flow. The approach has the advantage of enabling remote sensing methods such as radar or optical-based detection, with the potential for mapping currents over a larger area (multiple $\mathrm{km}^{2}$) compared with point measurements. In addition, waves are most sensitive to currents near the free surface, precisely the regime where other conventional methods such as ADCP struggle. The vast majority of wave-based near-surface current measurements reported in the literature have used radar, including high frequency (HF) radar (e.g., Crombie, 1955; Young & Rosenthal, 1985; Stewart & Joy, 1974; Ha, 1979; Fernandez et al., 1996; Teague et al., 2001; Shrira et al., 2001) and more recently X-band radar systems (e.g., Young & Rosenthal, 1985; Gangeskar, 2002; Lund et al., 2015; Campana et al., 2016, 2017; Lund et al., 2018), also in some cases to reconstruct the bathymetry (e.g., Hessner et al., 2014; Hessner & Bell, 2009). Optical methods have also been used to a lesser extent (Dugan et al., 2001; Dugan & Piotrowski, 2003; Laxague et al., 2017; Horstmann et al., 2017; Laxague et al., 2018).

Though wave-based current measurements offer several distinct advantages compared to other methods, they have a number of inherent challenges. Firstly, determining the current profile as a function of depth without stringent a priori assumptions as to the functional form requires the ability to measure waves over a spectrum of wavelengths and directions. The quality of results is thus dependent on the sea state (Lund et al., 2015; Campana et al., 2016, 2017). Secondly and more fundamentally, the determination of the current depth profile from wave dispersion measurements is a mathematically ill-posed inverse problem. The inferred current profile is not mathematically unique, and noise in wave measurements gets amplified in the inversion process (Ha, 1979). As a result, a priori assumptions and constraints of the depth-dependence of the current profile based on physical intuition have typically been imposed.

Despite these difficulties, wave-based current measurements have been used in the field for many decades. The techniques involve reconstructing the near-surface current from measured alterations to the wave frequency, and are termed “inversion methods.” The most common and elementary methods involve determining a single current vector representative of a weighted average of the near surface flow, with other more recent methods reconstructing some degree of detail as to the depth-dependence of the flow. In reviewing the previously developed inversion methods, we first consider the dispersion relation for small-amplitude linear waves propagating atop a depth-varying flow, which can be approximated as:

[TABLE]

where $\omega_{\mathrm{DR}}$ is the wave frequency, $\omega_{0}$ the frequency in quiescent waters, $\mathbf{k}=\{k_{x},k_{y}\}$ the wavevector, $k=|\mathbf{k}|$, and $\tilde{\mathbf{c}}$ a wavenumber-dependent Doppler shift velocity due to the background current. The $z=0$ plane is the undisturbed water surface and the bottom is found at $z=-h$ with $h>0$. We shall mostly work in the deep water regime $kh\gtrsim\pi$ where $h=\infty$ can be assumed. As first shown by Stewart & Joy (1974), the Doppler shift can be approximated as a weighted average of the current profile as a function of depth as

[TABLE]

where $\mathbf{U}(z)=[U(z),V(z)]$ is the current profile. The finite depth version of the Stewart & Joy (SJ) approximation (2) was derived by Skop (1987) and extended by Kirby & Chen (1989). The weighting term decays exponentially with depth (in deep water), reaching a value of 0.2% of the surface value at a depth equal to half the wavelength ($kz=-\pi$). Short wavelengths are thus sensitive only to currents near the surface, whereas longer wavelengths are affected by currents at greater depths. The inversion method involves using values of $\tilde{\mathbf{c}}(k)$ obtained from experimental data to determine the unknown $\mathbf{U}(z)$.

A word of warning is warranted when referring to $\tilde{\mathbf{c}}$ as the “Doppler shift” as is conventional. While $\tilde{\mathbf{c}}$ occurs in (1) exactly as would a Doppler frequency shift resulting from Galileian transformation upon changing reference system, it should not be interpreted as such. A misunderstanding has arisen from this name that the same Doppler shift should also be added to the wave’s group velocity to account for the shear, but this is not correct as pointed out by Banihashemi et al. (2017). Rather, the group velocity remains $d\omega/dk$, for which taking the $k$-dependence of $\tilde{\mathbf{c}}$ into account is key. We shall follow the numenclatorial convention in the literature and refer to $\tilde{\mathbf{c}}$ as the Doppler shift velocity while bearing this in mind.

Various wave detection methods are sensitive to different spectral ranges of $k$ and have led to the development of a number of inversion methods. In the case of HF radar, the detected signal is dominated by resonant Bragg scattering, effectively measuring the Doppler velocity of a surface wave with a wavelength half that of the radar system. Data reported from single-frequency HF radar is often referred to as the surface current, yet more precisely it is a weighted average of the current profile from (2), as it essentially measures $\tilde{\mathbf{c}}(k_{\mathrm{HF}})$ ($k_{\mathrm{HF}}$ being the wavenumber of the resonant wave) without information concerning the depth-dependence. Depth-profile information can be obtained by using multiple radar frequencies (Stewart & Joy, 1974; Ha, 1979; Fernandez et al., 1996; Teague et al., 2001) which probe different resonant wavenumbers. Similarly, other detection methods such as X-band radar or optical techniques inherently measure a wide spectrum of wavelengths, thus evaluating (2) at many $k$-values and enabling the use of inversion methods to estimate the current depth-dependence.

Inversion methods of determining $\mathbf{U}(z)$ from a set of measured values of $\tilde{\mathbf{c}}_{i}=\{\tilde{c}_{x,i},\tilde{c}_{y,i}\}$ at discrete wavenumbers $k_{i}$ can be carried out separately for each velocity component, i.e. $U(z)$ can be found from $\tilde{c}_{x,i}$, and $V(z)$ from $\tilde{c}_{y,i}$. To ease the notation, in the following we outline the new inversion method using $U(z)$ and $\tilde{c}_{i}$ to denote the flow velocity and Doppler shift velocities, with the implicit understanding that they may correspond to either dimension in the horizontal plane. The subscript $i$ indicates that the respective variable takes on a discrete set of values as may be extracted from experimental data.

Assuming a given functional form to the current profile, one can assign effective depths to the measured Doppler velocities based on the wavenumber by finding the depth at which the Doppler velocity is equal to the current, i.e. $\tilde{c}_{i}=U{(}Z_{\mathrm{eff}}(k_{i}){)}$. For the commonly assumed case of a current profile which varies linearly with depth, $U(z)=Sz+U_{0}$, where $S$ is the vorticity and $U_{0}$ the surface current. By the approximation (2) the Doppler shift is approximated as

[TABLE]

The last form shows that assuming linear current, deep water and using the SJ approximation, the appropriate effective depth is

[TABLE]

(In other words $Z_{\mathrm{eff}}(k)$ is approximately $8\%$ of the wavelength.) A similar relation can also be derived for a logarithmic profile (Plant & Wright, 1980). We refer to the method of estimating $U(z)$ from a measured $\tilde{c}_{i}(k)$ using (3) or its sibling assuming a logarithmic profile as the effective depth method (EDM). The EDM has been used extensively in the literature for estimating near-surface shear currents (e.g., Teague et al., 2001; Lund et al., 2015; Laxague et al., 2017, 2018; Stewart & Joy, 1974; Fernandez et al., 1996).

A clear weakness of the EDM, however, is that it relies on assumptions as to the functional form of the depth dependence. Ha (1979) developed a method for inverting (2) directly based on a series of measured $\tilde{\mathbf{c}}$ values, which was further developed and applied to data from X-band radar by Campana et al. (2016). The method involves a Legendre quadrature approximation to the integral, with constraints on the curvature of the current profile as well as the distance from an initial guess in order to suppress the amplification of experimental noise. The method avoids initial assumptions as to form of the current profile and yields current estimates at greater depths. The reconstructed $U(z)$ has comparable accuracy relative to the EDM when compared against ADCP “truth” measurements.

We present a new inversion method which is completely free of parameters. The method, which is derived assuming deep water, uses the current profile obtained by the EDM, and fits it to a polynomial function. The method then makes use of a simple relation which follows from (2) to construct an improved estimate of the true profile $U(z)$ directly from the coefficients of the fit. The method is validated and tested on experimental data from a laboratory setup, where the background current velocity profile and wave spectrum could be well-controlled and characterized.

In the following we describe the new method in Section 2, and examine its performance also in finite water depth. Section 3 describes the experimental setup and analysis of the data, where an adapted version of a normalized scalar product (NSP) method is used to extract Doppler shifts from wave spectra. Section 4 demonstrates the use of the new inversion method on experimentally measured Doppler shifts, where in situ measurements of the current profile are used as “truth” measurements for validation. The performance of the method is evaluated by considering the fractional decrease in error of the depth profile achieved by the new inversion method compared to the EDM.

2 Polynomial effective depth method

From experimental data of the wave spectrum, a set of Doppler shift velocities $\tilde{\mathbf{c}}_{i}$ at unique wavevector magnitudes $k_{i}$ can be obtained by a number of methods such as least squares techniques (Campana et al., 2017; Senet et al., 2001), or a normalized scalar product (NSP) method (Huang et al., 2016; Huang & Gill, 2012; Serafino et al., 2010) used herein (described in section 3).

Assuming a polynomial current profile of the form $U(z)=\sum_{n=0}^{\infty}u_{n}z^{n}$ in deep water, evaluation of (2) yields the SJ approximation

[TABLE]

for the Doppler shift velocities. We notice that $(-2k)^{-1}$ is equal to the mapping function $Z_{\mathrm{eff,lin}}(k)$ used in the EDM assuming a linear frofile, equation (4). Using the EDM with this mapping the estimated current profile is

[TABLE]

Thus, the mapped profile $U_{\mathrm{EDM}}(z)$ is also of polynomial form with coefficients of the $n$-th order term differing by a factor $n!$ from those of the true profile $U(z)$. The estimated velocity profile $U_{\mathrm{EDM}}(z)$ will suffer from inaccuracies since the mapping function is not the correct one. The new inversion method, referred to hereafter as the polynomial effective depth method (PEDM), seeks to improve this by making use of the simple relationship between the coefficients in the series representation of $U_{\mathrm{EDM}}(z)$ and the true profile $U(z)$, namely that they differ by a factor $n!$.

Explicitly, the PEDM procedure consists of the following three steps:

For each of the measured values $\tilde{c}_{i}$, assign effective depths $z_{i}=-(2k_{i})^{-1}$ according to the EDM procedure of (3) and (4). 2. 2.

Obtain $U_{\mathrm{EDM}}(z)$ by fitting the set of points $\{z_{i},\tilde{c}_{i}\}$ to a polynomial of degree $n_{\mathrm{max}}$:

[TABLE]

where $u_{\text{EDM},n}$ are the coefficients obtained in the polynomial fit. 3. 3.

Then the improved PEDM estimate is

[TABLE]

Equation (8) follows immediately from a comparison of (6) and (7), where $u_{\text{EDM},n}=n!u_{n}$.

2.1 Theoretical limitations

Two notable potential complications arise: finite depth where (5) and (6) are no longer strictly valid, and realistic situations where errors in experimentally measured Doppler shifts, which in addition are measured at a finite range of wavenumbers, lead to errors in the fitted polynomial coefficients rapidly increasing for higher values of $n$.

2.1.1 Performance for finite depth

In the case of finite depth $h$, an explicit relation of the form of (6) cannot be derived since the mapping function $Z_{\mathrm{eff}}(k)=-(2k)^{-1}\tanh kh$ in finite depth cannot be solved with respect to $k$ analytically, but must be inverted numerically. The approximation (2), moreover, obtains a more complicated form less amenable to analytical treatment (Skop, 1987). To examine the effect of finite depth on the accuracy of the PEDM, we consider an exponential profile of the form $U(z)=U_{0}\exp(\alpha z)$, with \alpha$$=8\cdot\mathrm{min}[k_{i}]/\tanh(\mathrm{min}[k_{i}h]) to preserve the same functional form within the range of mapped depths regardless of the water depth, and $U_{0}$ the surface current.

We consider the implementation of the PEDM in finite depth with $n_{\mathrm{max}}=10$, simply using the finite depth mapping function in step 1 of section 2 to assign effective depths $z_{i}=-(2k_{i})^{-1}\tanh k_{i}h$. Steps 2-3 of the PEDM procedure were unchanged. The fractional depth-integrated root mean square (RMS) error between $U_{\mathrm{PEDM}}(z)$ and $U(z)$ was calculated for cases over a range of water depth values $\mathrm{min}[k_{i}]h$, with the results shown in Figure 1a. For all but the shallowest depths considered here, the deep water mapping function results in errors at the 1% level. For most realistic combinations of water depth and relevant wavenumbers, Figure 1a indicates that the finite depth mapping function and (6) yield sufficient accuracy.

2.1.2 Effect of limitations of measured Doppler shifts

As mentioned, the fact that $\tilde{c}_{i}(k)$ is measured for a finite range of wavenumbers will affect accuracy. This is true of any inversion method for reconstructing $U(z)$ from dispersion measurements.

To handle the realistic case of experimentally measured Doppler shifts at a finite range of wavenumbers, we extend the three-step process described in section 2 to a 6-step process (the first three steps of which are illustrated schematically in Figure 1b):

Fit the mapped Doppler shifts to a polynomial of order $n_{\mathrm{max}}$ to produce the profile $U_{\mathrm{EDM}}^{\mathrm{init}}(z)$ (Steps 1-2 in section 2), using the finite depth mapping function if appropriate. 2. 2.

Create additional velocity-depth pairs by linearly extrapolating up to the surface and down to cutoff depth $z_{c}$. The extrapolation is performed based on a linear fit to $U_{\mathrm{EDM}}^{\mathrm{init}}(z)$ over a depth interval $\delta z_{T}$ and $\delta z_{B}$ at the shallow and deep end of the regime of mapped depths respectively, denoted in Figure 1b. The extrapolated points are shown as the black squares. 3. 3.

Perform a second polynomial fit on the expanded set of points (also of order $n_{\mathrm{max}}$) to produce the profile considered to be $U_{\mathrm{EDM}}$. 4. 4.

Scale polynomial coefficients defining $U_{\mathrm{EDM}}$ by $n!$ as in (8) to produce a profile $U_{\mathrm{PEDM}}^{\mathrm{init}}(z)$. 5. 5.

Create a new set of linearly extrapolated points down to $z_{c}$ based on a linear fit to $U_{\mathrm{PEDM}}^{\mathrm{init}}(z)$ over a depth region $\delta z_{B}/2$ at the deep end of the regime of mapped depths. Extrapolation is not performed up to the surface (thus differing from step 2). 6. 6.

Perform a final polynomial fit on the set of points including $U_{\mathrm{PEDM}}^{\mathrm{init}}(Z_{\mathrm{eff}}(k_{i}))$ and the extrapolated points in Step 5, to produce $U_{\mathrm{PEDM}}$.

The final current profile may be dependent on the parameters $n_{\mathrm{max}}$, $\delta z_{T}$, $\delta z_{B}$, as well as $z_{c}$, and a method for choosing optimal values of these parameters is necessary. To proceed, we note that when the exact form of the current profile $U(z)$ is considered, the Doppler shifts calculated using (2) or another suitable approximation method will agree with the measured values save for experimental measurement errors. The process of calculating the Doppler shifts given a prescribed current profile we refer to as the “forward problem.” Though the accuracy of (2) and its finite depth counterpart (Skop, 1987) is likely sufficient, we use a direct integration method of arbitrary accuracy due to Li & Ellingsen (2019) to evaluate the Doppler shifts to avoid this unnecessary source of error. We define an RMS difference between the measured Doppler shifts and those calculated by the forward problem ($\tilde{c}_{F,i}$) as

[TABLE]

where the overbar represents an average over all wavenumbers. For accurate evaluation of $\tilde{c}_{F,i}$, the cutoff depth was chosen as $z_{c}=2(\mathrm{min}[k_{i}])^{-1}$ (four times the deepest mapped depth), being set to the water depth in cases where the bottom was shallower than $z_{c}$.

Values of $n_{\mathrm{max}}$, $\delta z_{T}$, and $\delta z_{B}$ were in practice chosen to minimize $\epsilon_{\mathrm{RMS}}$ to in a sense find the most probable current profile in the presence of experimental noise.

3 Experimental and Data Analysis Methods

We test and evaluate the accuracy of the PEDM on experimental data by measuring wave spectra of waves propagating atop a controlled background shear flow generated in a small-scale laboratory setup, shown in Figure 2. The current depth-profile of the shear flow is measured by particle image velocimetry (PIV), which can be used as “truth” data to compare against the profiles obtained by the PEDM.

The setup consists of a pump which drives laminar flow over a 2x2 meter transparent plate, where different shear profiles can be obtained by various methods of flow conditioning. One method consists of a sequence of honeycomb structures and a curved wire mesh (Dunn & Tavoularis, 2007), which distorts the streamlines of the flow producing a profile with peak velocity at the surface, and decreasing with depth with approximately constant shear. The surface current and near-surface shear strength can be controlled by adjusting the water depth and pump power. It is noted that strong shear near the bottom due to the boundary layer is also created, yet for the depths ($\sim 8-10$ cm) and wavelengths we consider the influence of the boundary layer on wave dispersion is negligible. Another method is to make use of a region of flow where the water surface is nearly stagnant (at rest in the laboratory frame of reference) which occurs near the downstream end of the system due to the formation of a Reynolds ridge from surface contaminants (Scott, 1982). The region exhibits strong near-surface shear as the incoming flow dips beneath the stationary viscoelastic surface layer to form a surface boundary layer. The upstream extent of this stagnation region can be increased by the insertion of a horizontal bar in the downstream end as shown in Figure 2. A laboratory coordinate system is defined as shown in Figure 2, with the $x$, $y$ and $z$-axes aligned with the streamwise, spanwise and vertical dimensions respectively.

The depth profile of the shear flow was measured at varying locations in the streamwise and spanwise directions using a planar PIV setup with high power light-emitting diodes (LEDs) as the illumination source similar to the system of Willert et al. (2010). Emission from the LED’s (Luminus PB-120) was approximately collimated in one dimension to produce a planar light sheet using either a fiber bundle splayed out into a linear array and a cylindrical lens, or a thin rectangular slit mounted above the LED array. The water was seeded with 40 $\mu$m diameter polystyrene spheres (Microbeads AS), and particle images were acquired by a camera mounted out of the plane as shown. Image pairs were processed to obtain the streamwise velocity as a function of depth. The setup could be translated to perform measurements at different positions in both horizontal dimensions.

Waves were created using a vertical piston wavemaker mounted at the upstream end of the setup. The wavemaker was run at variable frequencies from 1 to 4 Hz as a function of time, 10 s at each constant frequency in steps of 0.1 Hz, to produce a sufficiently wide spectrum in frequency-wavevector space. The waves were measured using a synthetic Schlieren (SS) method (Moisy et al., 2009), consisting of a random dot pattern mounted below the transparent bottom plate, and viewed from above by a camera mounted $\sim$ 2 m optical path length from the free surface. The gradient of the free surface, $\nabla\eta(x,y,t)\equiv[\eta_{x}(x,y,t),\eta_{y}(x,y,t)]$, can be found by digital image correlation (DIC), comparing camera frames of the dot pattern beneath a perturbed free surface to that of an unperturbed reference frame. Uncertainty in the measured gradients was estimated to be 0.001 based on analysis of images taken with a still water surface. Typical measured root mean squared (RMS) gradients of the waves were between 0.02 – 0.1 in magnitude, resulting in a relative uncertainty of 5% or less.

The frequency-wavevector spectrum of the wave gradient field in a 10 s time window roughly corresponding to a given driven wavemaker frequency was calculated as

[TABLE]

where $P_{x}^{l}$ and $P_{y}^{l}$ are the three dimensional discrete Fourier transforms in spatial and temporal dimensions of the surface gradient components obtained directly from the SS method, which are first multiplied with a spatiotemporal windowing filter prior to transformation,

[TABLE]

where $L_{x}$ and $L_{y}\sim 0.5$ m are the physical lengths of the spatial domain, $T=10$ s the extent of the temporal domain, and $\sigma_{m}=1/4$. The spatiotemporal domain is assumed to be centered around $\{x,y,t\}=0$ such that $F(x,y,t)$ is peaked in the center of the domain. The spectra $P^{l}$ for each time window were summed together to produce a single spectrum $P=\sum_{l}P^{l}$ containing all frequency spectral components. For the purposes in this work, the fact that the wave spectrum is defined with the free surface gradient instead of the free surface elevation is insignificant, since the gradient field has the same periodicity in space and time as the surface elevation.

Assuming small wave-steepness, maximum values of the gradient spectrum $P$ are concentrated on the linear dispersion surface $\omega_{\mathrm{DR}}(k_{x},k_{y})$, which was assumed to be the sum of two components, a quiescent water term and a term due to the subsurface flow (1). The quiescent water dispersion relation $\omega_{0}(k)$ is of the form

[TABLE]

with $g$ the gravitational constant, $\sigma$ the surface tension constant, and $\rho$ the water density. The surface tension coefficient depends on the level of contamination of the water, and was determined by analyzing the wave spectrum recorded with the pump turned off using a pneumatic wavemaker discharging bursts of air at controlled frequencies of 5-10 Hz. A set of frequency-wavenumber pairs $\{k_{i},\omega_{i}\}$ were extracted by finding the peak wavenumber $k_{i}$ of the spectrum along various directions in $\mathbf{k}$ space for a given frequency $\omega_{i}$. The set of points was then fit to (12) with $\Gamma\equiv\sigma/\rho$ the fitting parameter. For the stagnation region flows, contaminants become concentrated in the viscoelastic surface layer, and thus a notably different value of the surface tension coefficient may result when compared to quiescent waters where the contaminants disperse over the whole water channel surface. To obtain a representative value of the surface tension in the stagnation region, we insert horizontal bars dipping just below the surface at the upstream and downstream boundaries of the measurement region and spanning the entire width of the channel in the $y$-direction prior to turning the pump off. The bars prevent the spreading of the surface contaminants over the entire channel region when the pump is turned off. Within the stagnation region there is in fact a gradient in surface tension in the streamwise direction, necessary to balance the surface shear stress of the fluid (Harper & Dixon, 1974). Using values of the viscosity in clean water and the maximum measured surface shear based on profiles measured by PIV, we estimate the variation of the surface tension coefficient $\sigma$ to be 0.008 $\mathrm{Nm}^{-1}$ across the measurement area, or $8\times 10^{-6}$ $\mathrm{m}^{3}\mathrm{s}^{-2}$ in the value of $\Gamma$, a relative variation of $\sim 20$%. We assume the measurements of $\Gamma$ using the method described above to be representative of the spatially averaged value within the measurement region. The effect of the inaccuracy thus introduced on our results will be discussed shortly.

The process here of determining the surface tension coefficient $\Gamma$ is specific to the small-scale laboratory setup, as in most practical cases in the field the length scales of the measured waves are in the regime where surface tension forces can be neglected. In cases when investigating short wavelengths in the ocean (e.g., Laxague et al., 2017), a reasonable estimate to the surface tension coefficient and density can be assumed a priori.

Both (1) and (12) describe wave propagation assuming fluid properties ($\Gamma$, $h$, and $\mathbf{U}$) to be invariant across horizontal spatial dimensions. However, for the case of the stagnation region flows, both $\Gamma$ and $U(z)$ vary across the streamwise dimension, due to the surface shear stress balance and the development of the surface boundary shear layer respectively. To quantify the effect these variations have on wave dispersion, we calculate the difference in phase velocities for a wave propagating at the upstream versus downstream ends of the measurement region. For the case of surface tension, we assume $\Gamma$ to vary by $8\times 10^{-6}$ $\mathrm{m}^{3}\mathrm{s}^{-2}$, and for the surface boundary layer, the difference between the minimum and maximum values of the measured streamwise velocity measured in upstream versus downstream positions for the strongest shear current. The results are shown in Figure 3a as a function of wavenumber for waves propagating downstream (similar trends occur for upstream propagating waves). The variation of the current profile results in a greater variation in phase velocities ($\sim$ 20 mm/s) across the measurement region compared to surface tension gradients where the variation is $\leq$10 mm/s for the wavenumber range shown. The values in Figure 3a place a bounds on potential variations and uncertainties of the extracted wave Doppler shifts $\tilde{\mathbf{c}}(k)$, though it is expected that Doppler shifts will be representative of the spatially averaged values across the measurement region. For current profiles produced with the curved mesh configuration, significantly less variation across the measurement region is expected once the shear profile has reached a stable state within the measurement area, and there is in this case no streamwise gradient in surface tension.

The Doppler shift velocities as a function of wavenumber were extracted by analyzing the gradient spectrum spectrum $P$. The range of wavenumbers to consider was chosen based on the azimuthally-averaged wave number spectrum:

[TABLE]

where $\mathbf{k}=\{k\cos\theta,k\sin\theta\}$ is defined in polar coordinates here where $\theta$ is the angle in the $x,y$ plane from the positive $x$-axis. The spectra of waves atop the four current profiles considered here are shown in Figure 3b as a function of wavenumber, scaled by their maximum value. The wavenumber range for extraction of Doppler shifts was chosen to be wavenumbers where $S(k)$ was greater than 0.1 of the peak value for wavenumbers less than the peak value, and greater than 0.02 of the peak value for wavenumbers larger than the peak value. The minimum wavenumber was $\sim$20 $\mathrm{rad}\cdot\mathrm{m}^{-1}$ for all profiles, and the maximum between $\sim$120-190 $\mathrm{rad}\cdot\mathrm{m}^{-1}$. A set of wavenumbers $k_{i}$ was specified spanning minimum to maximum values in steps of $2\pi/(10L_{x})$.

For each wavenumber $k_{i}$, Doppler shift velocities were extracted by considering the signal $P(\mathbf{k},\omega)$ on a cylindrical surface of constant wavenumber magnitude $k_{i}$ in $(\mathbf{k},\omega)$ space, and using an NSP method (Huang et al., 2016; Huang & Gill, 2012; Serafino et al., 2010). The cylindrical surface as well as the dispersion surface from (1) is shown in Figure 4a for the case of a depth-uniform current. The method works to effectively determine the frequency of intersection $\omega_{\mathrm{DR}}(k_{i},\theta)$ as a function of $\theta$ between the cylindrical surface and the dispersion relation surface (which corresponds to peak values of $P$), where the wavevector arguments of $\omega_{\mathrm{DR}}$ are expressed in polar coordinates. From (1), it is apparent that in quiescent waters ($\tilde{\mathbf{c}}(k)=0$) the frequency of intersection is independent of azimuth angle, whereas in the presence of a current there is an additional oscillating component with amplitude and phase determined by $\tilde{\mathbf{c}}(k)$, as seen in Figure 4a as the dashed curve.

We proceed by finding Doppler velocity components $\tilde{c}_{x,i}$ and $\tilde{c}_{y,i}$ at wavenumber $k_{i}$. First, we define a characteristic function

[TABLE]

where $a=(\sigma_{m}T)^{-2}$ is defined based on the Gaussian width in Fourier space given the applied spatial Gaussian filter $F$ defined in equation (11). Dependence on $\tilde{c}_{x,i}$ and $\tilde{c}_{y,i}$ is implicitly included in $\omega_{\mathrm{DR}}$. In addition, we consider the second harmonic spectral components $\{2\mathbf{k},2\omega\}$ and define a modified spectrum

[TABLE]

where $P_{i}^{\prime}$ is then scaled such that the minimum value is zero. The signal at the higher harmonic is due to the weak non-linearity of the measured surface waves as well as non-linearities in the SS measurement system (Senet et al., 2001). Assuming the spectral peak associated with the second harmonic has comparable spectral width as the fundamental harmonic, the contribution to the peak from the second term in (15) would have a smaller width in $\theta$-$\omega$ space due to the factor two in the argument of $P$. Including the second harmonic may thus increase the sensitivity to currents by making the peak of $P_{i}^{\prime}$ more localized. Example values of $P_{i}^{\prime}$ on cylindrical surfaces of constant wavenumber are shown in Figure 4b and c for $k_{i}$ = 75 and 125 $\mathrm{rad}\cdot\mathrm{m}^{-1}$ respectively, as a function of $\theta$ and $\omega$ for waves atop a shear current. In both cases, the peak frequency as a function of $\theta$ displays a clear oscillatory trend due to the presence of shear as expected.

We find the Doppler shift velocities by maximizing the scalar product $N$ between $G$ and $P_{i}^{\prime}$:

[TABLE]

where $\langle...\rangle$ indicates a double integral over $\theta$ and $\omega$ (the same integral as (13)). To avoid local maxima other than those associated with the dispersion relation, $N$ is first evaluated on a grid of points spanning expected values of the Doppler shift velocity components, with the Doppler shifts corresponding the maximum value of $N$ used as an initial guess for further optimization. The resulting curves $\omega_{\mathrm{DR}}(k_{i}\cos\theta,k_{i}\sin\theta)$ from the fitting routine are shown as the dashed lines in Figure 4b-c. The dotted lines show the frequency in quiescent waters. As can be seen, there is a distinct departure in the peak signal as a function of angle that is captured by the NSP fit, but inconsistent with the quiescent water frequency as it should be. The Doppler shifts as a function of wavenumber are expected to display a smooth functionality based on (2), and values were removed using an outlier filter. Both components were fit to a first order polynomial to produce functions $\tilde{c}_{x}^{O}(k)$ and $\tilde{c}_{y}^{O}(k)$. Outliers were identified by considering the set $\{\tilde{c}_{x,i}-\tilde{c}_{x}^{O}(k_{i})\}$ (and the equivalent for the $y$-direction) and data lying more than 1.5 times the interquartile range below the first quartile and the same interval above the third quartile were removed.

4 Results and Discussion

To validate and examine the accuracy of the PEDM, we apply it to Doppler shifts measured in the laboratory setup with the current profile measured by PIV used as “truth” measurements to compare against. We consider experimental data for waves atop 4 different shear flows, referred to as profiles a-d), shown in Figure 5. Profiles a-c) are in the stagnation region at different flow rates which lead to varying near surface shear strengths and curvature. Profile d) was produced using the curved wire mesh, and had weaker surface shear strength and near-constant vorticity with depth. The parameters including the measured surface tension coefficient $\Gamma$ are given in Table 1. The velocity was not measured for the bottom $\sim$1-2 cm depth where the bottom boundary layer was located.

The measured Doppler shifts using the NSP method described in section 3 for the four profiles are shown in Figure 6. Doppler shifts $\tilde{c}_{x}^{F}(k)$ calculated with theory assuming the measured PIV profile are shown as the dashed curves. As no mean flow in the $y$-direction was expected, the true values of the $y$-components of the Doppler shifts $\tilde{c}_{y,i}$ were assumed to be zero at all depths. Differences between experiment and theory are $\leq 1$ cm/s over most wavenumbers, except for a 1-2 cm/s bias for profile a). The reason for the bias is not known, yet could be a result of a greater streamwise variation in the shear profile given that the pump power was greatest for this profile.

The PEDM was implemented as described in section 2.1.2 for each component of the Doppler shifts separately. Current profiles $U_{\mathrm{PEDM}}(z)$ were calculated with 900 combinations of the parameters $n_{\mathrm{max}}$, $\delta z_{T}$, and $\delta z_{B}$: $10$ values of $\delta z_{T}$ and $\delta z_{B}$ each, ranging from $0.5$-$4$ mm and $1$-$20$ mm, respectively, and $9$ values of $n_{\mathrm{max}}$ ranging from $2$-$10$. For each combination, the RMS difference $\epsilon_{\mathrm{RMS}}$ between the measured Doppler shifts and those from the forward problem with $U_{\mathrm{PEDM}}(z)$ was evaluated. Profiles where the initial polynomial fit $U_{\mathrm{EDM}}(z)$ was not monotonic were discarded. The combination of parameters that gave the lowest value of $\epsilon_{\mathrm{RMS}}$ were used to produce a profile that was presumed to be the most probable estimate.

The monotonic assumption was based on the fact that the Doppler shifts (of which $U_{\mathrm{EDM}}(z)$ is based) can be viewed as a weighted average of the current depth-profile, thus resulting in a large degree of smoothing of oscillations in the true profile when considering $U_{\mathrm{EDM}}(z)$ obtained from the mapped depths. Over a finite range of wavenumbers, it is assumed that the true Doppler shifts are monotonic for most all realistic current profiles, and that profiles $U_{\mathrm{EDM}}(z)$ that are not monotonic result from errors in the Doppler shifts. It is however important to note that the monotonic assumption here does not also constrain the profile $U_{\mathrm{PEDM}}(z)$, given the scaling of the polynomial coefficients.

The process of calculating $U_{\mathrm{PEDM}}(z)$ profiles and evaluating $\epsilon_{\mathrm{RMS}}$ for the 900 combinations of PEDM parameters with roughly 100 wavenumber-Doppler shift pairs took approximately 6 minutes on an Intel®CoreTM i7-4770 3.40 GHz processor with 32 GB of RAM. However, the vast majority of time was spent evaluating $\epsilon_{\mathrm{RMS}}$ using the direct integration method. It is noted that for cases where all wavenumbers can be assumed to be in deep water and the approximation accuracy of (2) is deemed sufficient, (5) may be used to evaluate $\tilde{c}_{F,i}$ directly from the PEDM polynomial coefficients. When using (5), the same process took only 16 s.

The results of applying the PEDM to the $x$-components of the Doppler shifts for the four profiles are shown in Figure 7. The black curve denotes the current profile as measured by PIV, the average over the spatial locations within the wave measurement area with the error bars denoting the maximum and minimum values measured by PIV over the spatial locations. Profiles $U_{\mathrm{EDM}}(z)$ and $U_{\mathrm{PEDM}}(z)$ using the optimal set of parameters are shown as the dash-dotted and dashed curves respectively, along with the mapped Doppler shifts. For profiles a-c), the PEDM is a clear improvement over the EDM with notably increased accuracy over most all depths. Given the relatively strong curvature of the profiles, the assumption of a linear profile that was inherent in the mapping function is not valid here, and the mapped Doppler shifts deviate notably compared to the measured current profile. The deviation is greatest for profile a) and successively decreases for profiles b) and c) which is expected based on the weakened curvature of these profiles. For profile d) where the true profile has near-constant vorticity, the assumption of a linear profile is largely valid and the PEDM offers negligible improvement in accuracy over the EDM as may be expected. The shaded regions are discussed shortly.

To evaluate the improvement in accuracy of the PEDM, we calculate the depth-integrated RMS difference $\Delta U_{\mathrm{RMS}}$ between $U_{\mathrm{PEDM}}(z)$ or $U_{\mathrm{EDM}}(z)$ and the profile measured by PIV over the range of mapped depths. The results are summarized in Table 2, along with the optimal PEDM parameters for each profile. The ratio shown in the rightmost column is the degree of improvement in accuracy achieved by the PEDM relative to the EDM. An improvement of $>3\times$ is achieved for profiles a-c), with a maximum improvement of 5.1$\times$ for profile b). For profile d), the PEDM is marginally less accurate than the EDM, yet the absolute value of $\Delta U_{\mathrm{RMS}}$ remains small compared to the other profiles. For all profiles, the PEDM achieves a depth-integrated RMS absolute accuracy $<10$ mm/s relative to the PIV profiles.

4.1 Dependence on PEDM parameters

By using the combination of parameters $n_{\mathrm{max}}$, $\delta z_{T}$, and $\delta z_{B}$ that give the minimum $\epsilon_{\mathrm{RMS}}$ value, the values are thus set algorithmically during the running of the PEDM “algorithm” rather than as a required input determined prior to it. Thus from a user perspective the method is made effectively parameter free as we will now explain. It is noted that the same parameters are necessary in the use of the EDM as well, in creating a smooth velocity profile to fit the set of mapped Doppler shifts.

We examine the dependence of the results on the choice of the PEDM parameters by calculating $\Delta U_{\mathrm{RMS}}$ for each combination of parameters for both $U_{\mathrm{EDM}}(z)$ and $U_{\mathrm{PEDM}}(z)$, and plotting $\Delta U_{\mathrm{RMS}}$ against $\epsilon_{\mathrm{RMS}}$ as is shown in Figure 8 for the four current profiles. Also shown are results assuming a depth-uniform profile ($n_{\mathrm{max}}=0$) and constant shear ($n_{\mathrm{max}}=1$) which are independent of the choice of $\delta z_{T}$ and $\delta z_{B}$. It is noteworthy that $\Delta U_{\mathrm{RMS}}$ cannot be evaluated in realistic situations where “truth” measurements do not exist, so a criteria for choosing the optimal set of PEDM parameters to achieve a small value of $\Delta U_{\mathrm{RMS}}$ is desired based on metrics such as $\epsilon_{\mathrm{RMS}}$ that may be readily evaluated purely from the wave spectral data. The parameter combinations resulting in the minimum value of $\epsilon_{\mathrm{RMS}}$ are outlined with the open green squares in Figure 8, corresponding to the profiles $U_{\mathrm{PEDM}}$ shown in Figure 7.

Ideally, there would be a strong correlation between small values of $\epsilon_{\mathrm{RMS}}$, which can be calculated from the experimental data only, and $\Delta U_{\mathrm{RMS}}$ for which it is our goal to minimize. In Figure 8a-b) for the profiles with the strongest curvature, there is noticeable correlation for the smallest values of $\epsilon_{\mathrm{RMS}}$. For those cases, various values of $\epsilon_{\mathrm{RMS}}$ all yield values of $\Delta U_{\mathrm{RMS}}$ that are a significant improvement over the EDM cases (shown as the circles). It is notable that for profile b) where the PEDM profile with lowest value of $\epsilon_{\mathrm{RMS}}$ yielded a 5.1$\times$ reduction in $\Delta U_{\mathrm{RMS}}$ relative to the EDM, other points near the minimum $\epsilon_{\mathrm{RMS}}$ value still give a $\sim$3X or greater improvement in accuracy (the same being true for profile a)). For profiles c) and d), there is significantly less correlation between $\epsilon_{\mathrm{RMS}}$ and $\Delta U_{\mathrm{RMS}}$. Nonetheless, for profile c), the minimum value of $\epsilon_{\mathrm{RMS}}$ yields a value of $\Delta U_{\mathrm{RMS}}$ that is notably less than that of the EDM and constant shear case. For profile d) there is no significant difference in $\Delta U_{\mathrm{RMS}}$ between the EDM, PEDM, and constant shear cases considering the smallest values of $\epsilon_{\mathrm{RMS}}$, which may be expected given the approximately linear form of the current profile. For all cases, there is a distinct improvement in accuracy relative to the depth-uniform assumption. In addition, for all profiles the EDM displayed a similar level of accuracy relative to the case of constant shear, which is reasonable given that the same assumption was inherent to the EDM. Furthermore, profiles a) and b) with the greatest degree of curvature display the largest improvement over the constant shear case considering the lowest values of $\epsilon_{\mathrm{RMS}}$.

Choosing the optimal set of PEDM parameters based on $\epsilon_{\mathrm{RMS}}$ in a sense can be considered to yield the most probable current profile, i.e. the profile that agrees to the greatest degree with the experimentally measured Doppler shifts. However, given experimental noise it is useful to examine the variation in current profiles for parameter combinations that yield values of $\epsilon_{\mathrm{RMS}}$ near the minimum value, as those profiles may be considered nearly as probable. We calculate the bounds on the range of current values as a function of depth considering all profiles where $\epsilon_{\mathrm{RMS}}$ is within 10% of the minimum value, and show these bounds as the shaded regions in Figure 7. For the stagnation region profiles a-c), the spread is narrowest for a-b) which may be expected based on the stronger correlation between $\epsilon_{\mathrm{RMS}}$ and $\Delta U_{\mathrm{RMS}}$ as shown in Figure 8, where small values of $\epsilon_{\mathrm{RMS}}$ yield a smaller spread in the values of $\Delta U_{\mathrm{RMS}}$. For profile c), the spread in $\Delta U_{\mathrm{RMS}}$ is much greater, and less accurate profiles with near constant shear are included. As Figure 8c shows, the lowest value of $\epsilon_{\mathrm{RMS}}$ is much closer to that of the EDM even though the improvement in $\Delta U_{\mathrm{RMS}}$ is very significant. Had the threshold for the shaded region been set lower, the least good, near-linear profiles would be excluded.

Another potential reason for the increased spread in profile c) is the fact that the measured Doppler shifts appear slightly less smooth as a function of wavenumber when compared to profiles a-b). Furthermore, for profile a) where the measured Doppler shifts displayed a bias relative to those calculated from theory yet are relatively smooth as a function of wavenumber, the PEDM results in a very narrow spread around the most probable current profile that also has a corresponding bias towards reduced current strength near the surface compared to the PIV profile.

The same procedure and data analysis is applied to the $y$-components of the measured Doppler shifts and shown in Figures 9 and 10. As there was expected to be no current in this direction for all cases, the results represent the case of a depth-uniform profile in a moving reference frame. As expected, there is negligible improvement in accuracy using the PEDM relative to the EDM. The results serve as further important confirmation that the PEDM results do not deviate significantly from the results of the EDM in cases where the assumptions of a linear profile are valid. As shown in Figure 10, assumption of constant shear results in roughly the minimum value of $\epsilon_{\mathrm{RMS}}$, with only a slight increase in $\Delta U_{\mathrm{RMS}}$ relative to the depth-uniform current assumption. Due to experimental noise, results for both the EDM and PEDM result in slightly sheared current profiles, yet absolute values of $\Delta U_{\mathrm{RMS}}$ remain $<$1 cm/s for all parameter combinations in the vicinity of the minimum $\epsilon_{\mathrm{RMS}}$ value. Note that the range of values of the horizontal current strength axis in Figure 9 is reduced compared to Figure 7.

4.2 Scalability and Applicability of the Results

Given the small scale of the laboratory setup and the use of a different method to measure the wave spectrum than what may be used in the field, some discussion of the scalability and applicability of the results reported herein is warranted.

The absolute accuracy achieved herein with the PEDM is related to the scale of the setup, as well as the characteristics of the wave spectrum. The more pertinent metric is the fractional improvement in accuracy relative to the EDM, which is expected to be scalable to larger measurement setups and different techniques of measuring the wave spectrum. The relative improvement using the PEDM is related to the form of the current profile. In cases where the profile is approximately linear over the range of depths, limited improvement is expected since the approximation to which the EDM’s mapping function was based is valid. In cases where the current profile has greater curvature near the surface, the PEDM is found to yield a greater fractional improvement in accuracy. The PEDM thus acts in a sense to improve the estimate to the current profile where possible, while performing similarly with the EDM otherwise. Note that the shape of the lab current profiles in the stagnation region, profiles a-c) in Figure 5, are representative of a scaled-down surface shear layer such as may be produced in the wind-swept ocean Ekman layer or in a river delta plume such as reported by Kilcher & Nash (2010). They differ in shape only by a constant subtraction of the deep-water velocity which corresponds to a constant offset in Doppler shifts.

4.2.1 Resolution

The absolute accuracy of the Doppler shifts is fundamentally determined by, among other factors, the size of the measurement domain $L_{x}$ which sets the resolution in $k$-space, \Delta k$$=2\pi/L_{x}. Herein, the extraction of the Doppler shift velocities was performed by evaluating the wave spectrum on a surface in spatiotemporal Fourier space with wavenumber $k$ kept constant, requiring interpolation between the available discrete values of $\{k_{x},k_{y}\}$. A smaller value of $\Delta k$ reduces errors due to interpolation and also decreases the spectral leakage from neighboring wavenumber components. In an attempt to bound the uncertainties in Doppler shifts caused by interpolation we define a velocity shift $\Delta c_{\Delta k}$ so that

[TABLE]

$\Delta c_{\Delta k}(k)$ is thus the depth-uniform current velocity that causes the linear dispersion surface to move by approximately one pixel in $k$-space for the relevant constant frequency $\omega_{0}(k)$. The values of $\Delta c_{\Delta k}$ over the range of wavenumbers where Doppler shifts were extracted are shown in Figure 11 (see also Figure 6 for the Doppler shift wavenumber range). Another source of uncertainty in the Doppler shift involves the spread of the spectrum in frequency space, related to $\Delta\omega=2\pi/T$, where $T$ is the total measurement period. Again, we transform this quantity to a velocity:

[TABLE]

which is also shown in Figure 11. Given the sizes of our measurement domain in space and time, $\Delta c_{\Delta k}$ is nearly an order of magnitude greater than $\Delta c_{\Delta\omega}$ over the range of relevant wavenumbers, indicating that resolution in $k$-space is the main contribution to uncertainties in the Doppler shifts. Examining the figure gives an estimate to the upper bounds to the uncertainties that can be expected in the Doppler shifts, and similarly the reconstructed profiles due to the finite spectral resolution. Comparing the values of $\Delta c_{\Delta k}$ to the values of $\Delta U_{\mathrm{RMS}}$ from the PEDM, it is evident that a great degree of sub-pixel resolution is achieved using the NSP and PEDM methods: $\Delta U_{\mathrm{RMS}}$ is less than the minimum value of $\Delta c_{\Delta k}$ for all current profiles, being orders of magnitude less than values of $\Delta c_{\Delta k}$ for the lower wavenumbers.

We note that values of $\Delta c_{\Delta k}$ and $\Delta c_{\Delta\omega}$ for full-scale measurements in the ocean using for example X-band radar are typically within an order-of-magnitude of the values shown in Figure 11, assuming spatial domain size $L_{x}\sim 750$ m, $T\sim 10$ min, and wavenumbers in the range of $0.05-0.3$ $\mathrm{rad}\cdot\mathrm{m}^{-1}$ as is common (e.g. Lund et al., 2015). Thus, though values of $\Delta U_{\mathrm{RMS}}$ from measurements in the ocean at large scales are expected to be larger than those reported here, it is not expected that the errors will increase by orders-of-magnitude.

4.2.2 Scalability

Consider now how the small-scale experimental setup scales up to an oceanographic scale. First, it is obvious that the one effect which does not scale up, is that of surface tension, which is utterly negligible at the wavelengths measurable with e.g. X-band radar. In our experiment we do observe Bond number $\rho g\lambda^{2}/\sigma\lesssim\mathcal{O}(1)$ at the shortest wavelengths, yet the majority of our spectrum lies in the gravity wave regime, thus being physically directly comparable. This said, the PEDM method is not sensitive to whether or not the dispersion relation has capillary corrections at high $k$, and so the stringency of our testing is little altered by this.

Assuming wavelengths to lie in the gravity wave regime, and assuming essentially infinite depth as is approximately true of our experiment, the system scales in the following way. Now only a single nondimensional group remains, a shear-Froude number based on three physical parameters: a typical wavelength of the spectrum, $g$, and a suitably defined depth-averaged shear. A suitable definition is

[TABLE]

referred to as $\delta$ by Ellingsen & Li (2017). $\langle S\rangle_{k}$ is the depth averaged shear along $\mathbf{k}$ suitably weighted for wave number $k$. Full similarity can be obtained if, by scaling up the velocity profile to oceanographic scale, the range of important $k$-values in the wave spectrum yields the same values of $\mathrm{Fr}_{S}$. Let’s assume $U(z)$ is the lab current, and an oceanographic current of the same shape is $U_{O}(z)=u^{*}U(\delta z)$ with $\delta$ a small parameter describing the slower variation with depth and $u^{*}$ the fraction of the velocities at $z=0$. To probe the velocity profile into the depth in a similar manner as before, a lower wave number (i.e. longer wavelength) $k^{\prime}=\delta k$ is required. On the whole we obtain $\mathrm{Fr}_{S}\to u^{*}\sqrt{\delta}\mathrm{Fr}_{S}$. In other words, similarity is in order if $u^{*}\sqrt{\delta}\sim\mathcal{O}(1)$.

Our most strongly sheared velocity profile, in Figure 7a), resembles in shape and magnitude a very strong oceanographic velocity profile, such as that can be found in the Columbia River delta (Kilcher & Nash, 2010), if we let $\delta=1/500$ and $u^{*}=12$, for example, resulting in $u^{*}\sqrt{\delta}\sim 0.54$ and shear-Froude numbers of the same order of magnitude. Wavelengths $500$ times those of the lab are reasonable for waves in the area, between $8$ and $80$ m for the wave numbers of Figure 6. Hence we conclude that, while the strongest shear tested in the lab is a little stronger than can be expected of a particularly strong scaled-up equivalent, it is a satisfactory test of the PEDM theory in realistic settings. Given the ease of high quality flow measurements, scaled-down lab experiments thus offer an ideal test-bed for studies of ocean wave propagation on shear currents.

We now comment on the range of depths at which the near-surface current profile is estimated. The depth range is determined directly by the range of mapped depths, and hence the range of wavenumbers in the measured spectrum. Though the choice of the mapping function is in a sense arbitrary, we argue the choice is reasonable based on intuition considering (2). At a depth $(2k)^{-1}$ the cumulative integral of the weighting function $2k\mathrm{e}^{2kz}$ is 0.63, i.e. a wave is influenced by roughly comparable amounts by currents at greater vs. shallower depths, indicating a reasonable choice of the depth assignment for most current profiles. Given the rapidly decreasing sensitivity of waves to currents at greater depths, the polynomial fits of the PEDM can be considered to be an expansion of the near-surface current profile in the top layer of the water column, valid over the depth range of the mapped Doppler shifts. As is well-known with polynomial fits, large errors can result with extrapolation for prediction of currents at greater depths. In the laboratory experiments reported here, the depth range of the reconstructed flow is only a few centimeters, while in the ocean with wave spectra measured by X-band radar the depths may extend to tens of meters, given a roughly three orders of magnitude increase in the scale of the measured wavenumbers.

4.2.3 Wave spectrum measurement

For the laboratory results presented here, the wave spectrum was measured using a synthetic Schlieren method which measures directly the gradient components of the free surface, differing from methods that are practical for field measurements on a larger scale. However, for the purposes of inversion methods, all that is required is a signal that has the same periodicity in space and times as the wave spectrum. As mentioned in the introduction, various methods of measuring the wave spectrum in the radar and optical regime have already been used in reconstructing near surface currents. The choice of the wave spectral measurement method affects primarily the range of wavenumbers that are probed and is relatively inconsequential in terms of the inversion method process, affecting only the details of extraction of the Doppler shifts. A main difference between field measurement techniques such as X-band radar and the SS is that the mapping of free surface elevation to measured signal is, to a greater degree, nonlinear. The nonlinearities result in a signal at higher harmonics in the wave spectrum, yet the fundamental harmonic has the same periodicity in space and time as true wave component. Furthermore, the NSP method uses the signal at the second harmonic in determining the Doppler shifts. Thus, though the SS method employed in this work is impractical to be used in field measurements at larger scales, it can be viewed as an equivalent technique to those used in the field for the purposes here of fundamentally studying inversion methods.

4.2.4 Applications

The PEDM method may be applied to Doppler shifts extracted from wave spectra obtained by observation techniques readily available with today’s technology, such as X-band radar or optical images of the ocean surface, as well as potential future methods for remotely sensing the directional wave spectrum. The Doppler shifts may be extracted by a number of means such as least squares techniques or the NSP method described and further developed herein.

As demonstrated in figure 7, the PEDM offers greatest improvement in accuracy over the EDM in cases where the current profile has strong near-surface curvature within the range of mapped Doppler shifts. For the case of wave spectra measured by X-band radar where the mapped depths may typically be on the order of 2-10 m (e.g. Lund et al., 2015), current profiles with strong curvature are expected to occur in times of high winds, and at specific locations such as river deltas with strong shear currents driven by density differences in the fluid (e.g. Kilcher & Nash, 2010). Use of the PEDM to achieve a more accurate current depth-profile under such circumstances could result in improved characterization of submesoscale currents (Lund et al., 2018),improved estimates of wave steepness for predicting breaking waves (Zippel & Thomson, 2017), and improved mapping of shear currents for coastal engineering applications, for example. Under extreme sea states such as during hurricanes, improved accuracy in the reconstruction of remotely-sensed shear current profiles could allow for better prediction of wave and current forces on structures, where Dalrymple (1973) has shown that currents even with velocities small compared to the wave orbital velocities can result in a notable increase in the forces on structures. In the latter case, however, strong wave nonlinearity and imaging difficulties may make remote sensing difficult in practice.

For wave spectra measured using optical-based methods, the range of mapped depths is typically significantly shallower than for X-band radar data, in some cases resolving the top few centimeters of the water column where the current may have strong curvature even under moderate conditions (Laxague et al., 2018). The PEDM has the potential to improve the accuracy of the reconstruction in such cases, furthering applications such as studies of the air-sea interaction as well as the transport of contaminants near the ocean surface (Laxague et al., 2018).

In conditions where the current profile is approximately depth-uniform over the range of mapped depths, the PEDM is not expected to increase the accuracy of the reconstructed currents compared to the EDM or other existing methods which assume depth-uniform flow, yet figure 9 demonstrates that the PEDM gives essentially identical results in such cases, eschewing the need to employ different methods in different conditions. By being simple to employ and performing equally well or better than current methods, we propose that the PEDM can replace competing inversion methods in current use in most situations. The exception we can imagine is situations where calculation cost is a very severe restriction.

4.2.5 Limitations and challenges

As with all inversion methods, the absolute accuracy of the PEDM is affected by the wave spectrum bandwidth in terms of wavenumber and angular spread. Reconstruction of the depth profile of the flow places more stringent demands on the wave spectrum having a broader range of wavenumbers and directions, when compared to methods aimed at estimating a single (depth-uniform) velocity vector, given the additional fitting parameters associated with the PEDM method: the PEDM involves $n_{\mathrm{max}}+1$ polynomial coefficients for each horizontal dimension, whereas depth-uniform estimation requires only one. The need for a sufficient spectrum of waves to be present, however, is due to fundamental physics and will affect any method whereby currents are estimated from surface wave dispersion. If the currents have no surface imprint, clearly they simply cannot be inferred from surface measurement. Likewise, sufficient image quality is a fundamental requirement for all methods.

In addition, under some circumstances such as extreme sea states nonlinear wave interactions become more prevalent, in which case analysis of the wave spectrum becomes more complicated due to the presence of bound waves. The same complication has also been observed for moderate wave slopes in a wind wave tank (Laxague et al., 2017). Analyzing the wave spectrum to extract the Doppler shifts corresponding to currents when nonlinear wave interactions are prevalent requires further study.

The PEDM method, like other similar methods which it aspires to replace, assumes horizontally homogeneous currents. When the horizontal variation is not slow compared to all relevant wavelengths, such as will often be the case particularly in coastal areas, more advanced methods will be required, beyond the current state-of-the-art.

5 Conclusions

A new method for reconstructing near surface current profiles from measurements of the wave spectrum has been presented, demonstrated and carefully tested and compared to the state-of-the-art inversion method.

The method is easy to implement. It takes the present state-of-the art technique of assigning effective depths to measured Doppler shift velocities (the effective depth method, EDM) as its starting point. A polynomial fit is made to the EDM profile from whose coefficients a new velocity profile estimate of polynomial form is created via a simple derived relation. The resulting polynomial profile is an improved estimate to the true current profile compared to state-of-the-art methods such as the EDM as it does not make any a priori assumptions on the general shape of the profile, and involves very little added complexity.

Our new polynomial effective depth method (PEDM) was tested on data obtained from a laboratory setup where background currents of different depth profiles could be created in a controlled manner and measured independently using particle image velocimetry which was used as “truth” measurements. The laboratory setup is an ideal test-bed for further studies regarding remote sensing of near-surface shear currents given the large degree to which the current profile and wave spectrum can be controlled and the straightforward scalability of the results up to oceanic scales. The PEDM offers a $>3\times$ improvement in accuracy relative to the EDM for profiles with strong near-surface curvature. For cases where the true current profile has approximately constant shear, the assumptions upon which the EDM is based are fulfilled, and the PEDM offers limited improvement in accuracy. The estimate produced is then similar to that of the EDM in accuracy and shape, demonstrating the robustness of the method.

A simple criterion was developed to determine optimal values for parameters involved in the polynomial fits to achieve the most probable current profile estimate. The criterion depends on the measured Doppler shift data only, and thus the PEDM involves no free parameters. A novel adaptation of the normalized scalar product method (NSP) was developed to extract Doppler shifts from wave spectra at multiple wavenumbers, including the second harmonic of the spectrum.

The results indicate that the method can be applied to full scale field measurements to obtain higher accuracy in reconstructing near surface shear profiles from the wave spectrum, beneficial across a wide variety of oceanic applications.

Acknowledgements.

SÅE and AÅ  were funded by the Research Council of Norway (FRINATEK), grant no. 249740. Data and scripts used for the analysis are available at: https://doi.org/10.18710/8JBWCJ (doi: 10.18710/8JBWCJ).