Coastal Engineering

A novel empirical model for the streamwise velocity in oscillatory boundary layer flow, valid in the rough turbulent regime, is presented. The model consists of simple expressions that require only the free-stream velocity time-series and equivalent sand-grain roughness length of the bed to be known a priori . A frequency-independent attenuation and phase lead is assumed for all flow harmonics, expressions for which are extracted from data from previous experimental studies made in 3 different oscillatory flow tunnels. Only the oscillating component of the flow is considered in the model; steady streaming is neglected. Errors in kinematics predicted by the model are typically two orders of magnitude smaller than the maximum oscillatory velocity in the free-stream. Hence, it is well-suited to engineering application due to its simplicity and accuracy.


Introduction
Over the past few decades, considerable research effort has been directed toward the study of oscillatory boundary layer flows similar to those that occur at the seabed under sea waves in the coastal zone.Using dimensional analysis, Jonsson (1966) demonstrated that the boundary layer thickness  bl and friction factor  w are functions of the Reynolds number  =  ∕ and inverse relative roughness ∕  , where  is the maximum oscillatory velocity in the free-stream,  =  ∕ is the orbital semi-excursion,  = 2∕ is the flow angular frequency,  is the flow period,  is the fluid kinematic viscosity and   is the equivalent sand-grain roughness length (Nikuradse, 1933).Several laboratory studies have been previously conducted to investigate oscillatory boundary layer hydrodynamics; key studies are indicated in the -∕  regime diagram shown in Fig. 1.Most early experiments involved sinusoidal oscillatory flows.In the natural coastal environment, this assumption is simplistic because waves, and hence, the free-stream velocity time-series near the seabed, develop significant skewness and asymmetry as a result of shoaling processes (e.g.Malarkey and Davies, 2012).Many studies (e.g.King, 1991;Dibajnia and Watanabe, 1998;van der A et al., 2010;Silva et al., 2011) have demonstrated the importance of flow skewness and asymmetry for the transport of sediments.These sediment transport studies have been complemented by more recent research in which the effects of skewness and asymmetry on rough-wall boundary layer hydrodynamics have been experimentally investigated in oscillatory flow tunnel (OFT) facilities (van der A et al., 2011;Yuan and Madsen, 2014;O'Donoghue et al., 2021;Dunbar, 2022).An OFT is a U-tube facility in which a piston is used to generate an oscillatory flow.OFT facilities are well suited to the study of flows analogous to those that occur under real waves due to the capacity of a relatively small facility to generate large amplitude oscillatory flows.Additionally, the amplitude and phase of each harmonic component of the flow in an OFT can be precisely controlled.However, unlike the flow under surface gravity waves, the mean flow field in an OFT is horizontally uniform and lacks a vertical component.Xie et al. (2021) conducted a study of flows with  and ∕  values similar to van der A et al. (2011), Yuan and Madsen (2014), O'Donoghue et al. (2021) and Dunbar (2022) in a large wave flume.In this type of facility, the flow is naturally comprised of multiple harmonics but their relative amplitude and phase cannot be controlled.
The studies of van der A et al. (2011), Yuan and Madsen (2014), O'Donoghue et al. (2021), Xie et al. (2021) and Dunbar (2022) primarily cover the rough turbulent regime (see Fig. 1).This regime has been a particular focus for recent study because it represents the bottom boundary layer under storm-generated waves over a seabed comprised of sand or gravel, which is an especially important scenario for the transport of sediments.
Some simpler models for oscillatory boundary layer flow have been proposed.Several authors have used the eddy viscosity hypothesis of Boussinesq (1877) as the basis for such a model.Taking this approach, the equation of motion in the boundary layer is given by where  is time,  is ensemble-averaged streamwise velocity, subscript ∞ denotes a quantity in the free-stream,  is the vertical coordinate and  t is the eddy viscosity.Kajiura (1968), Grant and Madsen (1979), Brevik (1981) and Myrhaug (1982) proposed models where  t depends on  only and provided solutions for sinusoidal oscillatory flow.Trowbridge and Madsen (1984a,b), Davies (1986) and Gonzalez-Rodriguez and Madsen (2011) extended the work to  t that depends on both  and  for flows consisting of up to two harmonics in the free-stream.Yuan and Madsen (2015) developed a model that considered a larger number of flow harmonics.While the eddy viscosity approach gives useful results, fundamentally, these models have the weakness that they rely on assumed  t that are not always consistent with experimental data.Several authors (e.g.Jonsson and Carlsen, 1976;Sleath, 1987;Nielsen, 1985) pointed out that the shear stress and velocity gradient ∕ are generally not equal to zero at the same phase, leading to  t values that are negative or approach infinity at certain phases of the flow cycle (Nielsen, 1992).This behaviour leads to questions regarding the applicability of the eddy viscosity concept to oscillatory flows.Another approach to modelling turbulent oscillatory boundary layers was developed by Nielsen (1985Nielsen ( , 1992Nielsen ( , 2016) ) and Teng et al. (2022).For the classical problem of laminar oscillatory flow over a smooth wall (e.g.Stokes, 1851), if the first harmonic of free-stream velocity is given by  1,∞ = Re {  1   } , where  1 is the first harmonic orbital amplitude, the first harmonic velocity in the boundary layer is given by , where the complex velocity defect  1 =  −(1+)∕ and  = √ 2∕ is the Stokes length.Nielsen (1985) noted that ln | 1 | = arg( 1 ) in this case and demonstrated that this identity also holds reasonably well for non-laminar experimental data from Kalkanis (1964), van Doorn (1981, 1982, 1983) and Sleath (1982), suggesting the general form of  1 is Nielsen (1985Nielsen ( , 2016) where with universal vertical scale According to Nielsen (2016), this model holds when  < 2.5 × 10 4 and 17 < ∕  < 100.When  > 1.3 × 10 5 and ∕  > 476, ln It is worth noting that ln | 1 | = arg( 1 ) implies that the bottom phase lead of the first harmonic,  0 , is equal to 45 • (Nielsen, 1985).However, when  > 1.3 × 10 5 ,  0 typically takes on values of 15 • -30 • (e.g.van der A et al., 2011;Yuan and Madsen, 2014) which may explain why ln | 1 | ≠ arg( 1 ) in this range of Reynolds number.
Recently, Teng et al. (2022) proposed an improved version of the Nielsen (1985Nielsen ( , 1992Nielsen ( , 2016) ) model, applicable to a wider range of  and ∕  .They proposed a modified defect function that takes the form allowing flows for which ln | 1 | ≠ arg( 1 ) to be modelled, and provided empirical expressions for  1 ,  2 ,  1 and  2 based on experimental and numerical data across a wide range of  and ∕  .Their model generally performed well across a range of flow regimes for approximation of the first harmonic of streamwise velocity, in some cases outperforming a far more computationally intensive  −  turbulence model.However, the model of Teng et al. (2022) only predicts the first harmonic component of velocity.In practice, oscillatory flows under real waves have significant skewness and asymmetry, and multiple harmonic components are necessary to capture the flow shape.This article presents a new practical model for the kinematics of oscillatory flow in the rough turbulent regime that is thoroughly validated against experimental data.The model is empirically based on experimental data from several previous studies and is applicable to flows comprised of multiple harmonic components.The article is organised as follows: Section 2 introduces the experimental datasets used for model calibration.The model is presented in Section 3. In Section 4, the model is validated using experimental data.Finally, Section 5 contains concluding remarks.

Experimental datasets
In this article, experimental data from several previous OFT experiments (Jonsson and Carlsen, 1976;van der A et al., 2011;Yuan and Madsen, 2014;O'Donoghue et al., 2021;Dunbar, 2022) are used to develop and validate the empirical model.In the experiments, velocity profiles in a periodic oscillatory boundary layer over a fixed rough bed were measured using particle image velocimetry (PIV), laser Doppler anemometry (LDA), or a micropropellor.The data cover a range of  and ∕  , primarily within the rough turbulent regime, and include flows that consist of multiple harmonic components.The data come from measurements made in 3 different OFT facilities.A summary of the flow conditions is presented in Table 1.Note that the values of  and  reported in Table 1 may differ slightly to those given in the respective sources because of differences in the definitions used.  values in the table are the values reported in the respective sources.The Free-stream skewness  ∞ and asymmetry  ∞ are computed using and respectively, where  p =  −  is the oscillating component of velocity, overbar denotes a time average and  is the Hilbert transform.

Model formulation
The present work considers the problem of oscillating flow driven by an unsteady periodic pressure gradient in the vicinity of a rough boundary.A one-dimensional ensemble-averaged flow field (, ) is assumed.The flow described is practically identical to the flow inside an OFT.However, the steady streaming that normally occurs in an OFT is neglected.The magnitude of this streaming is typically <3% of max( ∞ ) (e.g.Yuan and Madsen, 2014).Hence, only the oscillating component of (, ),  p (, ), is considered.
Using a Fourier decomposition,  p (, ) can be written where   is the complex Fourier coefficient at frequency , * denotes complex conjugation and  is the number of harmonics that comprise the flow.Eq. ( 8) can be written in terms of real numbers as Hence, the oscillating component of the free-stream velocity time-series with zero mean current can be written as where are the amplitude and phase of the th harmonic of free-stream velocity, respectively.
Noting that the flow harmonics undergo an attenuation and phase shift within the boundary layer relative to the free-stream, the corresponding oscillatory velocity in the boundary layer can be written where and are the attenuation and phase lead of harmonic component  relative to the free-stream.In the classical smooth-wall laminar case, the exact analytical solutions for   and   are (Svendsen, 2006) and respectively.In the case of rough turbulent flow,   and   cannot be derived analytically due to the turbulence closure problem.Instead, empirical expressions for   () and   () can be extracted from experimental data.This approach is taken in the present work.

Assumption of frequency independence of 𝐾 𝑛 and 𝜙 𝑛
For discretely sampled data, the bispectrum of a time-series for frequency components   and   is given by (see Kim and Powers, 1979, for a derivation) where  denotes an ensemble average.Elgar and Guza (1985) and Elgar (1987) noted that for discretely sampled data,  and  are related to the real and imaginary parts of (  ,   ) by and respectively, where 1 ≤  ≤ , 1 ≤  ≤ ,  > ,  +  ≤ .By assuming each harmonic component of an oscillatory flow is subject to a frequency-independent attenuation  and phase shift  within the boundary layer relative to the free-stream (i.e.  =    ,∞ for all ), the bispectrum within the boundary layer is given by Hence, the following equations for  and  can be derived following Henderson et al. (2004) and Berni et al. (2013): Numerical results from a RANS model (Henderson et al., 2004), as well as measurements from wave flume facilities (Berni et al., 2013;Henriquez et al., 2014;Fromant et al., 2018) have shown a relationship between  at a single vertical position, , very close to the bed and  ∞ that agrees fairly well with Eq. ( 20) for transitional and turbulent flows, taking  as the phase lead of the first flow harmonic as close to  = 0 as possible.Fig. 2 compares vertical profiles of  and  from the nonsinusoidal experimental cases of Dunbar (2022) with  and  estimated using Eqs.( 20) and ( 21), taking  as  1 () extracted from the experimental measurements.The figure shows that there is good agreement between the data and Eqs. ( 20) and ( 21) throughout the boundary layer, not just at  ≈ 0. This suggests that the assumption of a frequency-independent  and  allows for accurate estimation of the shape of the oscillatory velocity time-series throughout the boundary layer.This assumption is adopted for the present practical model, taking  =  1 and  =  1 , to reduce the number of empirical expressions for   and   necessary for closure of the model from  each to just one each.Hence, Eq. ( 11) is simplified to 3.2.Empirical expressions for  1 () and  1 () Fig. 3 shows vertical profiles of  1 and  1 , obtained from all experimental flow cases in Table 1.Note that  1 ̸ → 0 as  → 0. This is due to the definition that  = 0 corresponds to a representative roughness crest level, and there is some flow below this level.This definition is adopted because of its practicality, since for a known roughness topology, the vertical position of the theoretical bottom location obtained from a logarithmically fitted velocity profile is not known a priori and would Fig. 3. Vertical profiles of  1 and  1 from all experimental cases in Table 1.
require an additional empirical equation, whereas a suitably defined representative roughness crest level can be obtained directly from the known roughness topology.Fig. 4 shows vertical profiles of  1 and  1 ∕ 0 where the vertical coordinate has been normalised by boundary layer thickness  bl as defined by Jensen et al. (1989), i.e. the  position at which  = max() at the phase when  ∞ = max( ∞ ).This definition of  bl has been shown to have a clear relationship with ∕  of the form  bl ∕  =  ( ∕  )  , where  and  are empirical constants, by many previous authors (e.g.Jonsson, 1966;Fredsøe and Deigaard, 1992;van der A et al., 2011;Yuan and Madsen, 2014).The plots show that each set of profiles collapse reasonably well onto one line when this normalisation is applied.There is some variation in the plot of  1 ∕ 0 .This is likely due to the uncertainty in determining the value of bottom phase lead,  0 , because LDA and PIV measurements are less reliable in the vicinity of a boundary as a result of laser reflection, and the precise definition of  = 0 differs slightly between the respective experimental sources.
Also shown on each plot is a line fitted by least-squares regression to the mean of all the experimental profiles linearly interpolated onto an evenly spaced  array.Expressions for  1 and  1 obtained from this fitting procedure are given by  1 = { 0.98 ŷ3 −0.77 ŷ2 +0.57ŷ+0.0079 ŷ3 −0.87 ŷ2 +0.58 ŷ+0.033 and where ŷ = ∕ bl .It is assumed that ŷ > 5 is in the free-stream; hence  1 = 1 and  1 ∕ 0 = 0 at ŷ > 5. where  c = 2 ac ∕ c following Silva et al. (2006) and van der A et al. (2011), where  c and  ac are the times from the zero-up crossing of  ∞ () to the zero-down crossing of  ∞ () and the time that  ∞ () = max[ ∞ ()], respectively.The bottom phase lead  0 can be estimated using the expression of Humbyrd (2012), which for rough turbulent flow is as follows ( 0 in degrees): Eqs. ( 22)-( 26) form a complete model for  p (, ) in a rough turbulent oscillatory boundary layer for given   and  ∞ .For flow over roughness with unknown   ,   can be estimated using   ≈ 2 50 , where  50 is median grain diameter, for roughness composed of fixed sand or gravel particles, or the expression of Flack et al. (2020), where   and   are the standard deviation and skewness of the roughness height field, respectively, can be used for roughness with known height field statistics.

Model validation
To test the performance of the model, a comparison of  p from the experimental data shown in Table 1 is made with  p computed using Eqs.( 22)-( 26).For each case shown in Table 1, input to the model is the measured  ∞ and the   reported by the authors for their respective experiments.For the comparison,  in Eq. ( 22) is set to 6.This value is excessive for the near-sinusoidal flows but is chosen to ensure that  p,∞ () is reproduced with sufficient accuracy for all flow cases.Fig. 5 compares normalised vertical profiles of  p interpolated onto 8 phases evenly distributed through the flow cycle.The figure shows that there is generally excellent agreement between the experimental data and the model, though in some cases a discrepancy can be observed.In all cases with significant  ∞ , such as ST200_ce and SK1549, discrepancies between the experimental and predicted profiles are apparent.This may suggest that the assumption of a frequencyindependent  and , which only becomes relevant for non-sinusoidal flows, could be a significant source of error in the model.However, the cases with significant  ∞ often exhibit very good agreement.For example, predicted and measured profiles in cases S605010c and FL320a_sa agree well despite  ∞ values of 0.28 and 0.47, respectively.It is worth noting that the experimental data used for model calibration contains more asymmetric than skewed flow cases, and most of the skewed flows have  ∞ values near 0.5, compared with the asymmetric flows that have more varied  ∞ values.Additionally, predicted profiles for some flows with near-zero values of  ∞ and Fig. 7. Comparison of experimental and predicted  p ( ′ ,  ′ ) for all flow cases in Table 1 at all  ′ and  ′ (dots).Both axes are normalised by the experimental value of  .The solid and dashed red lines correspond to perfect agreement and ±0.1 discrepancy, respectively.
∞ , such as cases test no. 2 and SP200_ce, also deviate noticeably from the measured profiles despite not being affected by the assumption of frequency-independent  and .Hence, the assumption of frequencyindependent  and  is not necessarily a primary source of error in the model.
Fig. 6 compares predicted and experimental normalised time-series of  p at ŷ = 1 and at a vertical position from the respective experimental dataset as close to ŷ = 0.1 as possible for 8 selected flow cases.The left and right columns show the comparison for near-sinusoidal and significantly non-sinusoidal flows, respectively.The 8 selected cases include flows for which the agreement between experiment and model seen in Fig. 5 is very good (S505010c, CSSR, S605010c, FL320a_sa) and for which a clear discrepancy can be observed (test no. 2, SP200_ce, ST200_ce, SK1549).Even for cases for which discrepancies are visible, there is still good agreement between the predicted and experimental time-series.
To undertake a more quantitative assessment of the overall model performance across all flow cases shown in Table 1, experimental and predicted  p values are compared.To ensure an equal number of comparisons are made for each flow case from Table 1, distributed evenly through the flow cycle and boundary layer,  p (, ) is linearly interpolated onto new vertical and temporal coordinates, 0.05 bl <  ′ ≤ min(max[], 5 bl ), and 0 ≤  ′ <  , with constant spacing,  ′ = min(max[], 5 bl )∕25 and  ′ =  ∕32, respectively.The lower limit of  = 0.05 bl is chosen to remove a small number of outliers that occur very close to  = 0 that result from slight differences in the definitions and methodology used to determine  = 0 between the respective experimental sources, and due to reduced LDA and PIV measurement accuracy in the vicinity of the bed; the upper limit of  = 5 bl is selected because above this level  1 = 1 and  1 = 0, so experimental and predicted values of  p are identical if  is sufficiently large.Fig. 7 shows the comparison of measured and predicted  p ( ′ ,  ′ ) for all the flows in Table 1.
The figure shows that there is excellent overall agreement for these flow conditions, with nearly all predictions following the 1:1 line closely.
The normalised error between predicted and experimental kinematics can be computed using  28) for all data points shown in Fig. 7.
that 25%, 50% and 75% of the values of  p ( ′ ,  ′ ) have a normalised error of < 0.01, < 0.02 and < 0.04, respectively.Hence, the model performs very well overall in predicting  p for a wide range of flow conditions primarily within the rough turbulent regime, throughout the boundary layer and flow cycle, with errors typically about two orders of magnitude smaller than  .
To obtain a more detailed insight into the model performance, the mean normalised error NE  is computed for each flow case using where   ′ and   ′ are the total number of  ′ and  ′ elements, respectively.Fig. 9

Conclusions
A novel empirical model for the oscillating streamwise velocity in a rough turbulent oscillatory boundary layer has been presented.A summary of key features of the model is as follows.
The model is based around the assumption of a frequency-independent attenuation  1 and phase lead  1 of each flow harmonic, an assumption that is demonstrated to be capable of accurately approximating the oscillating flow from experimental measurements.Empirical expressions for  1 and  1 are extracted from experimental data from  previous laboratory studies conducted in 3 distinct large-scale OFT facilities.The dataset covers 43 flow conditions, including flows comprised of multiple harmonic components.While the experimental data considered here are concentrated in the rough turbulent regime, in principle a similar empirical approach could be applied to experimental data from a range of flow regimes to obtain a more widely applicable model.
The model is simple to implement, requiring only the free-stream velocity time-series  ∞ () and equivalent sand-grain roughness length   as inputs.Orbital amplitude  c and flow frequency  can be calculated from  ∞ ().Amplitudes and phases of the harmonic components of  ∞ (),   and   , and orbital amplitude  1 , can be obtained using the fast Fourier transform (FFT) algorithm. 1 ,  1 and  p (, ) can then be estimated by substitution into Eqs.( 22)-( 26).
The model has been shown to perform very well in predicting  p (, ), with errors generally two orders of magnitude smaller than  .Additionally, the mean normalised error NE  of model predictions does not depend on , ∕  or | ∞ |.It is possible that flows with larger | ∞ | result in slightly reduced model accuracy, but further investigation is necessary to confirm if this is the case.Overall, the model is well-suited to engineering application due to its simplicity and accuracy.
The present model utilises data from OFT facilities to derive its empirical expressions.Since the flow in an OFT very closely approximates the bottom boundary layer under surface gravity waves, the model is likely to also allow accurate estimation of wave boundary layer kinematics; additionally, the small steady streaming that occurs in OFT facilities is neglected.The steady streaming under a surface gravity wave includes an additional contribution due to the streamwise inhomogeneity of the ensemble-averaged flow (e.g.Longuet-Higgins, 1953).The present model could, in principle, be linearly superposed with a model for the vertical profile of steady streaming to form a complete model for the ensemble-averaged flow in a rough turbulent oscillatory boundary layer under surface gravity waves.More experimental data are needed to test and develop the model for application to combined wave-current flow.
The new model is well-suited to forming the basis for an improved sediment transport model for rough turbulent oscillatory flow, since it offers practical predictions for the kinematics within the boundary layer which are generally not considered in existing sediment transport models.For this purpose, a prediction of intra-period bed shear stress,  0 (), which is not included in the present model, would also be necessary.An engineering predictor for  0 () in an oscillatory boundary layer such as that proposed by Larsen and Fuhrman (2019), which does not require

Fig. 6 .
Fig. 6.Comparison of normalised time-series of  p at ŷ = 1 (experimental: solid black line, predicted: dashed red line) and at the vertical position closest to ŷ = 0.1 (experimental: dash-dotted black line, predicted: dotted red line).

Fig. 8 .
Fig. 8 shows the cumulative distribution function (CDF) of NE for all flow cases shown inTable 1 at all  ′ and  ′ .The figure demonstrates

Fig. 9 .
Fig. 9. Scatter plots of NE  versus  (top-left), ∕  (top-right), | ∞ | (bottom-left) and | ∞ | (bottom-right) for all flow cases shown in Table1.Also shown are Pearson correlation coefficients  P between NE  and the -variable of each plot.
data from each flow condition consist of ensemble-averaged velocity (, ), with  = 0 corresponding to a representative roughness crest level (as defined in the respective sources), 0 ≤  <  and  = 0 corresponding to the phase of the zero-up crossing of  ∞ = (max[], ).
The boundary layer thickness  bl can be estimated using(van der A  et al., 2011)Comparison of normalised vertical profiles of  p from experimental data (solid black lines) and predicted by Eqs.(22)-(26) (dashed red lines) at 8 phases evenly distributed through the flow cycle for all flow conditions shown in Table1.Both the experimental and predicted values of  p are normalised by the experimental value of  and the vertical coordinates are normalised by  bl from experimental data.
shows NE  plotted against , ∕  , | ∞ | and | ∞ | to check if the performance is affected by any of these parameters.The Pearson correlation coefficients  P between NE  and the -variable of each plot are also shown.The data and  P values show that there is no significant correlation between NE  and , ∕  or | ∞ |.Hence, changing these parameters within the range of the data used for model calibration does not affect the reliability of the model predictions.The value of  P = 0.65 suggests that there may be a positive correlation between NE  and | ∞ |.The most obvious reason for model prediction accuracy to be affected by | ∞ | is the error introduced by the assumption of a frequency-independent  and  applied to all flow harmonics, since the influence of higher harmonics becomes important for flows with larger | ∞ |.However, if this is the case, it would be expected that a correlation between NE  and | ∞ | would also exist, which is not seen in Fig.9.Hence, it is plausible that the value  P = 0.65 seen in the figure is not statistically significant, but results from the limited number of flows with significant | ∞ | used for validation, the majority of which have | ∞ | values very close to 0.5.A larger number of flows with significant | ∞ | taking on a more varied range of values would be necessary to determine with more certainty whether | ∞ | truly affects the reliability of model predictions.

Table 1 .
Also shown are Pearson correlation coefficients  P between NE  and the -variable of each plot.