Next Article in Journal
Verification and Validation of CFD Based Form Factors as a Combined CFD/EFD Method
Previous Article in Journal
Social, Economic and Environmental Sustainability of Port Regions: MCDM Approach in Composite Index Creation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development and Validation of Quasi-Eulerian Mean Three-Dimensional Equations of Motion Using the Generalized Lagrangian Mean Method

by
Duoc Tan Nguyen
1,2,4,*,
Niels G. Jacobsen
3 and
Dano Roelvink
1,2,3
1
IHE Delft Institute for Water Education, Westvest 7, 2611 AX Delft, The Netherlands
2
Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628 CN Delft, The Netherlands
3
Deltares, Boussinesqweg 1, 2629 HV Delft, The Netherlands
4
Vietnam Institute of Seas and Islands, Hanoi 100000, Vietnam
*
Author to whom correspondence should be addressed.
J. Mar. Sci. Eng. 2021, 9(1), 76; https://doi.org/10.3390/jmse9010076
Submission received: 12 December 2020 / Revised: 6 January 2021 / Accepted: 8 January 2021 / Published: 13 January 2021
(This article belongs to the Section Ocean Engineering)

Abstract

:
This study aims at developing a new set of equations of mean motion in the presence of surface waves, which is practically applicable from deep water to the coastal zone, estuaries, and outflow areas. The generalized Lagrangian mean (GLM) method is employed to derive a set of quasi-Eulerian mean three-dimensional equations of motion, where effects of the waves are included through source terms. The obtained equations are expressed to the second-order of wave amplitude. Whereas the classical Eulerian-mean equations of motion are only applicable below the wave trough, the new equations are valid until the mean water surface even in the presence of finite-amplitude surface waves. A two-dimensional numerical model (2DV model) is developed to validate the new set of equations of motion. The 2DV model passes the test of steady monochromatic waves propagating over a slope without dissipation (adiabatic condition). This is a primary test for equations of mean motion with a known analytical solution. In addition to this, experimental data for the interaction between random waves and a mean current in both non-breaking and breaking waves are employed to validate the 2DV model. As shown by this successful implementation and validation, the implementation of these equations in any 3D model code is straightforward and may be expected to provide consistent results from deep water to the surf zone, under both weak and strong ambient currents.

1. Introduction

The interaction between waves and currents has been the subject of much research in recent decades. There are two representations of wave-averaged effects on the currents called “radiation stress” and “vortex force”. The concept of “radiation stress: was first introduced by Longuet-Higgins and Stewart [1] to explain the transfer of wave energy to a uniform current. This concept was used by Longuet-Higgins and Stewart [2] to study the changes in the mean surface level and the currents caused by gravity waves. The radiation stress concept has been successful in explaining phenomena such as wave “set-up”, “surf beats”, the steepening of the surface waves on adverse currents [3], and the generation of long-shore currents by oblique incident waves [4,5,6]. However, since “radiation stress” introduced by Longuet-Higgins and Stewart [1] is a two-dimensional horizontal tensor it is only practical for two-dimensional, depth-averaged models. In reality, the current is depth-dependent, so the vertical structure of the radiation stress should be specified.
Some scientists attempted to apply the “radiation stress” concept in three-dimensional models. Xie, Wu [7] applied radiation stress as a depth uniform body force in the Princeton ocean model, even though the radiation stress is caused by depth-varying wave velocity and hydrodynamic pressure. Therefore, their assumption for the vertical structure of radiation stress is not accurate. Xia, Xia [8] considered the vertical structure of radiation stress; however, the three-dimensional radiation stress formulation was derived from two-dimensional radiation stress.
The second representation of the wave and current interaction is expressed in terms of the vortex force. This representation was first developed by Craik and Leibovich [9] in the work of constructing a realistic theoretical model of steady Langmuir circulations. Their research focused on the near-surface layer to explain the generation of Langmuir currents as a result of the interaction between surface waves and wind-driven circulation through the action of a vortex force. Leibovich [10] extended this theory to allow vertical density stratification and slow time variation. McWilliams and Restrepo [11] developed a perturbation theory to obtain wave-averaged equations of motion. Their theory is based on the assumption of small wave slope and deep water. McWilliams, Restrepo [12] developed a system of mean equations of motion based on an asymptotic theory to account for the interaction of waves and currents. In this, the effects of waves on the current are expressed in terms of the vortex force formalism. However, the equations of McWilliams, Restrepo [12] are only valid when the ratio of mean current to the wave orbital velocity is a small quantity, and are only applicable outside the breaking zone. Newberger and Allen [13] developed a three-dimensional, hydrostatic model for surf zone applications, with applicability to linear waves interacting with a depth-uniform mean current. They divided the effect of waves on the mean currents into surface and body forces. The surface force represents the wave dissipation, and the body force represents the gradient of the Bernoulli head and vortex-force. The equations of McWilliams, Restrepo [12] were used by Uchiyama, McWilliams [14] for surf zone applications. In this, the non-conservative forcing by breaking waves, roller waves, bottom and surface streaming and wave-enhanced mixing are added through empirical formulas. Their equations were implemented in the COAWST (coupled ocean—atmosphere—wave—sediment transport) modeling system by Kumar, Voulgaris [15] with some modifications for empirical formulas of wave-induced forcing.
The relationship between “radiation stress” and “vortex force” representations was studied by Lane, Restrepo [16]. In this, the asymptotic assumption proposed by McWilliams, Restrepo [12] was used to look for the similarities and discrepancies of these two representations. They proved that these two representations are equivalent (Equations (38) and (39)). However, their work was only restricted to non-dissipative waves. All the theories mentioned above are expressed in an Eulerian-mean framework, though when finite-amplitude waves are present, the region between the wave trough and wave crest is not always filled by the fluid but by the air during part of the wave period. This poses a problem due to a large difference in density between the fluid and the air.
In the work of Mellor [17] and Mellor [18], a wave-following sigma-coordinate system was employed to couple the three-dimensional circulation models with wave models. The coupling included depth-dependent wave radiation stress terms. Their equations are inconsistent in the simple case of shoaling waves without energy dissipation [19]. Recently, Mellor [20] and Mellor [21] derived prognostic equations for Eulerian mean flow on sigma-coordinates. The three-dimensional momentum equations were inferred from the vertically integrated momentum equations by adding a term for which vertical integration is zero. Similar to the work of Xia, Xia [8], the inference of three-dimensional momentum equations from two-dimensional momentum equations is not straightforward. This inconsistency was also pointed out by Ardhuin, Suzuki [22]. Moreover, in the momentum equation of Mellor [20] and Mellor [21], there is a lack of a term related to the divergence of the vertical momentum flux [22].
The generalized Lagrangian mean (GLM) method was introduced by Andrews and McIntyre [23], hereafter referred to as AM. The basic idea of this method is to average over disturbance positions of the fluid particle. Therefore, the GLM method is valid from the bottom to the mean water surface even in conditions of finite-amplitude waves. This method provides a powerful foundation for the analysis of the wave–current interaction and gives a physical interpretation of the interaction between waves and currents. Based on the GLM method, AM developed a set of equations of mean motion in a general condition of the wave–current interaction. Their set of equations is complete and depends on thermodynamic properties such as entropy and enthalpy. However, the disturbance-related quantities, which include both wave-induced and turbulence-induced effects, are not explicitly represented through source terms. Their equations were employed by Leibovich [24] to derive Langmuir circulation equations under the assumption that the waves are dominated by their irrotational part. The GLM equations of AM were simplified by Dingemans [25] with the assumption of constant density, and removing all thermodynamic terms, yet leaving the disturbance quantities as implicit. Groeneweg [26] used an alternative method to obtain GLM equations, where the Reynolds-averaged Navier–Stokes (RANS) equations were rewritten in terms of GLM quantities. The mean quantities in RANS equations are obtained by applying the Eulerian mean method; therefore, his set of equations is not suitable to the region above the wave trough. His set of equations was implemented in the Delft3D-FLOW model by Walstra, Roelvink [27] with simplification for the wave-induced driving force. Ardhuin, Rascle [19] developed a practical set of equations of mean motion based on the work of Dingemans [25]. Their equations are written in term of quasi-Eulerian mean velocity u ^ i defined by:
u ^ i = u ¯ i L p i ,
where, u ¯ i L is the ith- component of GLM velocity, and p i is the ith-component of pseudomomentum defined by:
p i = ξ j x i { u j l + ( Ω × ξ ) j } ¯ ,
where, x i is the ith-component of position x, ξ is the disturbance displacement of the fluid particle, u j l is the jth-component of Lagrangian disturbance velocity, and Ω is the angular velocity of the Earth. In Equation (2) the summation convention for the indices is employed. This convention is also used throughout this paper with the indices from 1 to 3.
The equations of mean motion by Ardhuin, Rascle [19] are explicit in terms of the wave forcing and applicable outside the breaking zone; their equations provide qualitative results for surf zone applications [19]. This is due to the fact that the Stokes drift only approximates to pseudomomentum when the waves are irrotational and the mean flow is of second-order of the disturbance amplitude (AM).
In this paper, a set of equations of mean motion using the GLM method is developed. The equations are written in terms of quasi-Eulerian velocity defined as GLM velocity minus Stokes drift. The new equations are valid from offshore to coastal areas. Outside the surf zone and for non-dissipative waves, the new equations are identical to equations of Ardhuin, Rascle [19]; for dissipative waves, there are subtle differences. In the case of infinitesimal and conservative waves, the new equations reduce to the well-known classical Eulerian mean equations of motion. The new set of equations is validated with an adiabatic test, non-breaking waves propagating on a strong ambient current in a wave flume, breaking waves propagating over a barred profile in a wave flume, and obliquely incident breaking waves in a large-scale sediment transport facility (LSTF).
Figure 1 is a flowchart describing the methodology of this research. It helps readers understand the structure of this paper easier.

2. Derivation of the Quasi-Eulerian Mean Equations of Motion

The GLM method is an exact theory of nonlinear waves on a Lagrangian-mean flow proposed by AM. In the following, only some properties of the GLM operator are present. Full details of this method are given in the original paper of AM. The basic idea of the GLM method is to average over positions displaced by a certain disturbance. That is, if ξ ( x , t ) is the particle displacement of the fluid particle then the GLM of any quantity φ ( x , t ) is defined as:
φ ( x , t ) ¯ L = φ { x + ξ ( x , t ) , t } ¯ ,
where, the operator ( ) ¯ expresses a time average over a wave period. The term on the right-hand side of the above equation is a usual Eulerian mean operator. The following notation was employed throughout this study:
φ ξ ( x , t ) = φ { x + ξ ( x , t ) , t } ,
The Lagrangian disturbance quantity φ l is defined by (AM):
φ l ( x , t ) = φ ξ ( x , t ) φ ¯ L ( x , t ) ,
The quasi-Eulerian mean quantity φ ¯ and quasi-Eulerian disturbance quantity φ are defined, respectively, by (AM):
φ ¯ ( x , t ) = φ ¯ L ( x , t ) φ ¯ S ( x , t ) ,
φ = φ φ ¯ ,
where, φ ¯ S is the Stokes correction of the quasi-Eulerian mean quantity φ ¯ .
Assuming that any quasi-Eulerian disturbance can be decomposed into wave and turbulence components, such as:
φ = φ ˜ + φ t ,
where, φ ˜ and φ t are the wave and turbulent quantities, respectively. Moreover, the turbulent and wave quantities are assumed to be uncorrelated. That is, for any φ and ψ :
φ ˜ ψ t ¯ = 0 ,

2.1. Derivation of Quasi-Eulerian Mean Equations of Motion

2.1.1. Derivation of Quasi-Eulerian Mean Momentum Equation

In this work, the fluid is assumed incompressible ( ρ = c o n s t ) and the dependence of hydrodynamic processes on thermodynamic terms is neglected (Assumption 1). Let us start with the momentum equation for the total flow written in the Eulerian framework. The ith equation is expressed by:
u i t + u j u i x j + 2 ( Ω × u ) i + Φ x i + 1 ρ p x i + X i = 0 ,
where, i and j represent for the spatial directions (i and j run from 1 to 3), the angular velocity of the Earth Ω is assumed constant, Φ ( x , t ) is the potential of the gravitational force, p is pressure, and X is a function of non-wave dissipative forcing.
Evaluating Equation (10) at the disturbance position of the fluid particle Ξ = x + ξ to obtain:
( u i t + u j u i x j ) ξ + ( 2 ( Ω × u ) i ) ξ + ( Φ x i ) ξ + ( 1 ρ p x i ) ξ + ( X i ) ξ = 0 ,
Equation (11) is valid from the bottom to the free water surface. Assuming that the gravitational acceleration g is constant then:
Φ x i = δ i 3 g ,
where, δ i 3 is the Kronecker delta function given by:
δ i j = { 1   if   i = j 0   otherwise ,
Taking the time-average of Equation (11) (and after some manipulation) to obtain the following momentum equation written in term of GLM quantities (see Appendix A for detail):
D ¯ L u ¯ i L + 2 ( Ω × u ¯ L ) i + δ i 3 g + 1 ρ p ¯ x i + X ¯ i L = 1 ρ x j ( ξ j p x i ) ¯ 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k + O ( ε 3 ) ,
where, ε = | ξ | is a small parameter in the order of disturbance displacement amplitude, and D ¯ L is the Lagrangian mean material derivative defined as:
D ¯ L = / t + u ¯ L . ,
where, = ( x , y , z ) is the gradient operator.
In the above equation, both wave and current-induced turbulent effects are involved in the quasi-Eulerian disturbance and GLM terms. For example, quasi-Eulerian disturbance pressure p includes wave-induced pressure p ˜ and turbulence-induced pressure p t . There is no available theory to calculate such quasi-Eulerian disturbance terms. Therefore, it is necessary to separate wave and turbulent terms from the quasi-Eulerian disturbance. In the following, Equation (14) is used to develop a quasi-Eulerian mean momentum equation in the GLM framework. The goal of this exercise is that the wave-induced velocity, the turbulence, and the mean velocity are separated.
Using the definition of quasi-Eulerian mean quantity in Equation (6), Equation (14) can be rewritten as:
D ¯ L u ¯ i + 2 ( Ω × u ¯ ) i + δ i 3 g + 1 ρ p ¯ x i + X ¯ i = D ¯ L u ¯ i S 2 ( Ω × u ¯ S ) i X ¯ i S 1 ρ x j ( ξ j p x i ) ¯ 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k + O ( ε 3 ) ,
After some manipulation, the right-hand side of Equation (16) can be expressed as (see Appendix B):
D ¯ L u ¯ i S 2 ( Ω × u ¯ S ) i X ¯ i S 1 ρ x j ( ξ j p x i ¯ ) 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k = ( u i u j ¯ ) x j + u ¯ k S u ¯ i x k + O ( ε 3 ) ,
The first term on the right-hand side of Equation (17) can be decomposed into the wave and turbulent components, such as:
( u i u j ¯ ) x j = ( u ˜ i u ˜ j ¯ ) x j ( u i t u j t ¯ ) x j ,
where, u ˜ i and u i t are ith—components of wave and turbulent velocities, respectively. It is stressed that an assumption of no correlation between wave and turbulent quantities is employed in this step (see Equation (9)).
Equation (17) expresses the relationship between the Stokes drift, wave-induced pressure, and wave-radiation stress. Inserting Equations (17) and (18) into Equation (16) to obtain the following momentum equation in terms of quasi-Eulerian mean velocity:
D u ¯ i D t + 2 ( Ω × u ¯ ) i + δ i 3 g + X i = 1 ρ p ¯ x i ( u ˜ i u ˜ j ¯ ) x j + 1 ρ τ ¯ i j x j + O ( ε 3 ) ,
where, τ ¯ i j = ρ u i t u j t ¯ is the turbulent stress tensor.
The momentum Equation (19) includes four variables: three components of quasi-Eulerian mean velocity u ¯ and pressure p ¯ . Effects of waves and turbulence can be modeled as source terms. The momentum equation can be solved in combination with the continuity equation. The turbulent effect is expressed in the form of Reynolds turbulent stress, which can be calculated by using existing turbulent submodels, and the wave effect following an appropriate wave model.

2.1.2. Mass Conservation Equation

The mass conservation equation is necessary to close the set of equations of the mean motion. The mass conservation is simplified under the assumption of slow modulation of the waves (Assumption 2). For an incompressible fluid (Assumption 1), the mass conservation equation is expressed by (AM):
u ¯ L x + v ¯ L y + w ¯ L z = 1 2 ( 3 ξ j ξ k ¯ t x j x k + u ¯ l L 3 ξ j ξ k ¯ x l x j x k ) ,
In the following, the vertical GLM velocity is assumed a small quantity, i.e., w ¯ L = O ( ε ) . Using Assumption 2 the right-hand side of Equation (20) is simplified as:
1 2 ( 3 ξ j ξ k ¯ t x j x k + u ¯ l L 3 ξ j ξ k ¯ x l x j x k ) = t ( 1 2 2 ξ 3 2 ¯ z 2 ) + O ( ε 3 ) ,
Notice that quasi-Eulerian mean quantities are averaged over the wave period. Therefore, in the scale of the mean current, the right-hand side of Equation (21), generally, differs from zero. According to the definition of Stokes correction in AM, the term in brackets on the right-hand side of Equation (21) is the Stokes correction of the mean position of the fluid particle Z ¯ S , i.e.,
Z ¯ S = 1 2 2 ξ 3 2 ¯ z 2 + O ( ε 3 ) ,
In stationary waves, the temporal derivative of Z ¯ S is zero, so Equation (20) becomes:
u ¯ x + v ¯ y + w ¯ z = ( u ¯ S x + v ¯ S y + w ¯ S z ) ,
However, in nonstationary waves, the right-hand side of Equation (21) differs from zero and is also of second-order of the disturbance amplitude. In general, the continuity Equation (20) is rewritten as:
u ¯ x + v ¯ y + w ¯ z = Z ¯ S t ( u ¯ S x + v ¯ S y + w ¯ S z ) ,
Equation (24) indicates that the divergence of quasi-Eulerian mean velocity is compensated by the divergence of Stokes drift and the time variation of Stokes correction of the mean position of the fluid particle. Combined with the momentum Equation (19) we had a set of four independent equations in four unknowns as long as the wave and turbulent motions are described by an appropriate wave theory and a relation to the mean flow, respectively. In principle, these equations can be solved numerically.

3. Validation of Quasi-Eulerian Mean Equations of Motion

3.1. Model Implementation

In this part, a two-dimensional numerical model is developed based on the quasi-Eulerian mean equations of motion, which were developed in the previous part. The model is written for the variation of hydrodynamic properties in the x- and z-directions (2DV model). All hydrodynamic properties are assumed to be uniform in the y-direction. The accelerations of mean vertical velocity and dissipative forcing are neglected in the vertical momentum equation of mean motion (hydrostatic assumption for the mean flow). The wave properties such as wave amplitude, wave energy, and wave energy dissipation are calculated by wave energy balance equation and roller energy equation given by Roelvink and Reniers [28]. Besides, the horizontal variation of the mean atmospheric pressure at the water surface and Coriolis’ effect are assumed small and neglected.

3.1.1. DV Governing Equations

The quasi-Eulerian mean momentum equations in the 2DV model are given by:
u ¯ t + u ¯ u ¯ x + w ¯ u ¯ z υ m o l Δ u ¯ = g ζ ¯ x ( u ˜ 2 ¯ w ˜ 2 ¯ ) x u ˜ w ˜ ¯ z + 1 ρ ( τ ¯ 11 x + τ ¯ 13 z ) + O ( ε 3 ) ,
v ¯ t + u ¯ v ¯ x + w ¯ v ¯ z + X 2 = υ m o l Δ v ¯ = u ˜ v ˜ ¯ x v ˜ w ˜ ¯ z + 1 ρ ( τ ¯ 21 x + τ ¯ 23 z ) + O ( ε 3 ) ,
where, only molecular viscosity is considered as non-wave dissipative forcing and υ m o l is molecular viscosity.
The components of turbulence stress tensor are parameterized by:
τ ¯ 11 = ρ υ h T u ¯ x , τ ¯ 13 = ρ υ v T u ¯ z ,
τ ¯ 21 = ρ υ h T v ¯ x , τ ¯ 23 = ρ υ v T v ¯ z ,
where, υ h T and υ v T are horizontal and vertical turbulent viscosities, respectively. In this study, horizontal turbulent viscosity is assumed a constant υ h T = 1.0 × 10 3 m2/s, and vertical turbulent viscosity is assumed a constant-parabolic distribution, i.e.,
υ v T ( z ) = { κ δ ( 1 δ h ) u ¯ * , c i f ( z + h ) δ κ z ( 1 + z h ) u ¯ * , c o t h e r w i s e ,
where, κ = 0.041 is the Von Karman constant, u ¯ * , c = | τ ¯ c | / ρ is friction velocity, and τ ¯ c is the bed shear stress caused by the mean current.
In the condition of stationary waves, the quasi-Eulerian mean continuity equation in the 2DV model is:
u ¯ x + w ¯ z = Z ¯ S t u ¯ S x w ¯ S z ,
If the bed level is fixed then the depth-integrated continuity equation is simplified as:
ζ ¯ t + x h ζ ¯ L u ¯ d z = x ( h ζ ¯ L u ¯ S d x ) ,

3.1.2. Depth-Dependent Wave Radiation Stress in the 2DV Model

The following formulas are derived with the assumption that all the surface waves are uniform in the y-direction. Components of wave radiation stress tensor in horizontal momentum equations are defined as:
S x x = ρ ( u ˜ 2 ¯ w ˜ 2 ¯ ) S x z = ρ u ˜ w ˜ ¯ ,
S y x = ρ v ˜ u ˜ ¯ S y z = ρ v ˜ w ˜ ¯ ,
With the assumption of slow modulation of the waves, the shear components of radiation stress tensor can be estimated by a local linear wave theory. The horizontal components of the wave velocity are given by:
u ˜ = k 1 σ a k cosh k ( z + h ) sinh k h cos ( k 1 x + k 2 y σ t ) ,
v ˜ = k 2 σ a k cosh k ( z + h ) sinh k h cos ( k 1 x + k 2 y σ t ) ,
The vertical component of wave velocity is calculated from the continuity equation for the wave motion:
w ˜ = σ a sinh k ( z + h ) sinh k h sin ( k 1 x + k 2 y σ t ) k 1 k 2 σ a x sinh k ( z + h ) sinh k h cos ( k 1 x + k 2 y σ t ) ,
Then, shear components of the wave radiation stress tensor are:
S x x = ρ g k a 2 ( k 1 2 k 2 cosh 2 k ( z + h ) sinh 2 k h sinh 2 k ( z + h ) sinh 2 k h ) ,
S y x = ρ g a 2 k 1 k 2 k cosh 2 k ( z + h ) sinh 2 k h ,
The vertical distribution of normal components of wave radiation stress in dissipative waves was analyzed by Deigaard and Fredsøe [29]. Their study is restricted to shallow water waves, where the horizontal wave velocity is assumed to be depth-independent. This results in the linear variation of S x z and S y z with depth. Usually, the horizontal wave velocity is a depth-dependent quantity (e.g., wind waves in deep water), in which case more general formulation for S x z and S y z are required. In general, the normal components of wave radiation stress can be decomposed into conservative and decay parts, such as:
u ˜ w ˜ ¯ = u ˜ w ˜ ¯ C S + u ˜ w ˜ ¯ D C ,
v ˜ w ˜ ¯ = v ˜ w ˜ ¯ C S + v ˜ w ˜ ¯ D C ,
where, the subscripts CS and DC represent the conservative and decay parts of the normal components of wave radiation stress, respectively.
(a) Conservative part of the normal component of the wave radiation stress in weak ambient current.
When the ambient current is small in comparison with the near-bed orbital velocity, the conservative part of the normal component of the wave radiation stress can be calculated from Equations (34)–(36) to obtain:
u ˜ w ˜ ¯ C S = g k 1 2 a sinh ( 2 k ( z + h ) ) 2 k 2 sinh 2 k h a x ,
v ˜ w ˜ ¯ C S = g k 1 k 2 a sinh ( 2 k ( z + h ) ) 2 k 2 sinh 2 k h a x ,
The above formulas agree with the results obtained by You [30] and Groeneweg [26] when the incident angle of the wave is zero ( θ = 0 0 ).
(b) Conservative part of the normal component of the wave radiation stress in a strong ambient current.
As pointed out by Supharatid, Tanaka [31], and Nielsen and You [32], the ambient current has a significant impact on the vertical distribution of wave radiation stress. Therefore, Equations (41) and (42) are no longer suitable in the presence of a strong ambient current. The normal component of the wave radiation stress is enhanced by a factor C W R = ( C W R , 1 , C W R , 2 ) representing the effect of the ambient current. Equations (41) and (42) become:
u ˜ w ˜ ¯ C S = C W R , 1 g k 1 2 a sinh ( 2 k ( z + h ) ) 2 k 2 sinh 2 k h a x ,
v ˜ w ˜ ¯ C S = C W R , 2 g k 1 k 2 a sinh ( 2 k ( z + h ) ) 2 k 2 sinh 2 k h a x ,
For regular waves, the empirical factors C W R , 1 and C W R , 2 are calculated based on the formula, which was proposed by Nielsen and You [32], i.e.,
C W R , 1 = 1 + 100 u * ¯ σ a ( z + h ) D ,
C W R , 2 = 1 + 100 v * ¯ σ a ( z + h ) D ,
where, ( u * ¯ , v * ¯ ) is the friction velocity caused by waves and currents, a is the wave amplitude, and D = h + ζ ¯ the mean of total water depth.
As indicated by Ockenden and Soulsby [33], for a substantial part of the time, the bottom shear stresses caused by random waves exceed those of the corresponding regular waves. In this study, Formulas (43) and (44) are modified to apply for random waves as follows:
C W R , 1 = 1 + 100 2 u * ¯ σ a ( z + h ) D ,
C W R , 2 = 1 + 100 2 v * ¯ σ a ( z + h ) D ,
where, a = 0.5 H r m s with H r m s is the root mean square wave height. The empirical coefficient approximates unity when the ambient current is small in comparison with the near-bed orbital velocity. The friction velocity components caused by waves and currents are calculated by:
u * ¯ = τ ¯ b , 1 / ρ v * ¯ = τ ¯ b , 2 / ρ ,
where, τ ¯ b = ( τ ¯ b , 1 , τ ¯ b , 2 ) is the total bed-shear stress caused by waves and currents. Simply, the instantaneous total bed shear stress can be decomposed as:
τ b = τ w + | τ ¯ c w | ,
where, τ w is the wave-induced bed shear stress and τ ¯ c w is the bed-shear stress caused by the mean current in the presence of waves. In this work, τ w is calculated based on the formula introduced by Soulsby [34].
For monochromatic waves, shear stress τ ¯ c w is calculated by:
τ ¯ c w = 1 2 ρ f c w u ¯ b u ¯ b 2 + 1 2 | u o r b | 2 ,
where, f c w is the friction factor of the mean current in the presence of waves [35], u b = ( u b , v b ) is the near-bed horizontal velocity, and | u o r b | is the near-bed orbital velocity amplitude.
For random waves, the Formula (51) is modified based on the approximate practical formula of Feddersen, Guza [36], i.e.,
τ ¯ c w = 1 2 ρ f c w u ¯ b ( 1.16 s ) 2 + u ¯ b 2 ,
(c) Decay-related parts of the normal component of the wave radiation stress.
Outside the bottom boundary layer, the decay-related part of wave radiation stress gradient is caused by the dissipative forcing, i.e.,
u ˜ w ˜ ¯ D C z = F b r , 1 ( z ) ρ F m x , 1 ( z ) ρ ,
v ˜ w ˜ ¯ D C z = F b r , 2 ( z ) ρ F m x , 2 ( z ) ρ ,
where, F b r = ( F b r , 1 , F b r , 2 ) represents the effect of breaking wave and roller wave, and F m x = ( F m x , 1 , F m x , 2 ) represents the wave-induced mixing. In this work, the vertical distribution of the wave-induced forcing terms F b r and F m x is estimated by empirical formulas proposed by Uchiyama, McWilliams [14].
At the bottom, the wave energy is dissipated due to bottom friction. According to Longuet-Higgins [37]:
u ˜ w ˜ ¯ D C ( h ) = k 1 D f ρ σ ,
v ˜ w ˜ ¯ D C ( h ) = k 2 D f ρ σ ,
Defining F t o t as the total of wave-induced mixing and current-induced turbulent forcing, i.e.,
F t o t , 1 ρ = 1 ρ ( F m x , 1 + τ ¯ 11 x + τ ¯ 13 z ) ,
F t o t , 2 ρ = 1 ρ ( F m x , 2 + τ ¯ 21 x + τ ¯ 23 z ) ,
The total of wave-induced forcing caused by the conservative part of the wave radiation stress and breaking wave and roller wave-induced forcing F w is:
F w , 1 ρ = [ ( u ˜ 2 ¯ w ˜ 2 ¯ ) x + u ˜ w ˜ ¯ C S z ] + F b r , 1 ρ ,
F w , 2 ρ = [ u ˜ v ˜ ¯ x + v ˜ w ˜ ¯ C S z ] + F b r , 2 ρ ,
Then, the sum ( F w + F t o t ) represents the total effects of wave and turbulence on the mean current.

3.1.3. Bottom Boundary Layer Thickness in the Wave–Current Interaction Condition

In the condition of waves combined with current, Van Rijn [35] proposed the following formula:
δ = 0.2 A o r b ( A o r b / k n ) 0.25 ,
where, k n = 30 z 0 is the Nikuradse roughness.
However, Equation (61) does not account for the effect of near-bed mean current on the bottom boundary layer thickness. Therefore, it is only suitable if the near-bed mean current is small in comparison with the near-bed orbital velocity. When the near-bed current is significant and is comparable to the orbital velocity the following formula is proposed:
δ = 0.2 A o r b ( A o r b / k n ) 0.25 ( 1 + | u ¯ b | / | u o r b | ) ,
It is clear that when | u ¯ b | | u o r b | the Formula (62) reduces to Formula (61).

3.2. Numerical Approximation

The 2DV equations of mean motion are discretized based on the finite difference method on a fully staggered grid (C-grid). An implicit numerical scheme has been used to discretize the equations. Finally, the tridiagonal matrix algorithm (Thomas algorithm) has been used to solve these equations. In the model, the water level is approximated at the grid point ( i , k ) , the horizontal component of velocity at ( i + 1 / 2 , k ) , and the vertical component of velocity at ( i , k + 1 / 2 ) . The advection terms are approximated following the principles described in Stelling and Busnelli [38]. This method ensures the conservation of properties near large local gradient areas. The 2DV model developed in this study is a time-domain model starting from rest and simulating to equilibrium in all cases.

3.3. Adiabatic Test

The adiabatic test, described in Bennis, Ardhuin [39], is a seemingly simple but challenging test of the derived equations since any imbalance leads to strong spurious circulations. This test was applied firstly by Ardhuin, Rascle [19]. In this, a steady monochromatic small-amplitude wave propagates over a slope without dissipation. This test has an exact solution by solving Laplace’s equation for the instantaneous velocity potential with given bottom, surface, and lateral boundary conditions [19]. In the work of Ardhuin, Rascle [19], the adiabatic test was solved by using the NTUA-nl2 model (National Technical University of Athens numerical model) developed by Belibassakis and Athanassoulis [40]. The quasi-Eulerian mean current is depth-uniform.

3.3.1. Bathymetry

The bathymetry was symmetrical and varied slowly from 4 to 6 m in the x-direction, and was uniform in the y-direction (Figure 2). The maximum bottom slope was 2.6 × 10−2, and the reflection coefficient was R = 1.4 × 10 9 , so the reflected wave in the momentum balance could be neglected [19].

3.3.2. Boundary Conditions

At the boundary, a regular wave with a height of 1.02 m and a period of 5.26 s was imposed. This is also the wave that was used by Ardhuin, Rascle [19], and Bennis, Ardhuin [39] to test their models in the adiabatic condition. The mean water level at the outflow boundary is given by:
ζ ¯ = k a 2 2 sinh 2 k D ,
At the inflow boundary, the quasi-Eulerian mean velocity is vertical uniform and given by:
u ¯ ( z ) = 1 h + ζ ¯ L h ζ ¯ L u ¯ S ( z ) d z ,
At the outflow boundary, the Neumann boundary condition is applied, i.e.,
u ¯ x = u ¯ S x ,

3.3.3. Numerical Results

Figure 3 shows the spatial distribution of Stokes drift in the x-direction. It shows that the Stokes drift was constant over the horizontal bed and the magnitude of the Stokes drift increased with a decrease in water depth and vice versa.
The comparison of the mean water level calculated by the numerical model and calculated by the formula of Longuet-Higgins and Stewart [3] is given in Figure 4. It shows a perfect agreement between the two calculation methods.
In the adiabatic condition, the total forcing F t o t , 1 was zero. The vertical distribution of wave forcing term F w , 1 / ρ is presented in Figure 5. It shows that wave-induced forcing F w , 1 / ρ was zero when the waves propagated over a flat bed. On a sloping bed, this forcing was not nil and was distributed uniformly over depth. Then, the total forcing ( F w , 1 + F t o t , 1 ) was depth-uniform in the adiabatic condition.
When the wave propagated over a slope, the change of the wave height led to the change of Stokes drift. Due to the conservation of mass and momentum, the quasi-Eulerian mean velocity also changed. However, the vertical integration of total flow was still unchanged, and in this case, it equaled zero:
h ζ ¯ L u ¯ L d z = 0 ,
Since all dissipative forcing was absent, the quasi-Eulerian mean horizontal velocity was uniformly distributed over the vertical. However, the GLM velocity inherited the non-uniformity from the Stokes drift (Figure 6a). Figure 6b presents the quasi-Eulerian mean velocity calculated with the 2DV model. It proves that quasi-Eulerian mean equations of motion passed the adiabatic test.

3.4. Mean Current in the Presence of Non-Breaking Waves

In the experiment of Klopman [41], the vertical distribution of the mean current was measured in three different types of waves: monochromatic waves, bichromatic waves, and random waves. The experiments were performed for four conditions of ambient currents: currents only (CO), waves only (WO), waves following currents (WFC), and waves opposing currents (WOC). The wave height was chosen so that no wave breaking took place. Therefore, the bottom friction plays an important role in the vertical distribution of the mean current. In the following, the experimental data for random waves were employed to validate the 2DV numerical model.

3.4.1. Input Parameters

The experiment was performed in a wave flume that has a horizontal flat bottom. The flume was 45 m long, 1 m wide, and 0.5 m deep. The total discharge was kept constant: Q = 0 m3/s for the case of wave-only, and Q = 0.08 m3/s for the remaining cases. The properties of the random waves at the wave paddle are given in Table 1.
The flow velocities were measured at the center of the channel, i.e., 22.5 m from the wave paddle. Two laser-Doppler velocimetry flow meters (LDVs) were used to measure flow velocity components. The vertical distributions of Eulerian-mean velocities measured at the center of the flume are presented in Figure 6.
In the WO condition (Figure 7a), the wave propagated from the right-hand side to the left-hand side. It shows that the wave-induced streaming near the bed was in the same direction as the propagation of the surface waves. The horizontal mean velocity changes sign at a height of approximately 0.13 m from the bed [41]. Outside the bottom streaming layer, the mean velocity varied almost linearly. In Figure 7b, the vertical distribution of horizontal mean velocity is presented in three conditions: CO, WFC, and WOC. It shows that vertical profiles of the mean current changed significantly in the presence of surface waves. In the WFC condition, the velocity shear u ¯ / z was negative in the upper part of the water column (z/h > 0.4). In the WOC condition, the mean velocity decreased near the bed (z/h < 0.4), and increased near the surface (z/h > 0.4) in comparison with the current-only condition.
By linear extrapolation of the velocities in a semilogarithmic scale, Klopman [41] obtained the friction velocity u * 7.3 × 10 3 m/s. The vertical distribution of the Reynolds shear stress u t w t ¯ is present in Figure 8. The bottom shear stress was estimated by Klopman [41] about τ ¯ b / ρ = 4.6 × 1 0 5 m2/s2, corresponding to the friction velocity of u * = τ ¯ b / ρ 6.7 × 10 3 m/s. Then, the friction velocity calculated from the bed shear stress was slightly smaller than obtained from the velocity profile.

3.4.2. Boundary Conditions

In the model, the shear stress is assumed to vanish at the mean water surface, since non-breaking waves are considered, i.e.,
[ ( υ v + K f r ) u ¯ z ] z = ζ ¯ = 0 ,
where, υ v = υ v T + υ m o l . At the bottom, the bottom boundary condition is given by:
[ ( υ v + K f r ) u ¯ z ] z = h + δ = k D f ρ σ τ ¯ c w ρ ,
At the outflow boundary, the boundary condition for the mean water level and the mean current is given by:
ζ ¯ = k a 2 2 sinh 2 k h ,
u ¯ x = u ¯ S x ,
At the inflow boundary, the quasi-Eulerian mean velocity is given by:
u ¯ ( z ) = 1 h + ζ ¯ L Q ¯ L u ¯ S ( z ) ,
where, Q ¯ L is the mean of total discharge through the pipe.

3.4.3. The Numerical Results

The experiment of Klopman [41] is simulated by the 2DV model with spatial steps of Δ x = 0.15   m , Δ z = 0.0025   m , and a time step of Δ t = 0.5   s . In the experiment, due to small technical issues, there was uncertainty in the measured discharge (see Klopman [41] for more detail). However, these errors were not corrected in his document. Generally, the measured discharge is expressed as a total of the real discharge and error discharge, i.e.,
Q m e a s u r e d = Q r e a l + Q e r r ,
where, Q r e a l is the real discharge through the wave flume, Q m e a s u r e d is the measured discharge, and Q e r r is the error of flow discharge. In CO condition, it found that when Q e r r approximates to 0.003 m3s−1 (3.75% of the real discharge) a good agreement between numerical results and experimental data was obtained. In waves combined with current conditions, the error of flow discharge is assumed similar to the current only condition. In the WO condition, flow discharge Q e r r is zero by definition. In all tests, the bed roughness is kept constant z 0 = 4.0 × 10 5 m, corresponding to a Nikuradse roughness of k n = 1.2 × 10 3 m. In Table 2, the bottom boundary thickness is presented. In the condition of waves combined with current, bottom boundary thickness δ was calculated by two methods: the formula of Van Rijn [35] and its modified Formula (62). The results are presented in Table 2.
Table 3 presents the characteristics of the mean flow near the bed calculated by the 2DV model. It shows that bottom stress in conditions of waves combined with the current was much higher than that in conditions of WO and CO. This is because momentum mixing under wave-current interaction conditions was much higher than in conditions of both WO and CO (see for instance the discussion in Chapter 3 of [42]).
The vertical distribution of Reynolds turbulent viscosity is presented in Figure 9.
It shows that the viscosity in the WFC condition was bigger than in the WOC condition, and viscosity in waves combined with current conditions was bigger than in the CO condition. Moreover, the viscosity in the WO condition was much smaller than for other conditions.
In Figure 10, the conservative part of u ˜ w ˜ ¯ in different conditions of waves combined with the current was plotted. It clearly shows that the ambient current had a significant impact on the normal component of the wave radiation stress. Moreover, the conservative part of the normal component of wave radiation stress u ˜ w ˜ ¯ in the condition of the following waves was slightly bigger than that in the condition of the opposing waves.
In non-breaking waves, the wave-induced forcing term F w , 1 / ρ only represents the effect of the conservative part of the wave radiation stress. Figure 11a–c show the vertical distributions of forcing term F w , 1 / ρ and mixing term F t o t , 1 / ρ at the center of the wave flume. In all three tests, i.e., WO, WFC, and WOC, the forcing F w , 1 is completely compensated by total mixing F t o t , 1 at any water depth level. The total of these two forcing terms, i.e., ( F w , 1 + F t o t , 1 ) , is depth-uniform in the non-breaking wave condition.
Figure 12 shows the vertical distribution of the Reynolds turbulent stress at the center of the wave flume. Near the bed, the turbulent stress calculated by the 2DV numerical model was about τ ¯ b / ρ 54.4 × 10 6 m2/s2, and the corresponding friction velocity was 7.4 × 10 3 m/s (Table 4), which was in good agreement with the friction velocity obtained from the mean velocity profile by Klopman [41], i.e., 7.3 × 10 3 m/s.
The spatial distribution of the quasi-Eulerian mean velocity calculated with the 2DV model in the condition of CO is given in Figure 13. The (vertically uniform) inflow boundary was specified at position x = 0 m.
The comparison between numerical results and experimental data in the middle of the flume is given in Figure 14. The comparison is given in both the linear scale and semilogarithmic scale. It shows that the mean current profile calculated with the 2DV model closely followed the experimental data. The agreement was good not only in the upper part of the water column but also close to the bed.
In Figure 15, the spatial distributions of Stokes drift and quasi-Eulerian mean velocity in the WO condition are presented. In this, surface waves are imposed at x = 45 m and propagated towards the left. The Stokes drift was in the direction of the waves, whereas the quasi-Eulerian mean velocity was in the opposite direction. The magnitude of Stokes drift and quasi-Eulerian mean velocity decreased in the wave propagation direction due to the effect of bottom friction. However, the momentum transport of total flow through any vertical section was zero.
The comparison between the quasi-Eulerian mean velocity calculated with the 2DV model and experimental data in the WO condition is given in Figure 16. It shows a very good agreement between model results and experimental data in the whole vertical section even on both linear scale and semilogarithmic scale. The near-bed wave-induced streaming is in the wave propagation direction. The mean horizontal velocity varied near-linearly from above the bottom streaming layer up to the mean surface level.
The spatial distribution of quasi-Eulerian mean velocity in conditions of WFC and WOC are presented in Figure 17a,b. In the case of WFC, the magnitude of the mean velocity was not the biggest at the water surface but located inside the water column. In the case of WOC, the magnitude of velocity increased from the bottom to the water surface.
Figure 18a,b shows the comparisons between experimental data and numerical model results under the condition of waves combined with the current in a linear scale. It shows that the vertical profile of the mean velocity calculated by the 2DV model fit well with experimental data in the whole water column. The agreement holds for both WFC and WOC conditions.
In Figure 19a,b, the semilogarithmic scale is employed to see the agreement between the 2DV model results and experimental results. It shows that the agreement between the 2DV model results and experimental data was very good even in the area very close to the bed. Besides, the mean current velocity profiles using the boundary layer thickness formulation proposed by Van Rijn [35] (Equation (61)) were also plotted (dashed line). The result calculated by the modified formulation (62) was better than the result calculated by the formulation of Van Rijn [35], especially in the region close to the bed. It is noticed that the formula to calculate boundary layer thickness δ , i.e., Formula (62), is an extension of the formula proposed by Van Rijn [35] (Formula (61)). In Formula (38), the effect of near-bed current is accounted for when calculating δ . Table 2 presents boundary layer thickness δ in different conditions of waves combined with the current. In the wave-only condition, Formulas (38) and (37) are identical. However, the difference between them is significant in cases of a strong ambient current. Comparison presents in Figure 18 suggest that near-bed current should be accounted for to calculate bottom boundary layer thickness, especially in the case of a strong ambient current.

3.5. Breaking Waves Propagating in a Wave Flume

3.5.1. Bathymetry and the Wave Properties at the Boundary

The experiment of Boers [43] was carried out in a wave flume with dimensions of 40 m long, 1.05 m high, and 0.8 m wide. The bottom of the flume is filled with sand and mortar on the top layer. Two breaker bars are designed at the bottom. The still water level was fixed at the level z = 0 m. The bottom profile of the calculation area is depicted in Figure 20.
Two wave conditions with the highest and lowest wave height at the boundary in the experiment of Boers [43] were employed in this work (tests 1B and 1C). The properties of the waves at the offshore boundary are given in Table 4, where H s was the significant wave height of the waves.

3.5.2. Boundary Conditions

(a) Surface and bottom boundary conditions:
At the GLM water surface, the total shear stress is assumed to vanish, i.e.,
[ ( υ v + K f r + K b ) u ¯ z ] z = ζ ¯ L = 0 ,
where, υ v = υ m o l + υ v T
At the bed, the momentum dissipated by bottom friction is assumed to be transferred to the vertical mixing, i.e.,
[ ( υ v + K f r + K b ) u ¯ z ] z = h + δ = k 1 D f ρ σ τ ¯ c w , x ρ ,
(b) Lateral boundary conditions:
The quasi-Eulerian mean water level at the offshore boundary is:
ζ ¯ = k a 2 2 sinh 2 k D ,
The quasi-Eulerian mean velocity at the offshore boundary is:
u ¯ ( z ) = 1 h + ζ ¯ L h ζ ¯ L u ¯ S ( z ) d z ,
At the land boundary, the GLM velocity in the cross-shore direction is zero then:
u ¯ ( z ) = u ¯ S ( z ) ,

3.5.3. Model Validation

In the experiment, the cross-shore distribution of significant wave height was calculated from the zeroth-order spectral moment of surface elevation [43], that is:
H s = 4 0 f N S ( f ) d f ,
where, S ( f ) is the spectral energy density, and f N is Nyquist frequency of the measurements ( 10   Hz ).
The significant wave height in the 2DV numerical model is calculated from the wave energy balance equation. In Figure 21, the spatial distributions of significant wave height in tests 1B and 1C are presented. In this, the dots present the experimental data, and the solid line presents the numerical model results. It shows that the numerical model results fit well with measured data. In test 1B, the experimental data was slightly lower than the model results in the region from x = 5 m to x = 14 m. In test 1C, the recirculation was rather small then a better agreement between two kinds of data was obtained.
Figure 22a,b show the vertical distribution of forcing term F w , 1 / ρ and total mixing term F t o t , 1 / ρ at x = 22.9   m . It shows that the wave-induced forcing F w , 1 / ρ is completely balanced by total mixing F t o t , 1 / ρ . Total forcing ( F w , 1 + F t o t , 1 ) is depth—uniform in breaking wave conditions.
The comparison of the mean water level calculated by the 2DV numerical model with measured data is given in Figure 23. In both tests, the calculated mean water level fits well with measured data even in the breaking wave area.
The horizontal distribution of Stokes drifts and quasi-Eulerian mean velocities calculated by the 2DV model is presented in Figure 24. At the locations near the sandbars, it shows that the quasi-Eulerian mean velocity was increased in the shoreward direction. This is due to the increase of Stokes drift in the onshore direction. In the breaking wave area, it shows a large vertical shear of the quasi-Eulerian mean velocity.
The comparison between quasi-Eulerian mean velocities calculated by the 2DV numerical model and measured data along the wave flume is given in Figure 25 and Figure 26.
Overall, the 2DV model simulated quite well the vertical distribution of the mean velocity. In comparison with test 1B, the model results for the test 1C show a better agreement with experimental data, especially near the sandbars. It suggests that empirical formulas of Uchiyama, McWilliams [14] can be applied well for small-amplitude waves. For the waves of high-amplitude, these empirical formulas just give qualitative results in the breaking zone and further research on the wave forcing in breaking wave areas is needed.

3.6. Breaking Waves Propagating in a Large-Scale Facility

3.6.1. Laboratory Setup and Boundary Conditions

(a) Laboratory setup:
In this section, the 2DV model was employed to simulate the cross-shore and the longshore currents with the input data obtained from the large-scale sediment transport facility (LSTF). The design of the LSTF is presented in detail in the paper of Hamilton and Ebersole [44]. The facility has dimensions of approximately 30 m cross-shore by 50 m longshore by 1.4 m deep. The concrete beach has a longshore dimension of 31 m and a cross-shore dimension of 21 m, with a slope of 1:30. Unidirectional long-crested waves were generated with four piston-type wave generators. An active pumping and recirculation system comprised of 20 independent pumps and pipelines were used to control the cross-shore distribution of mean longshore current. Pumping rates were adjusted iteratively to converge toward the proper setting, based on the measurements along the beach.
The wave conditions at the offshore boundary are given in Table 5. The TMA spectrum (self-similar spectral shape) with the width parameter of 3.3 was employed. It is an extension of deep water spectrum JONSWAP for applications in shallow water.
(b) Surface and bottom boundary conditions:
At the mean water surface, the total shear stress is assumed to vanish, i.e.,
[ ( υ v + K f r + K b ) u ¯ z ] z = ζ ¯ L = 0 ,
[ ( υ v + K f r + K b ) v ¯ z ] z = ζ ¯ L = 0 ,
where, υ v = υ m o l + υ v T .
The bottom boundary condition in the cross-shore direction is given by:
[ ( υ v + K f r + K b ) u ¯ z ] z = h + δ = k 1 D f ρ σ τ ¯ c w , x ρ ,
According to Longuet-Higgins [4], the bed shear stress in the longshore direction can be calculated based on the local rate of energy dissipation. Therefore:
τ ¯ c w , y = k 2 [ ( 1 α r ) D w + D r + D f ] σ ,
From Equations (52) and (82) the bottom boundary condition in the longshore direction is given by:
v ¯ b = 2 k 2 [ ( 1 α r ) D w + D r + D f ] ρ σ f c w ( 1.16 s ) 2 + | u ¯ b | 2 ,
(c) Lateral boundary conditions:
At the offshore, boundary conditions for quasi-Eulerian mean water level and quasi-Eulerian mean velocity components are:
ζ ¯ = k a 2 2 sinh 2 k D ,
u ¯ = 1 h + ζ ¯ L h + δ ζ ¯ L u ¯ S d z ,
v ¯ = u ¯ tan θ ,
At the land boundary, the GLM velocity in the x-direction is zero. Besides, the no-slip boundary condition is assumed for quasi-Eulerian velocity in the y-direction. Then:
u ¯ ( z ) = u ¯ S ( z ) ,
v ¯ ( z ) = 0 ,

3.6.2. Numerical Results and Discussion

Figure 27 shows the comparison of significant wave height between experimental data and numerical results calculated by the 2DV model. The comparison shows a very good agreement between experimental data and model results.
In Figure 28, a comparison of the mean water level between experimental data and the 2DV model result is given. There was a small difference of about a few millimeters outside the breaking zone. The difference between these two kinds of data was bigger in the area closer to the coastline. It might be due to the recirculation system of the facility or because the experimental data was for a closed basin, so the volume of water in the setup was taken from offshore.
(a) Cross-shore direction:
Figure 29 presents the vertical distribution of wave forcing term F w , 1 / ρ and total mixing F t o t , 1 / ρ at the location x = 10.9   m . It shows that momentum caused by wave forcing was completely compensated by total mixing. The total of these two forcing terms was depth-uniform.
Figure 30 shows an overview of the spatial distribution of cross-shore velocity calculated by the 2DV model. It shows that outside the wave breaking zone, the vertical distribution of the quasi-Eulerian mean velocity was almost uniform. However, inside the wave breaking zone, the mean velocity shows a strong vertical shear.
The LSTF data was also employed by Teles, Pires-Silva [45] to validate the TELEMAC 3D model. In Figure 31, comparisons of cross-shore mean velocities between experimental data, results of the TELEMAC 3D model, and results of the 2DV model were carried out at six cross-shore locations.
The results obtained from the 2DV model fit quite well with the experimental data, and much better than the results obtained by the TELEMAC 3D model in most of the cross-sections. Closer to the surface and the land boundary the difference between experimental data and 2DV model results became bigger. The difference might be due to the use of empirical formulas for the wave-induced forcing, and the return flow in the facility. Improvement of these formulas making use of the large body of literature on the return flow profiles is recommendable but outside the scope of this study.
(b) Longshore direction:
Figure 32 shows the vertical distribution of the wave-induced forcing F w , 2 / ρ and total mixing F t o t , 2 / ρ in the longshore direction at x = 10.9   m . It shows that the total mixing F t o t , 2 is balanced by wave-induced forcing F t o t , 2 . The total forcing ( F w , 2 + F t o t , 2 ) is depth-uniform.
The cross-shore distribution of the longshore velocity calculated by the 2DV numerical model is presented in Figure 33. According to the results, the longshore velocity increases shoreward. The maximum magnitude of the longshore velocity reached approximately 0.35 m/s at the position x 14 m. After this point, the longshore velocity was decreased toward the shoreline.
In Figure 34, a comparison between experimental data and numerical results at one-third of the water depth above the bottom is given. It shows that the results obtained by the 2DV model agreed well with the experimental data and matched the observed cross-shore distribution better than the results obtained by the TELEMAC 3D model.
The comparison of the longshore velocities at different vertical cross-sections is given in Figure 35. The results of the TELEMAC 3D model show quite good agreement at four locations from x = 9.5 m to x = 13.9 m. However, at the location near the offshore x = 7.1 m and the location near the coastline x = 15.3 m the differences between experimental data and results from the TELEMAC 3D model were significant. In contrast, the results of the 2DV model show good agreement with experimental data at all six cross-sections.

4. Discussion

Following Van Rijn [35], the bottom boundary layer thickness δ was calculated by the Formula (61). In this, the effect of mean current on bottom boundary layer thickness was neglected. However, interactions of current with wave boundary layer and waves with current boundary layer led to the enhancement of current-induced friction, wave energy dissipation, and bed shear stress. Then, both waves and current should be considered in calculating bottom boundary thickness. In this paper, the formula given by Van Rijn [35] was modified to take into account the effect of ambient current. The use of the modified formula obtained a better agreement with experimental data than the use of the original formula of Van Rijn [35].
In the work of Nielsen and You [32], an empirical coefficient C W R was proposed to account for the effect of strong ambient current on the vertical distribution of wave radiation stress components. However, that coefficient is only applied to regular waves. In this work, a modified formula was proposed to apply for random waves. With that modification, the mean current profiles calculated by the 2DV model agreed well with experimental data obtained by Klopman [41].
In breaking wave conditions, dissipative wave forcing terms were estimated by using empirical formulas introduced by Uchiyama, McWilliams [14]. The cross-shore velocity profiles obtained by the 2DV model showed good agreement with experimental data presented in the works of Boers [43] and Hamilton and Ebersole [44]. However, there are differences between model results and observed data near the breaking points. It requires further research on the vertical distribution of wave breaking induced forcing. In the longshore direction, the effect of breaking wave-induced forcing was small. Then, the longshore current data calculated by the 2DV model fitted very well with experimental data.
The LSTF test was also employed by Teles, Pires-Silva [45] to validate the TELEMAC 3D model. That model was developed based on the work of Ardhuin, Rascle [19], and Bennis, Ardhuin [39]. As indicated by Ardhuin, Rascle [19], their equations are only expected to give qualitative results for surf zone applications. The comparison presented in the previous part also showed that the 2DV model obtained a better agreement with observed data than TELEMAC 3D model.

5. Conclusions

In this paper, a new set of equations of motion written in terms of quasi-Eulerian mean velocity was developed based on the GLM method. The new equations were valid from the bottom to the mean water surface even in the presence of finite-amplitude waves. These equations are practical for a wide range of applications from deep water to shallow water areas. All terms in the equations are expressed to the second-order of the wave amplitude. Both non-wave forcing and wave-induced forcing terms are under consideration. When the wave height is infinitesimal, the quasi-Eulerian mean equations of motion reduce to the classical Eulerian mean equations of motion. In cases of density stratification, the buoyancy effect should be included as external forcing.
In the work of Longuet-Higgins and Stewart [2], the effects of waves on the mean current are expressed in terms of the traditional radiation stress concept. It is a depth-integrated term and is only suitable for the depth-averaged equations of motion. In 3D problems, that concept is no longer suitable. In this paper, a depth-dependent wave radiation stress concept was introduced. Then, the effects of surface waves on the mean current could be specified at any level of the water column. The vertical distribution of the mean current was described in detail with the use of a depth-dependent wave radiation stress tensor. With the use of an empirical coefficient C W R , the depth-dependent radiation stress could be calculated in the condition of a strong ambient current. It is noticed that the vertical integration of the new depth-dependent wave radiation stress coincided with the traditional radiation stress.
A 2DV numerical model was developed to validate the new quasi-Eulerian mean equations of motion. The turbulent viscosity and wave-induced mixing were calculated by a simple submodel with the assumption of constant-parabolic distribution. The model shows a better comparison with experimental data than previous studies. The model passed the well-known adiabatic condition suggested by Ardhuin, Rascle [19]. It does not produce spurious velocities when the waves propagate on a sloping bottom. As a result, vertical uniform distribution of quasi-Eulerian mean horizontal velocity was obtained even on the slope.
Subsequently, the experiment of Klopman [41] for random waves was employed to validate the 2DV model in the condition of non-breaking waves combined with a current. The comparison between experimental data and numerical results showed very good agreement in the whole vertical section even in areas that were very close to the bed. In the condition of breaking waves, two experiments presented in the dissertation of Boers [43] were used. It was shown that the new equations performed very well for the experiment of smaller wave height (test 1C). In the experiment of bigger wave height (test 1B), the vertical distribution of the mean velocity near the breaking point gave qualitatively correct results. In the comparison of longshore current in the LSTF test [44], a good agreement was found not only for depth-averaged longshore current but also for depth varying longshore current. For cross-shore currents, there was a difference between model results and experimental data in the wave breaking zone. The difference between the experimental data and model results in tests 1B and LSTF was likely to be due to the empirical formulas for wave breaking that were strong simplifications of the complex breaking process. Further tuning of such formulations against a large number of datasets on wave decay and generated longshore and cross-shore currents was recommended.
Finally, with the use of quasi-Eulerian mean variables, the new set of equations of motion can be easily implemented to existing 3D models developed based on the classical Eulerian mean method. The implementation is straightforward and does not require much effort. It can improve significantly the results of simulating coastal processes such as coastal sediment transport, transport of plastic, and other pollutants such as oil slicks.

Author Contributions

D.T.N.: Ph.D. fellow; D.R.: supervisor; N.G.J.: co-supervisor. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded in a part by Deltares Strategic Research Fund.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

Experimental data presented in this study was obtained from: (1) Klopman [41]: non-breaking waves propagate in a wave flume with a flat bed. (2) Boers [43]: breaking waves propagate in a wave flume with a sloping bed. (3) Hamilton and Ebersole [44]: breaking waves propagate in the LSTF facility.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of the GLM Momentum Equation

In the following, the momentum equation of the mean motion written in terms of the GLM velocity was derived. It is the result of applying the GLM operator to the momentum equation of motion of the total flow. The equation obtained in this part is equivalent to the GLM momentum equation expressed in AM. The GLM momentum equation was used to develop a quasi-Eulerian mean momentum equation.
Averaging Equation (11) over a wave period to obtain:
( u i t + u j u i x j ) ξ ¯ + ( 2 ( Ω × u ) i ) ξ ¯ + ( Φ x i ) ξ ¯ + ( 1 ρ p x i ) ξ ¯ + ( X i ) ξ ¯ = 0 ,
According to AM, the first term on the left-hand side of Equation (A1) can be rewritten as:
( u i t + u j u i x j ) ξ ¯ = u ¯ i L t + u ¯ j L u ¯ i L x j ,
Since Ω is constant, the second term on the left-hand side of Equation (A2) becomes:
( 2 ( Ω × u ) i ) ξ ¯ = 2 ( Ω × u ¯ L ) i ,
Since only gravitational potential is considered in the term Φ ( x , t ) , the third term in the left-hand side of the Equation (A1) can be expressed as:
( Φ x i ) ξ ¯ = ( δ i 3 g ) ξ ¯ = δ i 3 g ,
When . u = 0 then according to AM:
. ξ = O ( ε 2 ) ,
Using Taylor expansion the pressure gradient term in the Equation (A1) can be expressed as:
1 ρ ( p x i ) ξ ¯ = 1 ρ p ¯ x i + 1 ρ x j ( ξ j p x i ) ¯ + 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k + O ( ε 3 ) ,
Inserting Equations (A2)–(A4), (A6) into Equation (A1) to obtain the momentum equation in the GLM framework:
D ¯ L u ¯ i L + 2 ( Ω × u ¯ L ) i + δ i 3 g + 1 ρ p ¯ x i + X ¯ i L = 1 ρ x j ( ξ j p x i ) ¯ 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k + O ( ε 3 ) ,

Appendix B. Derivation of the Radiation Stress Tensor in the GLM Framework

In this appendix, the relationship between the three-dimensional radiation stress gradient and other terms in the GLM momentum Equation (A7) is presented. From this, the GLM momentum Equation (A7) can be written in terms of quasi-Eulerian mean velocity.
For any quantity φ , AM obtained the following relationship between Lagrangian disturbance and quasi-Eulerian disturbance:
φ l = φ + ξ j φ ¯ x j + O ( ε 2 ) ,
Using (A8) the first term on the right-hand side of Equation (14) can be expressed as:
1 ρ x j ( ξ j p x i ¯ ) = 1 ρ x j [ ξ j ( p x i ) l ¯ ] 1 ρ x j ( ξ j ξ k ¯ 2 p ¯ x k x i ) + O ( ε 3 ) ,
Subtracting 14 from 11 to obtain the following equation for the evolution of disturbance motion:
D ¯ L u i l + ( 2 ( Ω × u ) i ) l + ( δ i 3 g ) l + ( X i ) l + ( 1 ρ p x i ) l = 0 ,
Since the gravitational acceleration g is assumed constant, the third term on the left-hand side of the Equation (A10) is neglected. Multiply Equation (A10) with ξ j and then take the spatial derivative / x j to obtain:
1 ρ x j [ ξ j ( p x i ) l ¯ ] = x j ( ξ j D ¯ L u i l ¯ ) x j ( 2 ξ j ( Ω × u l ) i ¯ ) x j ( ξ j X i l ¯ ) ,
The first term in the right-hand side of Equation (A11) is decomposed as:
x j ( ξ j D ¯ L u i l ¯ ) = x j ( D ¯ L ξ j u i l ¯ ) x j ( u i l D ¯ L ξ j ¯ ) ,
The first term on the right-hand side of A12 can be expressed as:
x j ( D ¯ L ξ j u i l ¯ ) = D ¯ L ( ξ j u i l x j ¯ ) + ( ξ j u i l ¯ ) x k u ¯ k L x j + O ( ε 3 ) ,
From the definition of Lagrangian disturbance, the second term on the right-hand side of A12 becomes:
x j ( u i l D ¯ L ξ j ¯ ) = ( u i l u j l ¯ ) x j ,
Substitute Equations (A13) and (14) into the Equation (A12):
x j ( ξ j D ¯ L u i l ¯ ) = D ¯ L ( ξ j u i l x j ¯ ) + ( ξ j u i l ¯ ) x k u ¯ k L x j ( u i l u j l ¯ ) x j + O ( ε 3 ) ,
In the following, all terms on the right-hand side of Equation (A15) are rewritten in terms of quasi-Eulerian disturbance based on the relationship (A8). Then, the first term in the right-hand side of Equation (A15) can be rewritten as:
D ¯ L ( ξ j u i l x j ¯ ) = D ¯ L ( ξ j u i x j ¯ ) + D ¯ L [ x j ( ξ j ξ k ¯ u ¯ i x k ) ] = D ¯ L ( ξ j u i x j ¯ ) + x j [ D ¯ L ( ξ j ξ k ¯ u ¯ i x k ) ] ξ j ξ k ¯ x l u ¯ i x k u ¯ l x j ξ j ξ k ¯ u ¯ l x j 2 u ¯ i x k x l + O ( ε 3 ) ,
Similarly, the second term on the right-hand side of (A15) becomes:
( ξ j u i l ¯ ) x k u ¯ k L x j = ( ξ j u i ¯ ) x k u ¯ k L x j + ( ξ j ξ l ¯ ) x k u ¯ i x l u ¯ k x j + ξ j ξ l ¯ u ¯ k x j 2 u ¯ i x k x l = x k ( ξ j u i ¯ u ¯ k x j ) + ( ξ j ξ l ¯ ) x k u ¯ i x l u ¯ k x j + ξ j ξ l ¯ u ¯ k x j 2 u ¯ i x k x l + O ( ε 3 ) ,
The third term on the right-hand side of Equation (A15) is expressed by:
( u i l u j l ¯ ) x j = ( u i u j ¯ ) x j + x j [ u i ξ l ¯ u ¯ j x l + u j ξ k ¯ u ¯ i x k + ξ k ξ l ¯ u ¯ i x k u ¯ j x l ] + O ( ε 3 ) ,
Substituting Equations (A16)–(A18) into Equation (A15):
x j ( ξ j D ¯ L u i l ¯ ) = ( u i u j ¯ ) x j x j [ u j ξ k ¯ u ¯ i x k + ξ k ξ l ¯ u ¯ i x k u ¯ j x l ] + D ¯ L ( ξ j u i x j ¯ ) + x j [ D ¯ L ( ξ j ξ k ¯ u ¯ i x k ) ] + O ( ε 3 ) ,
Similarly, the second and the third terms in the right-hand side of (A11) can be expressed in term of quasi-Eulerian disturbance quantities, such as:
x j ( 2 ξ j ( Ω × u l ) i ¯ ) = 2 ε i m n Ω m ξ j u n l x j ¯ + O ( ε 3 ) = 2 ε i m n Ω m ξ j u n x j ¯ + x j ( 2 ε i m n Ω m ξ j ξ k ¯ u ¯ n x k ) + O ( ε 3 ) ,
where, ε i m n is the Levi-Civita symbol defined by:
ε i m n = { + 1 if ( i , m , n )   is   an   even   permutation   of   ( 1 , 2 , 3 ) 1 if   ( i , m , n )   is   an   odd   permutation   of   ( 1 , 2 , 3 )     0 if   any   index   is   repeated ,
The third terms on the right-hand side of Equation (A11) is rewritten as:
x j ( ξ j X i l ¯ ) = ξ j X i l x j ¯ = ξ j X i x j ¯ + x j ( ξ j ξ k ¯ X ¯ i x k ) + O ( ε 3 ) ,
From Equations (A19), (A20), and (A22) the left-hand side of Equation (A11) can be expressed in term of quasi-Eulerian disturbance quantities as:
1 ρ x j [ ξ j ( p x i ) l ¯ ] = ( u i u j ¯ ) x j + x j [ u j ξ k ¯ u ¯ i x k + ξ k ξ l ¯ u ¯ i x k u ¯ j x l ] D ¯ L ( ξ j u i x j ¯ ) x j [ D ¯ L ( ξ j ξ k ¯ u ¯ i x k ) ] 2 ε i m n Ω m ξ j u n x j ¯ x j ( 2 ε i m n Ω m ξ j ξ k ¯ u ¯ n x k ) ξ j X i x j ¯ x j ( ξ j ξ k ¯ X ¯ i x k ) + O ( ε 3 ) ,
Replacing Equation (A23) into Equation (A9):
1 ρ x j ( ξ j p x i ¯ ) = ( u i u j ¯ ) x j D ¯ L ( ξ j u i x j ¯ ) 2 ε i m n Ω m ξ j u n x j ¯ ξ j X i x j ¯ + x j [ u j ξ k ¯ u ¯ i x k + ξ k ξ l ¯ u ¯ i x k u ¯ j x l ] x j [ u ¯ i x k D ¯ L ( ξ j ξ k ¯ ) ] x j [ ξ j ξ k ¯ D ¯ L ( u ¯ i x k ) ] x j ( 2 ε i m n Ω m ξ j ξ k ¯ u ¯ n x k ) x j ( ξ j ξ k ¯ X ¯ i x k ) 1 ρ x j ( ξ j ξ k ¯ 2 p ¯ x k x i ) + O ( ε 3 ) ,
The sixth term on the right-hand side of Equation (A24) can be expressed as:
u ¯ i x k D ¯ L ( ξ j ξ k ¯ ) = u ¯ i x k [ ξ j D ¯ L ξ k ¯ + ξ k D ¯ L ξ j ¯ ] = u ¯ i x k ( ξ j u k l ¯ + ξ k u j l ¯ ) = ξ j u k ¯ u ¯ i x k + ξ k u j ¯ u ¯ i x k + ξ j ξ l ¯ u ¯ k x l u ¯ i x k + ξ k ξ l ¯ u ¯ j x l u ¯ i x k + O ( ε 3 ) ,
Similarly, the seventh term on the right-hand side of Equation (A24) is expressed as:
ξ j ξ k ¯ D ¯ L ( u ¯ i x k ) = ξ j ξ k ¯ x k ( D ¯ L u ¯ i ) ξ j ξ k ¯ u ¯ l L x k u ¯ i x l + O ( ε 3 ) ,
From Equation (A26) the total of the last four terms in the right-hand side of Equation (A23) becomes:
x j [ ξ j ξ k ¯ D ¯ L ( u ¯ i x k ) + ( 2 ε i m n Ω m ξ j ξ k ¯ u ¯ n x k ) + ( ξ j ξ k ¯ X ¯ i x k ) + 1 ρ ( ξ j ξ k ¯ 2 p ¯ x k x i ) ] = x j [ ξ j ξ k ¯ x k ( D ¯ L u ¯ i + 2 ( Ω × u ¯ ) i + δ i 3 g + X ¯ i + 1 ρ p ¯ x i ) ] x j ( ξ j ξ k ¯ u ¯ l x k u ¯ i x l ) + O ( ε 3 ) ,
Since φ ¯ S = O ( ε 2 ) then from Equation (14), we have:
D ¯ L u ¯ i + 2 ( Ω × u ¯ ) i + δ i 3 g + X ¯ i + 1 ρ p ¯ x i = O ( ε 2 ) ,
Therefore, in the second-order of the accuracy of small disturbance amplitude Equation (A27) can be rewritten as:
x j [ ξ j ξ k ¯ D ¯ L ( u ¯ i x k ) + 2 ε i m n Ω m ( ξ j ξ k ¯ u ¯ n x k ) + ( ξ j ξ k ¯ X ¯ i x k ) + 1 ρ ( ξ j ξ k ¯ 2 p ¯ x k x i ) ] = x j ( ξ j ξ k ¯ u ¯ l x k u ¯ i x l ) + O ( ε 3 ) ,
From Equations (A25) and (A29) the sum of the last six terms in Equation (A24) can be expressed by:
x j [ u j ξ k ¯ u ¯ i x k + ξ k ξ l ¯ u ¯ i x k u ¯ j x l ] x j [ u ¯ i x k D ¯ L ( ξ j ξ k ¯ ) ] x j [ ξ j ξ k ¯ D ¯ L ( u ¯ i x k ) + 2 ε i m n Ω m ( ξ j ξ k ¯ u ¯ n x k ) + ( ξ j ξ k ¯ X ¯ i x k ) + 1 ρ ( ξ j ξ k ¯ 2 p ¯ x k x i ) ] = x j ( ξ j u k ¯ u ¯ i x k ) + O ( ε 3 ) ,
Substituting Equation (A30) into Equation (A24) to obtain the following relationship:
1 ρ x j ( ξ j p x i ¯ ) = ( u i u j ¯ ) x j D ¯ L ( ξ j u i x j ¯ ) 2 ε i m n Ω m ξ j u n x j ¯ ξ j X i x j ¯ x j [ ξ j u k ¯ u ¯ i x k ] + O ( ε 3 ) ,
Following the definition of Stokes correction (AM) the second term on the right-hand side of Equation (A31) can be expressed as:
D ¯ L ( ξ j u i x j ¯ ) = D ¯ L u ¯ i S D ¯ L ( 1 2 ξ j ξ k ¯ 2 u ¯ i x j x k ) + O ( ε 3 ) ,
The second term on the right-hand side of Equation (A32) is expressed as:
D ¯ L ( 1 2 ξ j ξ k ¯ 2 u ¯ i x j x k ) = 1 2 ξ j ξ k ¯ D ¯ L ( 2 u ¯ i x j x k ) + 1 2 2 u ¯ i x j x k D ¯ L ( ξ j ξ k ¯ ) = 1 2 ξ j ξ k ¯ [ 2 x j x k ( D ¯ L u ¯ i ) 2 u ¯ l x j x k u ¯ i x l u ¯ l x k 2 u ¯ i x j x l u ¯ l x j 2 u ¯ i x k x l ] + 1 2 ( ξ j u k ¯ + ξ j ξ l ¯ u ¯ k x l + ξ k u j ¯ + ξ k ξ l ¯ u ¯ j x l ) 2 u ¯ i x j x k + O ( ε 3 ) = 1 2 ξ j ξ k ¯ [ 2 x j x k ( D ¯ L u ¯ i ) ] 1 2 ξ j ξ k ¯ 2 u ¯ l x j x k u ¯ i x l + ξ j u k ¯ 2 u ¯ i x j x k + O ( ε 3 ) ,
The last three terms on the right-hand side of Equation (A31) become:
2 ε i m n Ω m ξ j u n x j ¯ = 2 ( Ω × u ¯ S ) i 2 ε i m n Ω m ( 1 2 ξ j ξ k ¯ 2 u ¯ n x j x k ) + O ( ε 3 ) ,
ξ j X i x j ¯ = X ¯ i S 1 2 ξ j ξ k ¯ 2 X ¯ i x j x k + O ( ε 3 ) ,
x j [ ξ j u k ¯ u ¯ i x k ] = ξ j u k x j ¯ u ¯ i x k + ξ j u k ¯ 2 u ¯ i x j x k = u ¯ k S u ¯ i x k 1 2 ξ j ξ l ¯ 2 u ¯ k x j x l u ¯ i x k + ξ j u k ¯ 2 u ¯ i x j x k + O ( ε 3 ) ,
Replacing Equations from A32 to A36 into Equation (A31):
1 ρ x j ( ξ j p x i ¯ ) + 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k = ( u i u j ¯ ) x j D ¯ L u ¯ i S 2 ( Ω × u ¯ S ) i X ¯ i S u ¯ k S u ¯ i x k + O ( ε 3 ) ,
Or equivalently:
D ¯ L u ¯ i S + 2 ( Ω × u ¯ S ) i + X ¯ i S + 1 ρ x j ( ξ j p x i ¯ ) + 1 2 ρ ξ j ξ k ¯ 3 p ¯ x i x j x k = ( u i u j ¯ ) x j u ¯ k S u ¯ i x k + O ( ε 3 ) ,
Equation (A38) shows the relationship between the evolution of Stokes drift and the disturbances of a fluid particle. This is the basis to develop quasi-Eulerian mean equations of motion.

References

  1. Longuet-Higgins, M.S.; Stewart, R.W. Changes in the form of short gravity waves on long waves and tidal currents. J. Fluid Mech. 1960, 8, 565–583. [Google Scholar] [CrossRef]
  2. Longuet-Higgins, M.S.; Stewart, R. Radiation stress and mass transport in gravity waves, with application to ‘surf beats’. J. Fluid Mech. 1962, 13, 481–504. [Google Scholar] [CrossRef]
  3. Longuet-Higgins, M.S.; Stewart, R. Radiation Stresses in Water Waves; A Physical Discussion, with Applications. In Deep Sea Research and Oceanographic Abstracts; Elsevier: Amsterdam, The Netherlands, 1964; pp. 529–562. Available online: https://www.sciencedirect.com/science/article/abs/pii/0011747164900014 (accessed on 10 January 2021).
  4. Longuet-Higgins, M.S. Longshore currents generated by obliquely incident sea waves: 1. J. Geophys. Res. 1970, 75, 6778–6789. [Google Scholar] [CrossRef]
  5. Bowen, A.J. The generation of longshore currents on a plane beach. J. Mar. Res. 1969, 27, 206–215. [Google Scholar]
  6. Thornton, E.B. Variation of longshore current across the suft zone. Coast. Eng. Proc. 1970, 1, 18. [Google Scholar] [CrossRef] [Green Version]
  7. Xie, L.; Wu, K.; Pietrafesa, L.; Zhang, C. A numerical study of wave-current interaction through surface and bottom stresses: Wind-driven circulation in the South Atlantic Bight under uniform winds. J. Geophys. Res. Oceans (1978–2012) 2001, 106, 16841–16855. [Google Scholar] [CrossRef] [Green Version]
  8. Xia, H.; Xia, Z.; Zhu, L. Vertical variation in radiation stress and wave-induced current. Coast. Eng. 2004, 51, 309–321. [Google Scholar] [CrossRef]
  9. Craik, A.; Leibovich, S. A rational model for Langmuir circulations. J. Fluid Mech. 1976, 73, 401–426. [Google Scholar] [CrossRef] [Green Version]
  10. Leibovich, S. On the evolution of the system of wind drift currents and Langmuir circulations in the ocean. Part 1. Theory and averaged current. J. Fluid Mech. 1977, 79, 715–743. [Google Scholar] [CrossRef]
  11. McWilliams, J.C.; Restrepo, J.M. The wave-driven ocean circulation. J. Phys. Oceanogr. 1999, 29, 2523–2540. [Google Scholar] [CrossRef]
  12. McWilliams, J.C.; Restrepo, J.M.; Lane, E.M. An asymptotic theory for the interaction of waves and currents in coastal waters. J. Fluid Mech. 2004, 511, 135–178. [Google Scholar] [CrossRef]
  13. Newberger, P.; Allen, J.S. Forcing a three-dimensional, hydrostatic, primitive-equation model for application in the surf zone: 1. Formulation. J. Geophys. Res. Ocean. 2007, 112. [Google Scholar] [CrossRef] [Green Version]
  14. Uchiyama, Y.; McWilliams, J.C.; Shchepetkin, A.F. Wave–current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone. Ocean Model. 2010, 34, 16–35. [Google Scholar] [CrossRef]
  15. Kumar, N.; Voulgaris, G.; Warner, J.C.; Olabarrieta, M. Implementation of the vortex force formalism in the coupled ocean-atmosphere-wave-sediment transport (COAWST) modeling system for inner shelf and surf zone applications. Ocean Model. 2012, 47, 65–95. [Google Scholar] [CrossRef]
  16. Lane, E.M.; Restrepo, J.M.; McWilliams, J.C. Wave–current interaction: A comparison of radiation-stress and vortex-force representations. J. Phys. Oceanogr. 2007, 37, 1122–1141. [Google Scholar] [CrossRef] [Green Version]
  17. Mellor, G. The three-dimensional current and surface wave equations. J. Phys. Oceanogr. 2003, 33, 1978–1989. [Google Scholar] [CrossRef]
  18. Mellor, G.L. The depth-dependent current and wave interaction equations: A revision. J. Phys. Oceanogr. 2008, 38, 2587–2596. [Google Scholar] [CrossRef]
  19. Ardhuin, F.; Rascle, N.; Belibassakis, K.A. Explicit wave-averaged primitive equations using a generalized Lagrangian mean. Ocean Model. 2008, 20, 35–60. [Google Scholar] [CrossRef] [Green Version]
  20. Mellor, G. A combined derivation of the integrated and vertically resolved, coupled wave–current equations. J. Phys. Oceanogr. 2015, 45, 1453–1463. [Google Scholar] [CrossRef]
  21. Mellor, G. On theories dealing with the interaction of surface waves and ocean circulation. J. Geophys. Res. Ocean. 2016, 121, 4474–4486. [Google Scholar] [CrossRef]
  22. Ardhuin, F.; Suzuki, N.; McWilliams, J.C.; Aiki, H. Comments on “A Combined Derivation of the Integrated and Vertically Resolved, Coupled Wave–Current Equations”. J. Phys. Oceanogr. 2017, 47, 2377–2385. [Google Scholar] [CrossRef]
  23. Andrews, D.; McIntyre, M. An exact theory of nonlinear waves on a Lagrangian-mean flow. J. Fluid Mech. 1978, 89, 609–646. [Google Scholar] [CrossRef] [Green Version]
  24. Leibovich, S. On wave-current interaction theories of Langmuir circulations. J. Fluid Mech. 1980, 99, 715–724. [Google Scholar] [CrossRef] [Green Version]
  25. Dingemans, M.W. Water Wave Propagation over Uneven Bottoms: Linear Wave Propagation; World Scientific: Singapore, 1997; Volume 13. [Google Scholar]
  26. Groeneweg, J. Wave-Current Interactions in a Generalized Lagrangian Mean Formulation. Ph.D. Thesis, Delft University of Technology, Delft, The Netherland, 1999. [Google Scholar]
  27. Walstra, D.; Roelvink, J.; Groeneweg, J. Calculation of wave-driven currents in a 3D mean flow model. In Coastal Engineering 2000; 2001; pp. 1050–1063. Available online: https://ascelibrary.org/doi/abs/10.1061/40549(276)81 (accessed on 13 January 2021).
  28. Roelvink, J.; Reniers, A. A Guide to Modeling Coastal Morphology; World Scientific: Singapore, 2011; Volume 12. [Google Scholar]
  29. Deigaard, R.; Fredsøe, J. Shear stress distribution in dissipative water waves. Coast. Eng. 1989, 13, 357–378. [Google Scholar] [CrossRef]
  30. You, Z.-J. On the vertical distribution of< ũ w ~ > by FJ Rivero and AS Arcilla: Comments. Coast. Eng. 1997, 30, 305–310. [Google Scholar] [CrossRef]
  31. Supharatid, S.; Tanaka, H.; Shuto, N. Interactions of waves and current (Part I: Experimental investigation). Coast. Eng. Jpn. 1992, 35, 167–186. [Google Scholar] [CrossRef]
  32. Nielsen, P.; You, Z.-J. Eulerian mean velocities under non-breaking waves on horizontal bottoms. In Coastal Engineering 1996; 1997; pp. 4066–4078. Available online: https://ascelibrary.org/doi/abs/10.1061/9780784402429.314 (accessed on 10 January 2021).
  33. Ockenden, M.; Soulsby, R. Sediment Transport by Currents Plus Irregular Waves. 1994. Available online: https://eprints.hrwallingford.com/347/ (accessed on 13 January 2021).
  34. Soulsby, R. Dynamics of Marine Sands: A Manual for Practical Applications. Thomas Telford, 1997. Available online: https://www.infona.pl/resource/bwmeta1.element.elsevier-69102d77-8480-3ca9-adc0-a09c289eb654 (accessed on 13 January 2021).
  35. Van Rijn, L.C. Principles of Fluid Flow and Surface Waves in Rivers, Estuaries, Seas and Oceans; Aqua Publications: Amsterdam, The Netherlands, 2011; Available online: https://www.aquapublications.nl/Contentsbook1.pdf (accessed on 13 January 2021).
  36. Feddersen, F.; Guza, R.; Elgar, S.; Herbers, T. Velocity moments in alongshore bottom stress parameterizations. J. Geophys. Res. Ocean. 2000, 105, 8673–8686. [Google Scholar] [CrossRef]
  37. Longuet-Higgins, M. The mechanics of the surf zone. In Theoretical and Applied Mechanics; Springer: Berlin/Heidelberg, Germany, 1973; pp. 213–228. [Google Scholar] [CrossRef]
  38. Stelling, G.S.; Busnelli, M.M. Numerical simulation of the vertical structure of discontinuous flows. Int. J. Numer. Methods Fluids 2001, 37, 23–43. [Google Scholar] [CrossRef]
  39. Bennis, A.-C.; Ardhuin, F.; Dumas, F. On the coupling of wave and three-dimensional circulation models: Choice of theoretical framework, practical implementation and adiabatic tests. Ocean Model. 2011, 40, 260–272. [Google Scholar] [CrossRef] [Green Version]
  40. Belibassakis, K.; Athanassoulis, G. Extension of second-order Stokes theory to variable bathymetry. J. Fluid Mech. 2002, 464, 35–80. [Google Scholar] [CrossRef] [Green Version]
  41. Klopman, G. Vertical Structure of the Flow Due to Waves and Currents: Laser-Doppler Flow Measurements for Waves following or Opposing a Current. Delft Hydraulics, Report No. H840.32, Part 2. 1994. Available online: https://repository.tudelft.nl/islandora/object/uuid:181f5c9c-f66f-4103-abe6-c4633137e351 (accessed on 10 January 2021).
  42. Fredsøe, J.; Deigaard, R. Mechanics of Coastal Sediment Transport; World Scientific: Singapore, 1992; Volume 3. [Google Scholar]
  43. Boers, M. Surf Zone Turbulence. 2005. Available online: https://repository.tudelft.nl/islandora/object/uuid:71428ad5-0658-49f7-b83e-0deb26b668cf (accessed on 10 January 2021).
  44. Hamilton, D.G.; Ebersole, B.A. Establishing uniform longshore currents in a large-scale sediment transport facility. Coast. Eng. 2001, 42, 199–218. [Google Scholar] [CrossRef]
  45. Teles, M.J.; Pires-Silva, A.A.; Benoit, M. Three dimensional modelling of wave-current interactions in the coastal zone with a fully coupled system. Geography. 2013. Available online: https://www.researchgate.net/profile/Maria_Teles2/publication/315066486_THREE_DIMENSIONAL_MODELLING_OF_WAVE-CURRENT_INTERACTIONS_IN_THE_COASTAL_ZONE_WITH_A_FULLY_COUPLED_SYSTEM/links/58c9539f4585157512329679/THREE-DIMENSIONAL-MODELLING-OF-WAVE-CURRENT-INTERACTIONS-IN-THE-COASTAL-ZONE-WITH-A-FULLY-COUPLED-SYSTEM.pdf (accessed on 10 January 2021).
Figure 1. A flowchart describing the methodology of the research.
Figure 1. A flowchart describing the methodology of the research.
Jmse 09 00076 g001
Figure 2. Bathymetry of the computational area.
Figure 2. Bathymetry of the computational area.
Jmse 09 00076 g002
Figure 3. Spatial distribution of Stokes drift.
Figure 3. Spatial distribution of Stokes drift.
Jmse 09 00076 g003
Figure 4. Distribution of mean water level ζ ¯ .
Figure 4. Distribution of mean water level ζ ¯ .
Jmse 09 00076 g004
Figure 5. Distribution of F w , 1 / ρ ( m / s 2 ).
Figure 5. Distribution of F w , 1 / ρ ( m / s 2 ).
Jmse 09 00076 g005
Figure 6. Vertical distribution of the generalized Lagrangian mean (GLM) velocity and quasi-Eulerian mean velocity.
Figure 6. Vertical distribution of the generalized Lagrangian mean (GLM) velocity and quasi-Eulerian mean velocity.
Jmse 09 00076 g006
Figure 7. Vertical distribution of the Eulerian-mean velocity.
Figure 7. Vertical distribution of the Eulerian-mean velocity.
Jmse 09 00076 g007
Figure 8. Vertical distribution of Reynolds shear stress u t w t ¯ in the CO condition.
Figure 8. Vertical distribution of Reynolds shear stress u t w t ¯ in the CO condition.
Jmse 09 00076 g008
Figure 9. Vertical distribution of turbulent viscosity ( υ v T ).
Figure 9. Vertical distribution of turbulent viscosity ( υ v T ).
Jmse 09 00076 g009
Figure 10. Vertical distribution of the wave radiation stress component u ˜ w ˜ ¯ C S .
Figure 10. Vertical distribution of the wave radiation stress component u ˜ w ˜ ¯ C S .
Jmse 09 00076 g010
Figure 11. Vertical distribution of wave-induced forcing terms.
Figure 11. Vertical distribution of wave-induced forcing terms.
Jmse 09 00076 g011
Figure 12. Reynolds turbulent stress u t w t ¯ in the CO condition.
Figure 12. Reynolds turbulent stress u t w t ¯ in the CO condition.
Jmse 09 00076 g012
Figure 13. Spatial distribution of quasi-Eulerian mean velocity in the CO condition.
Figure 13. Spatial distribution of quasi-Eulerian mean velocity in the CO condition.
Jmse 09 00076 g013
Figure 14. Vertical distribution of quasi-Eulerian mean velocity in the CO condition.
Figure 14. Vertical distribution of quasi-Eulerian mean velocity in the CO condition.
Jmse 09 00076 g014
Figure 15. Spatial distribution of mean velocity fields in the WO condition.
Figure 15. Spatial distribution of mean velocity fields in the WO condition.
Jmse 09 00076 g015
Figure 16. Vertical distribution of quasi-Eulerian mean velocity in the WO condition.
Figure 16. Vertical distribution of quasi-Eulerian mean velocity in the WO condition.
Jmse 09 00076 g016
Figure 17. Spatial distribution of quasi-Eulerian mean velocity.
Figure 17. Spatial distribution of quasi-Eulerian mean velocity.
Jmse 09 00076 g017
Figure 18. Vertical distribution of quasi-Eulerian mean velocity in the linear scale.
Figure 18. Vertical distribution of quasi-Eulerian mean velocity in the linear scale.
Jmse 09 00076 g018
Figure 19. Vertical distribution of quasi-Eulerian mean velocity in the semilogarithmic scale.
Figure 19. Vertical distribution of quasi-Eulerian mean velocity in the semilogarithmic scale.
Jmse 09 00076 g019
Figure 20. The bottom profile.
Figure 20. The bottom profile.
Jmse 09 00076 g020
Figure 21. Distribution of significant wave height H S .
Figure 21. Distribution of significant wave height H S .
Jmse 09 00076 g021
Figure 22. Distribution of wave-induced forcing terms at x = 22.9   m .
Figure 22. Distribution of wave-induced forcing terms at x = 22.9   m .
Jmse 09 00076 g022
Figure 23. Distribution of mean water level ζ ¯ .
Figure 23. Distribution of mean water level ζ ¯ .
Jmse 09 00076 g023
Figure 24. Distribution of Stokes drift (a,b) and quasi-Eulerian mean velocity (c,d).
Figure 24. Distribution of Stokes drift (a,b) and quasi-Eulerian mean velocity (c,d).
Jmse 09 00076 g024aJmse 09 00076 g024b
Figure 25. Vertical distribution of horizontal mean velocity in test 1B.
Figure 25. Vertical distribution of horizontal mean velocity in test 1B.
Jmse 09 00076 g025aJmse 09 00076 g025b
Figure 26. Vertical distribution of horizontal mean velocity in test 1C.
Figure 26. Vertical distribution of horizontal mean velocity in test 1C.
Jmse 09 00076 g026aJmse 09 00076 g026bJmse 09 00076 g026c
Figure 27. Distribution of significant wave height H S .
Figure 27. Distribution of significant wave height H S .
Jmse 09 00076 g027
Figure 28. Distribution of the mean water level ζ ¯ .
Figure 28. Distribution of the mean water level ζ ¯ .
Jmse 09 00076 g028
Figure 29. Vertical distribution of wave-induced forcing terms.
Figure 29. Vertical distribution of wave-induced forcing terms.
Jmse 09 00076 g029
Figure 30. Distribution of cross-shore velocity.
Figure 30. Distribution of cross-shore velocity.
Jmse 09 00076 g030
Figure 31. Vertical profiles of cross-shore velocity.
Figure 31. Vertical profiles of cross-shore velocity.
Jmse 09 00076 g031
Figure 32. Vertical distribution of wave-induced forcing terms.
Figure 32. Vertical distribution of wave-induced forcing terms.
Jmse 09 00076 g032
Figure 33. Distribution of the longshore velocity.
Figure 33. Distribution of the longshore velocity.
Jmse 09 00076 g033
Figure 34. Longshore velocity at z = 2 h / 3 .
Figure 34. Longshore velocity at z = 2 h / 3 .
Jmse 09 00076 g034
Figure 35. Vertical profiles of the longshore velocity at different locations.
Figure 35. Vertical profiles of the longshore velocity at different locations.
Jmse 09 00076 g035aJmse 09 00076 g035b
Table 1. Wave properties at the paddle.
Table 1. Wave properties at the paddle.
Wave TypeTp (s)Hrms (m)h (m)
Random1.70.10.5
Table 2. Bottom boundary thickness δ in different waves and current conditions.
Table 2. Bottom boundary thickness δ in different waves and current conditions.
ConditionsFormula (62)
( × 10 3   m )
Van Rijn [35]
( × 10 3   m )
CO0.10.1
WO1.31.3
WFC5.33.6
WOC4.93.6
Table 3. Characteristics of the near-bed mean flow.
Table 3. Characteristics of the near-bed mean flow.
Conditions | u ¯ b |
( × 10 2   m / s )
| τ ¯ b |
( × 10 2   kgm 2 / s 2 )
u
( × 10 2   m / s )
CO8.165.40.74
WO0.90.210.13
WFC8.1037.671.94
WOC5.8624.61.56
Table 4. Wave properties at the offshore boundary.
Table 4. Wave properties at the offshore boundary.
ExperimentHs (m)Tp (s)
Test 1B0.2062.03
Test 1C0.1033.3.3
Table 5. Incident wave properties at the boundary.
Table 5. Incident wave properties at the boundary.
Wave TypeTp (s)Hs (m) θ ( 0 ) h (m)
Irregular2.50.225100.667
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nguyen, D.T.; Jacobsen, N.G.; Roelvink, D. Development and Validation of Quasi-Eulerian Mean Three-Dimensional Equations of Motion Using the Generalized Lagrangian Mean Method. J. Mar. Sci. Eng. 2021, 9, 76. https://doi.org/10.3390/jmse9010076

AMA Style

Nguyen DT, Jacobsen NG, Roelvink D. Development and Validation of Quasi-Eulerian Mean Three-Dimensional Equations of Motion Using the Generalized Lagrangian Mean Method. Journal of Marine Science and Engineering. 2021; 9(1):76. https://doi.org/10.3390/jmse9010076

Chicago/Turabian Style

Nguyen, Duoc Tan, Niels G. Jacobsen, and Dano Roelvink. 2021. "Development and Validation of Quasi-Eulerian Mean Three-Dimensional Equations of Motion Using the Generalized Lagrangian Mean Method" Journal of Marine Science and Engineering 9, no. 1: 76. https://doi.org/10.3390/jmse9010076

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop