The evaluation and consideration of the mean flow in wave evolution equations are necessary for the accurate prediction of fluid particle trajectories under wave groups, with relevant implications in several domains, from the transport of pollutants in the ocean to the estimation of energy and momentum exchanges between the waves at small scales and the ocean circulation at large scale. We derive an expression of the mean flow at a finite water depth, which, in contrast to other approximations in the literature, accurately accords with the deep-water limit at third order in steepness and is equivalent to second-order formulations in intermediate water. We also provide envelope evolution equations at fourth order in steepness for the propagation of unidirectional wave groups either in time or space that include the respective mean flow term. The latter, in particular, is required for accurately modeling experiments in water wave flumes in arbitrary depths.

## I. INTRODUCTION

As waves evolve on the ocean surface, they induce a mean flow to the fluid particles, particularly when nonlinear effects are taken into account. Near the surface, fluid particles experience a net horizontal displacement in the same direction as the water wave propagation, the so-called Stokes drift,^{1,2} that decays with depth. A return flow in the fluid column along vertical and horizontal directions guarantees the conservation of the water-mass transport, causing localized variations in the mean water level under wave groups and the associated propagation of infragravity waves^{3} at a finite water depth. An accurate description of the mean flow is, thus, necessary for the proper reconstruction of the fluid particle trajectories underneath water waves.^{4,5} Since the wave-induced mean flow is associated with the transport of energy, momentum, and other tracers,^{6} such as pollutants like plastic^{7} or offshore oil spill,^{8} it is relevant for environmental studies. Moreover, it has an impact on the general circulation of the ocean at large scales and its modeling,^{9–12} and on the nearshore circulation,^{3,13} in particular in shoaling regions,^{14,15} while it is affected by the presence of background shear currents.^{16,17}

Since Stokes drift and return flow are two phenomena occurring at second order in steepness, they are taken into account in the finite water depth nonlinear Schrödinger equation (NLS), which describes the evolution of the envelope of narrow-banded wave packets and can be obtained using a multi-scale development of water surface elevation and velocity potential at third order in steepness,^{18–21} or taking the narrow-banded limit of the Zakharov equation using an Hamiltonian approach.^{22,23} Indeed, at third order in steepness, it was shown^{20} that the mean flow term comes into play as a modification of the nonlinear coefficient, giving rise to a transition from focusing to defocusing regimes at the critical value $ k 0 h \u2248 1.363$, where *k*_{0} is the carrier wavenumber and *h* the depth. However, this contribution to the nonlinear coefficient disappears in the deep-water limit. At fourth order in steepness, the accurate formulation of the mean flow term needs to be accounted for, in both, finite and infinite depth waters.^{24}

Fourth-order terms in the NLS equation are necessary to explain features like asymmetrical evolution of spectra^{25} and asymmetries in the waveform,^{26–30} especially when inevitable wave focusing is at play.^{31} Several versions of high-order wave envelope evolution equations exist in the literature, which can be obtained when applying the multiple scales development using steepness and bandwidth with different orders of magnitude^{24,32–34} or the Hamiltonian approach.^{23,35–38} The main difference between all these equations is the treatment and approximation of the mean flow term.

Here, we will focus on narrow-banded unidirectional wave packets where steepness and bandwidth can be considered as parameters of the same order, like in the finite-depth developments described in Refs. 39 (denoted as Sed03 in the following) and 40 (Slu05). The challenge of such developments is the fact that they do not reduce to the Dysthe equation (Ref. 24, Dys79) in the deep-water limit, which has been shown to well reproduce wave tank experiments.^{28} We will show that this convergence depends on how the mean flow is approximated. Moreover, we will provide finite-depth envelope equations at fourth order in steepness in both space-like and time-like formulations, with correct limiting expressions in deep water. In particular, such unidirectional time-like equations that allow a continuous scaling from intermediate to the deep water limit are relevant for reproducing wave tank experiments with high accuracy in arbitrary depths. In contrast, directional sea states are typically found in the open ocean, and this restricts the applicability of the above unidirectional modeling equations.

The paper is organized as follows: In Sec. II, we will derive the expression of the mean flow to be inserted in the envelope equation at fourth order in steepness for the evolution in time. We will compare this expression with the ones already existing in the literature and show that it indeed correctly describes the mean flow in the whole range of water depths, i.e., from intermediate to deep water regimes. In Sec. III, we will provide an envelope equation for the water wave evolution in space, relevant for modeling, for instance water tank experiments, and the corresponding mean flow term, at fourth order in steepness. Again, we will perform the comparison for various water depth scenarios. Finally, in Sec. IV, we will summarize our findings.

## II. FOURTH-ORDER EQUATION: PROPAGATION IN TIME

*U*of a unidirectional progressive gravity wave packet propagating in time on the free surface of a homogeneous liquid with depth

*h*,

^{39–41}

^{,}

^{42}$ \eta ( X , t ) = 1 2 [ U ( X , t ) \u2009 exp \u2009 ( i ( k 0 X \u2212 \omega 0 t ) ) + c . \u2009 c . ]$, where

*ω*

_{0}and

*k*

_{0}are the angular frequency and wavenumber of the carrier wave, respectively.

^{39}in a reference frame moving with the group velocity $ x = X \u2212 c g t$,

*μ*is given in Eq. (A15), and the high-order nonlinear coefficients $ Q \u0303 41$ and $ Q \u0303 42$ in Eqs. (B4) and (B7), respectively. The nonlinear term $ \beta \u0302 D$ is a positive function (see Fig. 1), given in Eq. (A9), and since the dispersion coefficient $ \alpha \u0302$ is negative, it seems that the characterization of the focusing and defocusing regime is somehow lost in the present formulation.

_{g}*ε*, a dummy variable which is set to 1 at the end, that is helpful for grouping terms of the same order in steepness. The last term, depending on the zero harmonic of the velocity potential $ \varphi 0$, is the mean flow term that, from the multi-scale development, takes the following form [see Eqs. (35) and (57) in Sed03]:

*ν*given in Eq. (A16) and $ q \u0303 40 S$ in Eq. (B10). Its contribution can be included into the coefficients $ \beta \u0302 D , \u2009 Q \u0303 41 , \u2009 Q \u0303 42$ to obtain $ \beta \u0302 , \u2009 \beta \u0302 21 , \u2009 \beta \u0302 22$, respectively, in Eq. (1) [see steps from Eqs. (40) to (42) and from Eqs. (64) to (66) in Sed03]. Thus, the mean flow term is already taken into account in Eq. (1) at fourth order in steepness through the nonlinear coefficients.

We remark that from Eq. (3), in the deep-water limit, $ \u2202 \varphi 0 / \u2202 x \u2192 0$ since $ \nu \u2192 \u2212 \u221e$, while all the other coefficients are finite. This agrees with the developments at third order in steepness.^{20,43} Note that in the 1D shallow-water case, the mean flow term modifies the nonlinear term, see Eq. (2.21) in Ref. 20, while its contribution vanishes in the deep-water limit. However, the fact that the contribution of mean flow remains zero in deep water at higher-order approximation in steepness is not ideal, as we will discuss in the following subsection.

### A. Mean flow in the deep-water limit

^{24,32,44}and time-like 1D formulations.

^{28,45–48}Let us consider for the moment the propagation in time. The Dysthe equation reads

^{28,47}

*η*being a function of

*x*. The role of the mean flow Hilbert term in the evolution of pulsating wave packets

^{25}and narrow-banded irregular waves

^{27,28}has been shown in several experiments in deep water, and its presence is required for a correct modeling.

^{26}

^{32}suggested that the system can be closed by solving the equations for $ \varphi 0$ as a function of the envelope

*U*in the Fourier space. It is instructive to report here the explicit derivation of the Hilbert term since we will use analogous developments in the following. The zero harmonic of the velocity potential $ \varphi 0$ satisfies the Laplace equation in the entire water column with boundary conditions at the surface and at the bottom, giving a set of equations that constitutes the following Neumann problem:

*C*

_{1}independent of

*z*. This satisfies the bottom boundary condition for $ z \u2192 \u2212 \u221e$.

*C*

_{1}is obtained, and thus,

*x*at

*z*= 0 in the Fourier space is given by

### B. Multi-scale development and Hilbert term

^{24,32}Thus, there is the need to reconcile these results. This can be done as follows. In the derivation of the Hilbert term, we use the Laplace equation for $ \varphi 0$ that is the complete mean flow, and not just its approximation at second order in steepness as in Eq. (12). Indeed, the Laplace equation at third order for the mean flow is

*z*, this gives

*C*defined in Eq. (A17) and in Ref. 4, where $ D \u2032 = ( 1 + C F D ) \u2009 \omega 0 / ( 2 \sigma ) = \mu g \omega 0 / ( 8 \sigma 2 )$ (see definitions in Appendix A). Using Eq. (14) at

_{FD}*z*= 0, Eq. (15) reduces to

*D*is defined in Eq. (A13). Equation (16) can also be written, using Eq. (14), as

Thus, in the derivation of the Hilbert term, the complete mean flow (in the Laplace equation) is considered together with the mean flow at third order in steepness in the surface boundary condition. Consequently, a “hybrid” relation as described by Eq. (5) is obtained, which is different from Eq. (12) where only the first terms in the development of the mean flow at the surface are taken into account. In other words, the expression in Eq. (5) is inherently nonlocal. The same occurs in the nonlinear terms of the supercompact model,^{49} from which the Dysthe equation can be derived. The terms in Sed03 and Slu05 are instead local, as they only involve the surface, and not the entire water column. This is obviously an approximation, since the mean flow does involve a body of fluid, which is not immediately at the surface, as clearly shown in field measurements.^{3}

### C. Mean flow term in arbitrary depth

We now repeat the procedure used in Sec. II A for the case of intermediate water. We will consider two cases that differ based on the considered surface boundary condition: The Neumann problem is solved using the condition given by Eq. (17) in case 1, using Eq. (15) in case 2. Replacing the expression of $ \u2202 \varphi 0 / \u2202 x$ that is obtained in each case into the last term of Eq. (2) gives the final high-order NLS equation in an arbitrary finite depth and in the space-like form.

#### 1. Case 1

*C*and therefore,

*z*= 0, this gives

*x*at

*z*= 0 is given by

#### 2. Case 2

The surface boundary condition is now Eq. (15). Note that the relation $ \varphi 03 z = \u2212 h \varphi 01 x x$ has not been used to simplify the lhs of this equation, and both mean flow terms are, thus, considered being of the same order, see Eq. (13) in Ref. 4.

*C*(

*k*,

*t*):

*c*is the wave group speed at the carrier wavenumber. Performing analogous steps as in the previous case, the final expression for the Euler horizontal velocity in

_{g}*z*= 0 is

*x*on the rhs of Eqs. (25) or (28), considering that $ D = D \u2032 / ( 1 \u2212 c g 2 / ( g h ) )$ and using $ coth y = 1 / y + O ( y )$, we get in both cases 1 and 2 for small

*kh*numbers,

### D. Numerical comparisons

We compare the expressions for the horizontal velocity $ \u2202 \varphi 0 / \u2202 x$ listed in Table I with the sub-harmonic velocity potential $ \varphi 20$ at second order in steepness and its horizontal derivative calculated using the Dalzell analytical method [see the original paper, Ref. 50 (Dal99), and the explicit formulas reported in the Appendix of Ref. 51].

Case 1 | Eq. (25) | Third-order, nonlocal |

Case 2 | Eq. (28) | Third-order, nonlocal |

Dys79 (Ref. 24) | Eq. (5) | Third-order, nonlocal deep water |

Sed03 (Ref. 39) | Eq. (3) | Third-order |

Sed03 | Eq. (29) | Second-order |

Slu05 (Ref. 40) | Eq. (C10) | Third-order |

For benchmarking and validation purposes, we use the same parameters as in Ref. 51 (case C in their Table II), namely, a Gaussian (amplitude) spectrum *S*(*k*) with peak wavenumber $ k 0 = 0.0277$ m^{−1}, wavelength $ \lambda 0 = 2 \pi / k 0$, standard deviation of the spectrum given by symmetrical values $ k w = k w 1 = k w 2 = 0.27 k 0$, steepness $ \epsilon = 0.3$, and different values of the normalized depth $ k 0 h$ in the range $ [ 0.5 , 50 ]$, thus the case of finite and infinite water depth conditions. The angular frequency is calculated from $ \omega 0 2 = g k 0 tanh ( k 0 h )$. In particular, the surface elevation at first order in steepness given by the superposition of *N* = 30 waves is used to calculate the intensity of wave envelope, i.e., $ | U | 2$, and then, the horizontal velocity $ \u2202 \varphi 0 / \u2202 x$ through the different relations, as listed in Table I. We consider a focused wave group, composed of *N* individual sinusoidal wave components that are in phase at a single point in time and space^{52} (the origin in our case), and a random sea state by using a uniform distribution for the phases. An example of sea state realization is given in Appendix D.

Case 1 | Eq. (39) | Third-order, non-instantaneous |

Case 2 | Eq. (40) | Third-order, non-inst. |

Dys79 (Ref. 24) | Eq. (36) | Third-order, non-inst. deep water |

Sed03 (Ref. 39) | Eq. (38) | Third-order |

Sed03 | 1st term on rhs of Eq. (38) | Second-order |

Slu05 (Ref. 40) | Eq. (C11) | Third-order |

The results are summarized in Fig. 2 for waves focusing at *x* = 0, and in Fig. 3 for the case of waves superposition with random phases. In both cases, we see that the horizontal velocity calculated by Eq. (29) (gray line), which corresponds to the second-order expression in Sed03, reproduces well the Dalzell waveform (green line) only for $ k 0 h \u2264 2$ and cannot be applied in deep water where it gives a waveform amplitude much smaller than Dysthe's result (black line). On the contrary, Dysthe's expression should not be used for $ k 0 h < 10$ as in our considered cases, especially for large amplitude waveforms, where it does not attain Dalzell's accuracy. Although an elementary consideration of the dispersion relation suggests that $ k 0 h \u2248 5$ should be indistinguishable from $ k 0 h \u2192 \u221e$, the dynamics of the mean flow is consistently different from this limit up to $ k 0 h = 10$. The expressions for the horizontal velocity given by case 1 [Eq. (25)] and case 2 [Eq. (28)] behave in a similar way at all depths, providing in general an approximation that corresponds well to the Dalzell solution in intermediate waters, and to Dysthe in deep waters. It is also important to note that the expressions provided in Sed03 and Slu05 at third order in steepness [Eqs. (3) and (C10), respectively], give different results and do not correspond to the other models at any depth.

More detailed features can be deduced from Fig. 4, where the deviation operator with respect to the Dysthe expression, defined as $ D Dysthe = N \u2212 1 \u222b ( \varphi 0 x \u2212 \varphi 0 x Dysthe ) 2 \u2009 d x$, and to the Dalzell one, $ D Dalzell = N \u2212 1 \u222b ( \varphi 0 x \u2212 \varphi 0 x Dalzell ) 2 \u2009 d x$, with $ \varphi 0 x = \u2202 \varphi 0 / \u2202 x$, is shown as a function of the non-dimensional water depth $ k 0 h$. The integral is calculated over $ L = 80 \lambda 0$, and the normalization coefficient is $ N = ( \omega 0 D W \lambda 0 ) 2 L$, with $ \omega 0 D W$ calculated in deep water (DW). The Stokes series expansion of the velocity potential (and thus Dalzell's model) converges in shallow water^{53} if $ 3 \epsilon / ( 2 k 0 h ) 3 \u226a 1$, thus at low Ursell number, which implies in our case $ k 0 h \u226b 0.48$. This region is excluded in Fig. 4. It can be seen that the expression for the horizontal velocity given by case 2 [Eq. (28)] is accurate at second order at all depths, almost superposing to the second-order Dalzell solution, while the one given by case 1 [Eq. (25)] is the only expression that consistently converges to the Dysthe in the deep-water limit expression and is equivalent to case 2 for $ k 0 h < 5$. Note that the third-order corrections provided in Sed03 [Eq. (3)] and in Slu05 [Eq. (C10)] are different and do not extend the validity to deeper water regimes of the second-order expression [Eq. (29)], which is accurate for $ k 0 h < 2$ [see Fig. 4(b)]. This raises a warning on the quantitative accuracy of these third-order expressions.

## III. FOURTH-ORDER EQUATION: PROPAGATION IN SPACE

^{54}are easily included up to arbitrary order following Refs. 55 and 56, reducing the constraints in the bandwidth of the original Dysthe equation.

^{39}

*α*,

*α*

_{3}given in Appendix A and $ Q \u0303 41 , \u2009 Q \u0303 42$ in Appendix B. The dependence of the nonlinear coefficients on $ k 0 h$ is illustrated in Fig. 6. Note that $ \beta , \beta 21 , \beta 22 , B 21$, and $ B 22$ diverge quite strongly for $ k 0 h \u2192 0$, i.e., particularly in the defocusing regime. That said and as can be also analytically verified, all coefficients correctly reproduce their deep-water limit at $ k 0 h \u2192 \u221e$.

### A. Mean flow in time-like equations

*x*and

*t*separately, but only through their combination $ ( x \u2212 c g t )$, so that at the leading order, the Fourier transform in space, $ F x$, can be replaced by the Fourier transform in time, $ F t$. From Eq. (5), in the deep-water limit, it is possible to simply exchange the derivative with respect to space with that with respect to time and vice versa, since the Hilbert transform of the derivative is the derivative of the Hilbert transform, i.e., these two linear operators commute. As such,

^{28,45–48}

### B. Numerical comparisons

We now compare the expressions for $ \u2202 \varphi 0 / \u2202 t$ listed in Table II with the sub-harmonic velocity potential $ \varphi 20$ at second order in steepness and its time derivative calculated using the Dalzell analytical method in Refs. 50 and 51.

We use the same parameters as in Sec. II D, starting in this case by a Gaussian (amplitude) spectrum $ S ( \omega )$ peaked at $ \omega 0 = 1$ Hz.

The results for *N* = 30 waves focusing at *t* = 0 are shown in Fig. 7. We see that the second-order expression in Sed03 reproduces well the waveform only for $ k 0 h \u2264 1$, showing discrepancies with respect to Dalzell's solution at $ k 0 h = 2$ and not converging to Dysthe's solution in deep water. As for the space-like form, the expressions for $ \varphi 0 t$ given by case 1 [Eq. (39)] and case 2 [Eq. (40)] have a similar behavior at all depths, providing in general approximations that are as accurate as the Dalzell solution in intermediate waters, and as the Dysthe expression at deep waters. As before, the third-order expressions provided by Eqs. (38) and (C11) are different and do not correspond to the other models. Figure 8 shows the deviation operator with respect to the Dysthe expression, defined similarly to the spatial counterpart (see Sec. II D) as $ D Dysthe = N \u2212 1 \u222b ( \varphi 0 t \u2212 \varphi 0 t Dysthe ) 2 \u2009 d t$, and to the Dalzell one, $ D Dalzell = N \u2212 1 \u222b ( \varphi 0 t \u2212 \varphi 0 t Dalzell ) 2 \u2009 d t$ as a function of the non-dimensional water depth $ k 0 h$. The integral is now calculated over $ T = 80 \u2009 T 0$, and the normalization coefficient is $ N = ( \omega 0 \lambda 0 D W ) 4 T$, with $ \lambda 0 D W$ calculated at deep water. As for the space-like case, it can be seen that the expression given by case 2 [Eq. (40)] is accurate at second order at all depths as the second-order Dalzell solution, while the one given by case 1 [Eq. (39)] is the only expression that converges to the Dysthe term in the deep-water limit and is equivalent to case 2 for $ k 0 h < 3$.

## IV. CONCLUSION

During the evolution of surface gravity waves, fluid particles experience Stokes drift along the propagation direction of the waves and a return flow in vertical and horizontal directions, which closes the water-mass transport. Both processes give rise to a mean flow at the surface that needs to be accounted for in order to accurately predict the transport of pollutants^{6–8} (such as oil and microplastics) and, more in general, the impact of waves evolution at small scale on the ocean circulation at large scale.^{9–12} We have derived nonlocal (viz. non-instantaneous) expressions of the mean flow [Eqs. (25) and (39)] that correctly converge to the deep-water limit at third order in steepness, while being equivalent to second-order formulations in intermediate waters. We have included these expressions in an envelope evolution equation at fourth order in steepness in both, space-like [Eq. (2)] and time-like [Eq. (34)] formulations. We emphasize that the time-like form is relevant to study the evolution of unidirectional wave groups in space, thus for modeling experiments in water wave flumes at high accuracy in arbitrary depths. Future work will be focusing on the experimental validation of our results in different water depth regimes.

## ACKNOWLEDGMENTS

A.G., J.K., and M.B. acknowledge financial support from the Swiss National Science Foundation (Project No. 200020‐175697). D.E. acknowledges financial support from the Swiss National Science Foundation (Fellowship No. P2GEP2‐191480). A.C. acknowledges support from Kyoto University's Hakubi Center for Advanced Research.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Alexis Gomel:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (supporting); Writing – review & editing (equal). **Corentin Montessuit:** Formal analysis (equal); Investigation (equal). **Andrea Armaroli:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – review & editing (equal). **Debbie Eeltink:** Formal analysis (equal); Writing – review & editing (equal). **Amin Chabchoub:** Supervision (equal); Writing – review & editing (equal). **Jerome Kasparian:** Funding acquisition (lead); Supervision (equal); Writing – review & editing (equal). **Maura Brunetti:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (lead); Supervision (equal); Writing – original draft (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study and the MATLAB scripts are available from the corresponding author upon reasonable request.

### APPENDIX A: DISPERSION AND NONLINEAR COEFFICIENTS

*κ*, $ c g / c p = ( \sigma + \kappa ( 1 \u2212 \sigma 2 ) ) / ( 2 \sigma ) , \u2009 g h / c g 2 = ( c p 2 / c g 2 ) \kappa / \sigma , \u2009 2 \sigma 2 c g 2 / ( g h \u2212 c g 2 ) = 2 \sigma ( c g 2 / c p 2 ) / ( \kappa / \sigma \u2212 c g 2 / c p 2 )$ and

### APPENDIX B: NOTATION USED IN SEDLETSKY (2003) (Sed03)

^{39},

*Q*

_{41}and

*Q*

_{42}are given in Eqs. (67) and (68) in Ref. 39 and reported here for completeness,

From Eqs. (B8) and (B9), in the deep-water limit $ \kappa \u2192 \u221e , \u2009 \sigma \u2192 1 , \u2009 \nu \u2192 1 \u2212 4 \kappa , \u2009 \mu g \u2192 4 , \u2009 Q 41 \u2192 768 / ( 32 \xb7 16 ) = 3 / 2$, and $ Q 42 \u2192 128 / ( 32 \xb7 16 ) = 1 / 4$ recovering the Dysthe result for such terms.

### APPENDIX C: NOTATION USED IN SLUNYAEV (2005) (SLU05)

^{57},

### APPENDIX D: EXAMPLE OF SEA STATE REALIZATION

*S*(

*k*), given by

*k*> 0, where $ k 0 = 0.0277$ m

^{−1}, and $ k w = 0.27 \u2009 k 0$ is the dimensional bandwidth. Other spectra, like JONSWAP or Pierson–Moskowitz, can also be used. This yields the following surface elevation at first order in steepness:

*A*,

_{p}*x*, and

_{f}*t*are the amplitude, position, and time for the group at linear focus, with steepness given by $ \epsilon = A p k 0 = 0.3$. An example of sea state realization at second order in steepness, obtained using the Dalzell development,

_{f}^{50}is shown in Fig. 9 for the case of a focused wave group at

*x*= 0,

_{f}*t*= 0 for $ k 0 h = 1.5$. The power spectrum can be obtained as

_{f}^{58}$ P ( k ) = | \eta \u0302 | 2 / ( 2 d k )$, where $ \eta \u0302$ is the Fourier transform in space of the surface elevation, from which the significant wave height

*H*can be obtained as $ H s = 4 m 0 = 5.2$ m,

_{s}*m*

_{0}being equal to the area under the power spectrum curve.

## REFERENCES

*Theory and Applications of Ocean Surface Waves*

*Ocean Wave Dynamics*