An improved method for determining nearsurface currents from wave dispersion measurements
Abstract
A new inversion method for determining nearsurface shear currents from a measured wave spectrum is introduced. The method is straightforward to implement and starts from the existing stateoftheart technique of assigning effective depths to measured wavenumberdependent 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 userinput 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 improvement in accuracy relative to the stateoftheart for current profiles with significant nearsurface curvature. The experiments are dynamically similar to typical oceanographic flows such as winddrift profiles and our laboratory thus makes a suitable and eminently useful scale model of the reallife setting. Our results show that the new method can achieve improved accuracy in reconstructing nearsurface 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.
JGR:Oceans
Department of Energy and Process Engineering, Norwegian University of Science and Technology, Trondheim, Norway
Benjamin K. Smeltzerbenjamin.smeltzer@ntnu.no
Simple and more accurate method for reconstructing nearsurface current profiles from wave spectra
Laboratory setup where shear currents and waves can be wellcontrolled and measured
Novel adaptation of a method for extracting wavelengthdependent Doppler shifts from wave spectra
1 Introduction
Characterizing nearsurface ocean currents is of importance to a vast range of applications. At a fundamental scientific level, nearsurface 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 wavebody 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 opticalbased detection, with the potential for mapping currents over a larger area (multiple ) 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 wavebased nearsurface 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 Xband 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 wavebased 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 illposed 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 depthdependence of the current profile based on physical intuition have typically been imposed.
Despite the hardships, wavebased current measurements have been used in the field for many decades. The techniques involve reconstructing the nearsurface 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 depthdependence of the flow. In reviewing the previously developed inversion methods, we first consider the dispersion relation for smallamplitude linear waves propagating atop a depthvarying flow, which can be approximated as:
(1) 
where is the wave frequency, the frequency in quiescent waters, the wavevector, , and a wavenumberdependent Doppler shift velocity due to the background current. Presume is the undisturbed water surface and the bottom is found at with . We shall mostly work in the deep water regime where 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
(2) 
where 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% that at the surface at a depth equal to half the wavelength (). 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 obtained from experimental data to determine the unknown .
A word of warning is warranted when referring to as the “Doppler shift” as is conventional. While 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 , for which taking the dependence of into account is key. We shall follow the numenclatorial convention in the literature and refer to as the Doppler shift velocity while bearing this in mind.
Various wave detection methods are sensitive to different spectral ranges of 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 singlefrequency HF radar is often referred to as the surface current, yet is more precisely a weighted average of the current profile from (2), as it essentially measures ( being the wavenumber of the resonant wave) without information concerning the depthdependence. Depthprofile 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 Xband radar or optical techniques inherently measure a wide spectrum of wavelengths, thus evaluating (2) at many values and enabling the use of inversion methods to estimate the current depthdependence.
Inversion methods of determining from a set of measured values of at discrete wavenumbers can be carried out separately for each velocity component, i.e. can be found from , and from . To ease the notation, in the following we outline the new inversion method using and 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 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, i.e. . For the commonly assumed case of a current profile which varies linearly with depth, . By the approximation (2) the Doppler shift is approximated as^{1}^{1}1One should notice that the SJ approximation (2) is only accurate when ; Ellingsen & Li (c.f., e.g., 2017).
(3) 
The last form shows that assuming linear current, deep water and using the SJ approximation, the appropriate effective depth is
(4) 
(In other words is approximately of the wavelength.) A similar relation can also be derived for a logarithmic profile (Plant & Wright, 1980). Estimating from a measured using (3) or its sibling assuming a logarithmic profile we refer to as the effective depth method (EDM). The EDM has been used extensively in the literature for estimating nearsurface shear currents (e.g., Teague et al., 2001; Lund et al., 2015; Laxague et al., 2017; Laxague et al. , 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 values, which was further developed and applied to data from Xband 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 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, then makes use of a simple relation which follows from (2) to construct an improved estimate of the true profile 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 wellcontrolled 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 use of the new inversion method on experimentally measured Doppler shifts, where in situ measurements of the current profile are used as “truth” measurements to compare against. 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 at unique wavevector magnitudes 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 & Gill, 2012; Serafino et al., 2010) used herein (described in section 3).
Assuming a polynomial current profile of the form in deep water, evaluation of (2) yields the SJ approximation
(5) 
for the Doppler shift velocities. We notice that is equal to the mapping function used in the EDM assuming a linear frofile, equation (4). Using the EDM with this mapping the estimated current profile is
(6) 
Thus, the mapped profile is also of polynomial form with coefficients of the th order term differing by a factor from those of the true profile . The estimated velocity profile 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 cofficients in the series representation of and the true profile , namely that they differ by a factor .
Explicitly, the PEDM procedure consists of the following three steps:

From the measured values of , obtain using (3).

Fit the obtained values of to a polynomial of degree
(7) 
Then the improved PEDM estimate is
(8)
2.1 Theoretical limitations
Two notable potential complications arise: finite depth where (5) and (6) are no longer stictly 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 .
2.1.1 Performance for finite depth
In the case of finite depth , an explicit relation of the form of (6) cannot be derived since the mapping function in finite depth cannot be solved with respect to analytically. To examine the effect of finite depth on the accuracy of the PEDM, we consider an exponential profile of the form , with equal to a value twenty times that of the minimum wavenumber , and the surface current.
We consider the implementation of the PEDM in infinite depth with , simply using (8) with determined using the mapping function in deep water. The fractional depthintegrated root mean square (RMS) error between and was calculated for cases over a range of water depth values , with the results shown in Figure 1. For all but the shallowest depths considered here, the deep water mapping function results in errors at the 1% level. It is noted that for values of , the deep water mapping case results in depths being mapped below the bottom for some ; these points were discarded. For most all realistic combinations of water depth and relevant wavenumbers, Figure 1 indicates that the deep water mapping function and (6) yield sufficient accuracy.
2.1.2 Effect of limitations of measured Doppler shifts
As mentioned, the fact that is measured for a finite range of wavenumbers will affect accuracy. This is true of any inversion method for reconstructing from dispersion measurements.
To handle the realistic case of experimentally measured Doppler shifts at a finite range of wavenumbers, we perform a three step process: first fitting the mapped Doppler shifts to a polynomial of order , then creating additional velocitydepth pairs by extrapolating up to the surface and down to a cuttoff depth , and finally by performing a second polynomial fit on the expanded set of velocitydepth points. The extrapolation is performed based on the average shear of the initial polynomial fit current profile over a depth interval and at the shallow and deep end of the regime of mapped depths respectively. A set of velocities are calculated at intervals and up to and down respectively. The resulting polynomial coefficients from the second polynomial fit on the expanded set of points are considered the profile , and are then scaled by , whereby the above process is repeated with extrapolation only being performed at the deep end of the regime using a depth interval to produce a final set of polynomial coefficients, defining .
The final current profile may be dependent on the parameters , , , as well as , 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 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 () as
(9) 
For accurate evaluation of , the cuttoff depth was chosen as (four times the deepest mapped depth), being set to the water depth in cases where the bottom was shallower than .
Values of , , and were in practice chosen to minimize 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 smallscale laboratory setup, shown in Figure 2. The current depthprofile 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 nearsurface 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 ( 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 nearsurface 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 , and 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 lightemitting diodes (LEDs) as the illumination source similar to the system of Willert et al. (2010). Emission from the LED’s (Luminus PB120) 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 m 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 frequencywavevector 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 2 m optical path length from the free surface. The gradient of the free surface, , 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 frequencywavevector spectrum of the wave gradient field in a 10 s time window roughly corresponding to a given driven wavemaker frequency was calculated as
(10) 
where and 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,
(11) 
where and m are the physical lengths of the spatial domain, s the extent of the temporal domain, and . The spatiotemporal domain is assumed to be centered around such that is peaked in the center of the domain. The spectra for each time window were summed together to produce a single spectrum 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 wavesteepness, maximum values of the gradient spectrum are concentrated on the linear dispersion surface , 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 is of the form
(12) 
with the gravitational constant, the surface tension constant, and 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 510 Hz. A set of frequencywavenumber pairs were extracted by finding the peak wavenumber of the spectrum along various directions in space for a given frequency . The set of points was then fit to (12) with 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 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 viscocity in clean water and the maximum measured surface shear based on profiles measured by PIV, we estimate the variation of the surface tension coefficient to be 0.008 across the measurement area, or in the value of , a relative variation of %. We assume the measurements of 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 is specific to the smallscale 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 (, , and ) to be invariant across horizontal spatial dimensions. However, for the case of the stagnation region flows, both and 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 to vary by , 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 ( 20 mm/s) across the measurement region compared to surface tension gradients where the variation is 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 , 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 . The range of wavenumbers to consider was chosen based on the azimuthallyaveraged wave number spectrum:
(13) 
where is defined in polar coordinates here where is the angle in the plane from the positive 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 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 20 for all profiles, and the maximum between 120190 . A set of wavenumbers was specified spanning minimum to maximum values in steps of .
For each wavenumber , Doppler shift velocities were extracted by considering the signal on a cylindrical surface of constant wavenumber magnitude in space, and using an NSP method (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 depthuniform current. The method works to effectively determine the frequency of intersection as a function of between the cylindrical surface and the dispersion relation surface (which corresponds to peak values of ). From (1), it is apparent that in quiescent waters () 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 , as seen in Figure 4a as the dashed curve. Example values of on cylindrical surfaces of constant wavenumber are shown in Figure 4b and c for = 59 and 155 respectively, as a function of and for waves atop a shear current. In both cases, the peak frequency as a function of displays a clear oscillatory trend due to the presence of shear as expected.
We proceed by finding Doppler velocity components and at wavenumber . First, we define a characteristic function
(14) 
where is defined based on the Gaussian width in Fourier space given the applied spatial Gaussian filter defined in equation (11). Dependence on and is implicitly included in . In addition, to improve the sensitivity to currents we consider the second harmonic spectral components and define a modified spectrum
(15) 
where is scaled such that the minimum value is zero. The signal at the higher harmonic is due to the weak nonlinearity of the measured surface waves as well as nonlinearities in the SS measurement system (Senet et al., 2001).
We find the Doppler shift velocities by maximizing the scalar product between and :
(16) 
where indicates a double integral over and (the same integral as (13)). To avoid local maxima other than those associated with the dispersion relation, 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 used as an initial guess for further optimization. The resulting curves from the fitting routine are shown as the dashed lines in Figure 4bc. 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 and . Doppler shifts where (and the equivalent for the direction) were removed from the data set.
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 ad), shown in Figure 5. Profiles ac) 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 nearconstant vorticity with depth. The parameters including the measured surface tension coefficient are given in Table 1. The velocity was not measured for the bottom 12 cm depth where the bottom boundary layer was located.
Profile  Flow Type  Water Depth  Flow rate  

[mm]  s  
a  Stagnation region  95  0.021  
b  Stagnation region  95  0.017  
c  Stagnation region  95  0.014  
d  Curved mesh  80  0.014 
The measured Doppler shifts using the NSP method described in section 3 for the four profiles are shown in Figure 6. Doppler shifts calculated with theory assuming the measured PIV profile are shown as the dashed curves. As no mean flow in the direction was expected, the components of the Doppler shifts were assumed to be zero in theory. Differences between experiment and theory are cm/s over most wavenumbers, except for a 12 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 for each component of the Doppler shifts separately. Current profiles were calculated with 900 combinations of the parameters , , and : values of and each, ranging from  mm and  mm, respectively, and values of ranging from . For each combination, the RMS difference between the measured Doppler shifts and those from the forward problem with was evaluated. Profiles where the initial polynomial fit was not monotonic were discarded. The combination of parameters that gave the lowest value of 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 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 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 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 , given the scaling of the polynomial coefficients.
The results of applying the PEDM to the 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 errors bars denoting the maximum and minimum values measured by PIV over the spatial locations. Profiles and using the optimal set of parameters are shown as the dashdotted and dashed curves respectively, along with the mapped Doppler shifts. For profiles ac), 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 nearconstant 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 depthintegrated RMS difference between or 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 is achieved for profiles ac), with a maximum improvement of 13.7 for profile b). We consider this value to be potentially fortuitous, yet profiles a) and c) confirm the significant improvement in accuracy that may be expected for profiles with high degrees of nearsurface curvature. For all profiles, the PEDM achieves a depthintegrated RMS absolute accuracy mm/s relative to the PIV profiles.
Profile  

[mm]  [mm]  [mm/s]  [mm/s]  
a  8  0.5  17.9  32.7  9.3  3.5 
b  8  1.3  3.1  22.1  1.6  13.7 
c  10  1.7  5.2  14.2  3.8  3.7 
d  2  0.5  13.7  4.1  3.5  1.2 
4.1 Dependence on PEDM parameters
By using the combination of parameters , , and that give the minimum 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 for each combination of parameters for both and , and plotting against as is shown in Figure 8 for the four current profiles. Also shown are results assuming a depthuniform profile () and constant shear () which are independent of the choice of and . It is noteworthy that 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 is desired based on metrics such as that may be readily evaluated purely from the wave spectral data. The parameter combinations resulting in the minimum value of are outlined with the open green squares in Figure 8, corresponding to the profiles shown in Figure 7.
Ideally, there would be a strong correlation between small values of , which can be calculated from the experimental data only, and for which it is our goal to minimize. In Figure 8ab) for the profiles with the strongest curvature, there is noticeable correlation for the smallest values of . For those cases, various values of all yield values of 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 yielded a 13.7X reduction in relative the EDM, other points near the minimum value still give a 3X or greater improvement in accuracy (the same being true for profile a)). For profiles c) and d), there is significantly less correlation between and . Nonetheless, for profile c), the minimum value of yields a value of that is notably less than that of the EDM and constant shear case. For profile d) there is no significant difference in between the EDM, PEDM, and constant shear cases considering the smallest values of , 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 depthuniform 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 .
Choosing the optimal set of PEDM parameters based on 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 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 is within 10% of the minimum value, and show these bounds as the shaded regions in Figure 7. For the stagnation region profiles ac), the spread is narrowest for ab) which may be expected based on the stronger correlation between and as shown in Figure 8 where small values of yield a smaller spread in the values of . For profile c), the spread in is much greater, and less accurate profiles with near constant shear are included. As Figure 8c shows the lowest value of is much closer to that of the EDM even though the improvement in is very significant. Had the threshold for the shaded region been set lower, the least good, nearlinear profiles would be excluded.
Another potential reason for the increased spread in profile c) is the fact that the measured Doppler shifts appear less smooth as a function of wavenumber when compared to profiles ab). 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 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 depthuniform 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 , with only a slight increase in relative to the depthuniform current assumption. Due to experimental noise, results for both the EDM and PEDM result in slightly sheared current profiles, yet absolute values of remain 1 cm/s for all parameter combinations in the vicinity of the minimum 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 ac) in Figure 5, are representative of a scaleddown surface shear layer such as may be produced in the windswept 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 deepwater velocity which corresponds to a constant offset in Doppler shifts.
4.3 Error sources
The absolute accuracy of the Doppler shifts is fundamentally determined by, among other factors, the size of the measurement domain which sets the resolution in space, . Herein, the extraction of the Doppler shift velocities was performed by evaluating the wave spectrum on a surface in spatiotemporal Fourier space with wavenumber kept constant, requiring interpolation between the available discrete values of . Smaller values of reduces errors due to interpolation and also decreases the spectral leakage from neighboring wavenumber components. In an attempt to bound the errors in Doppler shifts caused by interpolation we define a velocity shift so that
(17) 
is thus the depthuniform current velocity that causes the linear dispersion surface to move by approximately one pixel in space for the relevant constant frequency . The values of 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 error in the Doppler shift involves the spread of the spectrum in frequency space, related to . Again, we transform this quantity to a velocity:
(18) 
which is also shown in Figure 11. Given the sizes of our measurement domain in space and time, is nearly an order of magnitude greater than over the range of relevant wavenumbers, indicating that resolution in space is the main contribution to errors in the Doppler shifts. Examining the figure gives an estimate to the upper bounds to the errors that can be expected in the Doppler shifts, and similarly the reconstructed profiles due to the finite spectral resolution. Comparing the values of to the values of from the PEDM, it is evident that a great degree of subpixel resolution is achieved using the NSP and PEDM methods: is less than the minimum value of for all current profiles, being orders of magnitude less than values of for the lower wavenumbers.
We note that values of and for fullscale measurements in the ocean using for example Xband radar are typically within an orderofmagnitude of the values shown in Figure 11, assuming m, min, and wavenumbers in the range of as is common (e.g. Lund et al., 2015). Thus, though values of 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 ordersofmagnitude.
4.4 Scalability
Consider now how the smallscale 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. Xband radar. In our experiment we do observe Bond number 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 , 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 shearFroude number based on three physical parameters: a typical wavelength of the spectrum, , and a suitably defined depthaveraged shear. A natural alternative is
(19) 
referred to as by Ellingsen & Li (2017). is the depth averaged shear along suitably weighted for wave number . Full similarity can be obtained if, by scaling up the velocity profile to oceanographic scale, the range of important values in the wave spectrum yields the same values of . Let’s assume is the lab current, and an oceanographic current of the same shape is with a small parameter describing the slower variation with depth and the fraction of the velocities at . To probe the velocity profile into the depth in a similar manner as before, a lower wave number (i.e. longer wavelength) is required. On the whole we obtain . In other words, similarity is in order if .
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 and , for example, resulting in and shearFroude numbers of the same order of magnitude. Wavelengths times those of the lab are reasonable for waves in the area, between and 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 scaledup equivalent, it is a satisfactory test of the PEDM theory in realistic settings. Given the ease of high quality flow measurements, scaleddown lab experiments thus offer an ideal testbed for studies of ocean wave propagation on shear currents.
We now comment on the range of depths at which the nearsurface 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 the cumulative integral of the weighting function 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 nearsurface current profile in the top layer of the water column, valid over the depth range of the mapped Doppler shifts. As is wellknown 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 Xband radar the depths may extend to tens of meters, given a roughly three orders of magnitude increase in the scale of the measured wavenumbers.
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 Xband 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.
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 stateoftheart.
The method is easy to implement. It takes the present stateofthe 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 stateoftheart 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 testbed for further studies regarding remote sensing of nearsurface 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 improvement in accuracy relative to the EDM for profiles with strong nearsurface 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 using the second harmonic to increase sensitivity to currents.
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://dataverse.no/dataset.xhtml?persistentId=doi%3A10.18710%2F8JBWCJ&version=DRAFT.References
 Banihashemi et al. (2017) Banihashemi, S., Kirby, J. T., and Dong, Z. (2012), Approximation of Wave Action Flux Velocity in Strongly Sheared Mean Flows, Ocean Modelling 116, 33–47.
 Campana et al. (2016) Campana, J., Terrill, E. J., and de Paolo, T. (2017), The development of an inversion technique to extract vertical current profiles from Xband radar observations. J. Atmos. Oceanic Technol., 33,20152028.
 Campana et al. (2017) Campana, J., Terrill, E. J., and de Paolo, T. (2017), A New Inversion Method to Obtain Upper–Ocean Current–Depth Profiles Using XBand Observations of Deep–Water Waves. J. Atmos. Oceanic Technol., 34, 957–970.
 Crombie (1955) Crombie D. D. (1955), Doppler spectrum of sea echo at 13.56 Mc./s. Nature 175, 681–682.
 Dalrymple (1973) Dalrymple, R. A. (1973), Water wave models and wave forces with shear currents, Tech. Rep. 20, Coastal and Oceanographic Engineering Laboratory, Uni. Florida.
 Dugan & Piotrowski (2003) Dugan, J.P. & Piotrowski, C.C. (2003), Surface current measurements using airborne visible image time series. Remote Sens. Environ., 84, 309319.
 Dugan et al. (2001) Dugan, J.P., Piotrowski, C.C. & Williams, J.Z. (2001), Water depth and surface current retriavals from airborne optical measurements of surface gravity wave dispersion. J. Geophys. Res., 106, 1690316915.
 Dunn & Tavoularis (2007) Dunn, W. & Tavoularis, S. (2007), The use of curved screens for generating uniform shear at low Reynolds numbers. Exp. Fluids, 42, 281290.
 Ellingsen & Li (2017) Ellingsen, S.Å. & Li, Y. (2017), Approximate dispersion relations for waves on arbitrary shear flows. J. Geophys. Res.: Oceans, 122 9889–9905.
 Fernandez et al. (1996) Fernandez, D.M., Vesecky, J.F. & Teague, C. (1996), Measurememts of upper ocean surface current shear with highfrequency radar. J. Geophys. Res., 101, 2861528625.
 Gangeskar (2002) Gangeskar, R. (2002), Ocean current estimated from Xband radar sea surface images. IEEE Trans. Geosci. Remote Sens. 40, 783–792.
 Ha (1979) Ha, E.C.. (1979), Remote sensing of ocean surface current and current shear by HF backscatter radar. Tech. Rep. D4151, Stanford University.
 Harper & Dixon (1974) Harper, J. F. & Dixon, J. N. (1974), The leading edge of a surface film on contaminated flowing water. Proc. 5th Australasian Conf. on Hydraulics and Fluid Mech., 499505.
 Hessner & Bell (2009) Hessner, K. & Bell, P. S. (2009), High resolution current and bathymetry determined by nautical Xband radar in shallow waters. OCEANS 2009EUROPE, 1–5.
 Hessner et al. (2014) Hessner, K., Reichert, K., Borge, J.C.N., Stevens, C. L., & Smith, M. J. (2014), Highresolution Xband radar measurements of currents, bathymetry and sea state in highly inhomogeneous coastal areas. Ocean Dyn., 64, 989998.
 Horstmann et al. (2017) Horstmann, J., Stresser, M., & Carrasco, R. (2017), Surface currents retrieved from airborne video. Proc. OCEANS 2017Aberdeen.
 Huang & Gill (2012) Huang, W. & Gill, E. (2012), Surface current measurement under low sea state using dual polarized Xband nautical radar. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., 5, 1868–1873.
 Kilcher & Nash (2010) Kilcher, L. F. & Nash, J. D. (2010), Structure and dynamics of the Columbia River tidal plume front. J. Geophys. Res.: Oceans, 115, C05590.
 Kirby & Chen (1989) Kirby, J. T. and Chen, T.M. (1989), Surface waves on vertically sheared flows: approximate dispersion relations. J. Geophys. Res. 94, 1013–1027.
 Kudryavtsev et al. (2008) Kudryavtsev, V., Shrira, V., Dulov, V. & Malinovsky, V. (2008), On the vertical structure of winddriven sea currents. J. Phys. Ocean. 38, 2121–2144.
 Laxague et al. (2017) Laxague, N. J. M., et al. (2017), Passive optical sensing of the nearsurface winddriven current profile. J. Atmos. Ocean. Technol. 34, 1097–1111.
 Laxague et al. (2018) Laxague, N. J. M., et al. (2018), Observations of nearsurface current shear help describe oceanic oil and plastic transport. Geophys. Rev. Lett. 45, 245249.
 Li & Ellingsen (2019) Li, Y. & Ellingsen, S.Å., A framework for modelling linear surface waves on arbitrary vertical shear currents and varying bathymetry. J. Geophys. Res.: Oceans (in press, DOI:10.1029/2018JC014390).
 Lund et al. (2015) Lund, B., Graber, H. C., Tamura, H., Collins III, C. O., and Varlamov, S. M. (2015), A new technique for the retrieaval of nearsurface vertical current shear from marine Xband radar images. J. Geophys. Res. 120, 8466–8486.
 Lund et al. (2018) Lund, B. et al. (2018), Nearsurface current mapping by shipboard marine Xband radar: a validation. J. Atmos. Ocean. Technol. 35, 1077–1090.
 Moisy et al. (2009) Moisy, F., Rabaid, M. & Salsac, K. (2009), A synthetic schlieren method for the measurement of the topography of a liquid surface. Exp. Fluids 36, 1021–1036.
 Plant & Wright (1980) Plant, W.J. & Wright, J.W. (1980), Phase speeds of upwind and downwind traveling short gravity waves. J. Geophys. Res. 85, 3304–3310.
 Scott (1982) Scott, J. C. (1982), Flow beneath a stagnant film on water: the Reynolds ridge. J. Fluid Mech., 116, 283–296.
 Senet et al. (2001) Senet, C. M., Seeman, J. & Ziemer, F. (2001), The nearsurface current velocity determined from image sequences of the sea surface. IEEE Trans. Geosci. Remote Sens. 39, 492–505.
 Serafino et al. (2010) Serafino, F., Lugni, C. & Soldovieri, F. (2010), A novel strategy for the surface current determination from marine Xband radar data. IEEE Geosci. Remote Sens. Lett. 7, 231–235.
 Shrira et al. (2001) Shrira, V., Ivonin, D. V., Broche, P. & Maistre, J. C. (2001), On remote sensing of vertical shear of ocean surface currents by means of a singlefrequency VHF radar. Geophys. Res. Lett. 28, 3955–3958.
 Skop (1987) Skop, R. A. (1987), Approximate dispersion relation for wavecurrent interactions. J. Waterway, Port, Coastal, and Ocean Eng. 113, 187–195.
 Stewart & Joy (1974) Stewart, R. J. & Joy, J. W. (1974), HF radio measurements of surface currents. Deep Sear Res. Oceanograph. Abs. 21, 10391049.
 Teague et al. (2001) Teague, C. C., Vescky, J. F. & Hallock, Z. R. (2001), a comparison of multifrequency HF radar and ADCP measurements of nearsurface currents during COPE3. IEEE J. Oceanic Eng. 26, 399–405.
 Terray et al. (1996) Terray, E. A., et al., (1996), Estimates of kinetic energy dissipation under breaking waves. J. Phys. Oceanogr., 26, 792–807.
 Willert et al. (2010) Willert, C., Stasicki, B., Klinner, J. & Moessner, S. (2010), Pulsed operation of highpower light emitting diodes for imaging flow velocimetry. Meas. Sci. Techn., 21, 075402.
 Young & Rosenthal (1985) Young, I. R. & Rosenthal, W. (1985), A threedimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents. J. Geophys. Res. 90, 1049–1059.
 Zippel & Thomson (2017) Zippel, S. and Thomson, J. (2017), Surface wave breaking over sheared currents: Observations from the Mouth of the Columbia River.J. Geophys. Res.: Oceans 122(4): 3311–3328 .