# A Weakly Nonlinear System for Waves and Sheared Currents over Variable Bathymetry

Julien Touboul, Veronica Morales-Marquez, Kostas Belibassakis

► To cite this version:

Julien Touboul, Veronica Morales-Marquez, Kostas Belibassakis. A Weakly Nonlinear System for Waves and Sheared Currents over Variable Bathymetry. Journal of Marine Science and Engineering, 2024, 12 (3), pp.509. 10.3390/jmse12030509. [https://doi.org/10.3390/jmse12030509](https://doi.org/10.3390/jmse12030509)

HAL Id: hal-04746284

https://hal.science/hal-04746284v1

Submitted on 24 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Weakly Nonlinear System for Waves and Sheared Currents over Variable Bathymetry

Julien Touboul<fnref n="1" />,*, Veronica Morales-Marquez<fnref n="2" /> and Kostas Belibassakis<fnref n="3" />

**Abstract:** The wave-current-seabed interaction problem is studied by using a coupled-mode system developed for modeling wave scattering by non-homogeneous, sheared currents in variable bathymetry regions. The model is based on a modal series expansion of wave velocity based on vertical eigenfunctions, dependent on local depth and flow parameters, including propagating and evanescent modes. The latter representation is able to accurately satisfy the wave flow continuity condition and the no-entrance boundary condition on the sloping parts of the seabed. A new derivation of a simplified nonlinear system is introduced using decomposition to a mean flow and a perturbative wave field. To force the system to consider incoming waves at the inlet, boundary knowledge of periodic, travelling nonlinear water waves over a flat bottom is required. For this purpose, specific solutions are derived using the semi-analytical method based on the stream function formulation, for cases of water waves propagating above linearly and exponentially sheared currents. Results obtained by the application of the CMS concerning the propagation of waves and currents—in particular, examples characterized by depth inhomogeneities—are presented and discussed, illustrating the applicability and performance of the method.

**Keywords:** wave-current interaction; sheared current; variable bathymetry; coupled modes

## 1. Introduction

The problem of water propagating on top of vortical flows has attracted the attention of many scientists within the water wave community. Indeed, it has recently been shown that waves interacting with an opposing shear current over variable bathymetry can have a significant effect on the focusing of the wave field; see Rey et al. [1]. Such situations can be observed in natural conditions, when strong currents interact with the bathymetry [2,3], or in cases of strong wind forcing on the ocean-atmosphere interface; see, e.g., [4], and the recent review by Zhan et al. [5]. The propagation of periodic waves in the presence of a uniform shear background current has been studied by using complete description of the flow in the framework of Euler equations by Nwogu [6] and by Zhang et al. [7] using RANS equations.

For both coastal safety and engineering purposes, including the estimation of loads on structures in nearshore and coastal regions, the development and application of efficient models for predicting the nonlinear, phase-resolved behavior of water waves and currents in such areas are important. Results in this direction, for several case studies, have been presented by various authors; see, e.g., Refs. [8–11]. In particular, wave- and current-dominated orthogonal flows over fixed rough beds have recently been considered by Faraci et al. [12], where a review of coastal research also focused on wave-current flows, either collinear or at an angle can be found.

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Touboul, J.; Morales-Marquez, V.; Belibassakis, K. A Weakly Nonlinear System for Waves and Sheared Currents over Variable Bathymetry. *J. Mar. Sci. Eng.* **2024**, *12*, 509. https://doi.org/10.3390/jmse12030509

Academic Editor: Anders Jensen
Knudby

Received: 14 February 2024

Revised: 13 March 2024

Accepted: 17 March 2024

Published: 19 March 2024

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">1 École Centrale Méditerranée, Aix Marseille Université, CNRS, IRPHE UMR 7342, 13384 Marseille, France; julien.touboul@centrale-med.fr</footnote>

<footnote n="2">2 University of Toulon, Aix Marseille University, CNRS/INSU, IRD, MIO UM 110, Mediterranean Institute of Oceanography, 83130 La Garde, France; veronica.morales-marquez@mio.osupytheas.fr</footnote>

<footnote n="3">3 School of Naval Architecture and Marine Engineering, National Technical University of Athens, Zografos, 15773 Athens, Greece</footnote>

<footnote>* Correspondence: kbel@fluid.mech.ntua.gr</footnote>
It is remarked that in many cases of wave-current-seabed interaction, computational difficulties are encountered in the analysis, since large areas often need to be considered, and the evolution of water waves propagating in such inhomogeneous media involves multiple-scale phenomena. In such cases, simplified models can be derived, allowing to separate the scales, solving the background flow on the one hand and the perturbative flow corresponding to water waves on the other hand. The latter can be solved by expanding its vertical structure on a selected basis consisting of modes and applying a Galerkin projection method to derive coupled equations with respect to the horizontal mode amplitudes. Indeed, by applying such a technique, a significant reduction in unknowns can be obtained, depending on the truncation made for the expansion. This approach has been widely studied in the absence of currents [13–17] or in the presence of vertically uniform currents, as in, e.g., [17–19]. When it comes to flows interacting with vorticity, Yang and Liu [20] recently introduced a Boussinesq-type formulation, projecting the velocity field on a polynomial basis. Also, an analysis of solitary wave and linear shear current interaction by using Green–Naghdi (GN) equations was described by Duan et al. [21].

Moreover, coupled-mode systems have been derived, involving propagating and evanescent modes, which are suitable to describe wave dispersion in intermediate water depth. Such models are applied to the scattering of water waves by horizontally shearing currents in variable bathymetry. In the case of wave propagation in the presence of vertically sheared currents, a generalized mild-slope equation was derived by Touboul et al. [22] by retaining only the propagation mode in the expansion, with a restriction to vertically linear currents. Extensions of the above approach to multi-mode systems for waves propagating over arbitrarily sheared currents were presented in [23], but they were essentially restricted to the framework of linear theory. The incorporation of non-linear effects in the above systems to describe the evolution of the flow in the framework of Euler equations was introduced in [24].

The present work is in continuation with the above developments. It aims to provide a weakly nonlinear, simplified system, in the framework of water waves propagating in the presence of currents, the latter being supposed to be known a priori. There are various motivations for proceeding so. The first to be cited is numerical efficiency. Indeed, describing a complex flow, constituted by both the current and the wave-induced flow, might require more modes to reach enough accuracy. Additionally, in order to tackle the multiscale aspects involved in coastal regions, there is a growing interest in coupling numerical models, relying on separate tools to resolve the dynamics of the slow scales (the currents) and the faster ones (the water waves), which the present work focuses on. A physical argument should also be cited. Disentangling the influence of water waves on the mean flow and the modifications enforced to the perturbative field by the mean flow is a necessary step to better understand wave-current interactions.

The structure of the paper is as follows: The derivation of the nonlinear CMS, separating the mean flow from the perturbative field, is presented in Section 2. The reduction of this model to a single mode allows us to obtain a simplified, two-equation system, aiming to model the dynamics of weakly nonlinear water waves propagating over varying bathymetry, on top of arbitrarily sheared currents.

The forcing of this model, at the wavemaking boundary, requires the knowledge of periodic, traveling nonlinear water waves in the presence of vertically sheared currents over a flat bottom. For this specific purpose, appropriate solutions are derived here by means of the stream function formulation, following the semi-analytical method initially introduced by Rienecker and Fenton [25] and extended to the specific case of water waves propagating above linearly sheared currents by Francius and Kharif [26]. A further extension of this approach, to more general cases of currents, is described in Section 3.

Finally, the numerical results illustrating the performance of the nonlinear CMS integrated in the time domain are presented in Section 4, and finally, conclusions and directions for future research are discussed in Section 5.
## 2. The Weakly Nonlinear Coupled-Mode System

The present work is based on the derivation initially presented in Ref. [24]. In the following subsections, a description of the problem under consideration is presented, briefly reminding readers of the basic theory and explaining its extension.

### 2.1. Decomposition of the Perturbative Field in a Local-Mode Series

As illustrated in Figure 1, we consider for simplicity the two-dimensional problem of weakly nonlinear water waves propagating over currents which vary with depth according to a general vertical profile. Both the bathymetry characterized by the depth function $h(x)$ and the vertical structure of the steady current $U_0(x, z)$ are considered to be known, where bold notation is used to indicate vector quantites. The flow domain $D(t)$ is bounded below by an impermeable boundary defined the depth function $z = -h(x)$, and it is bounded above by the free surface, which is modeled by the function $z = \eta(x, t)$, denoting the free surface elevation. It is assumed that $\eta(x, t) + h(x) > 0$. The velocity field is considered to result from the contribution of steady rotational currents and a perturbative flow describing water waves:

$$U(x, z, t) = U_0(x, z) + u(x, z, t). \qquad (1)$$

<figure><img src="image_4.png" /><figcaption>Figure 1. Waves propagating in variable bathymetry in the presence of sheared currents.</figcaption></figure>

The problem to be solved only considers the perturbative field, $u(x, z, t)$. Following this purpose, this field can be described by a spectral-type representation defined as a series expansion, with the general expression given by the summation of local modes (indexed by an integer $n$) as follows:

$$u(x, z, t) = \sum_{n=0}^{+\infty} u_n(x, t) Z_n(z; h(x), \eta(x, t)),$$

where $Z_n(z; h(x), \eta(x, t))$ denote vertical functions which will be separately defined below for the horizontal and the vertical flow components. The latter vertical modes are different for each component of the wave velocity and are parametrically dependent on the local depth function and the free surface elevation, having the property of accurately treating the continuity equation and the no-entrance bottom boundary condition. In addition, the system of local vertical modes is complete in the interval $z \in [-h(x), \eta(x, t)]$ for each horizontal $x$ position.
More specifically, if the perturbative velocity field is noted $u(x,z,t) = (u(x,z,t), w(x,z,t))$, the horizontal wave component is given by

$$u(x, z, t) = \sum_{n=0}^{+\infty} u_n(x, t) Z_n^{(1)}(z; h(x), \eta(x, t)) \quad (2)$$

where $Z_n^{(1)}(z; h(x), \eta(x, t))$ are defined by

$$Z_n^{(1)}(z; h, \eta) = \frac{\cos[k_n(z+h)]}{\cos[k_n(\eta+h)]}, \quad n = 0, 1, 2, 3, \dots \quad (3)$$

and the quantities $k_n = k_n(h(x), \eta(x, t))$, $n = 0, 1, 2, \dots$, are obtained as the roots of the following (dispersion-like) equation:

$$\mu_0 + k_n \tan[k_n(h + \eta)] = 0 \quad (4)$$

The first root is imaginary ($k_0 = i|k_0|$), while the positive index roots $k_n$, $n = 1, 2, \dots$, are real, with $k_n \approx \frac{n\pi}{(h + \eta)}$, $n \to \infty$. Function (3), together with the dispersion Equation (4), are obtained from regular eigenvalue problems in the local vertical interval $z \in [-h(x), \eta(x, t)]$, controlled by the frequency-type parameter $\mu_0$; see [16]. This directly ensures the completeness property of the vertical basis; see, e.g., [27,28].

The vertical component $w(x, z, t)$ of the wave field is modeled by a corresponding modal expansion based on another set of z-functions $Z_n^{(2)}(z; h(x), \eta(x, t))$, as follows:

$$w(x, z, t) = - \sum_{n=0}^{+\infty} \left( \frac{\partial u_n}{\partial x}(x; t) Z_n^{(2)}(z; h; \eta) + u_n(x; t) \frac{\partial Z_n^{(2)}}{\partial x}(z; h; \eta) \right) \quad (5)$$

The functions $Z_n^{(2)}(z; h, \eta)$ are defined as follows:

$$Z_n^{(2)}(z; h, \eta) = \frac{\sin[k_n(z+h)]}{k_n \cos[k_n(\eta+h)]}, \quad n = 0, 1, 2, 3, \dots \quad (6)$$

At the seabed, they fulfill $Z_n^{(2)}(z = -h; h, \eta) = 0$, as well as the condition that $\partial Z_n^{(2)} / \partial z(z) = Z_n^{(1)}(z)$ over the whole water column (see also Appendix A of Ref. [24]). Using the latter relations, in conjunction with the definition of $u(x, z, t)$, given by Equations (2) and (3), and $w(x, z, t)$, given by Equations (5) and (6), in the continuity equation $u_x + w_z = 0$ and in the bottom boundary condition $w(x, z = -h, t) + \partial h / \partial x u(x, z = -h, t) = 0$, it is easily seen that both these equations are automatically satisfied by the present modal-type expansion of the wave velocity components. This property greatly facilitates the treatment of the considered problem, since only the momentum and the free surface boundary conditions are left to be satisfied. Furthermore, the substitution of the above local-mode expansions in the momentum equations will lead to a nonlinear coupled-mode system of equations on the horizontal plane with respect to the unknown velocity amplitudes $u_n(x, t)$ and the wave elevation $\eta(x, t)$. Based on the fast convergence properties of the local mode series, the above coupled-mode system facilitates the numerical solution of the problem and will be discussed in the sequel.

## 2.2. Momentum Equations

As discussed in the previous subsection, the kinematics of the problem are analytically satisfied by means of the specific modal expansions of the wave velocity. We might now address the dynamics of the problem, by introducing this representation in the momentum
equations. We shall start with the Euler equations describing the total flow $\mathbf{U} = (U, W)$, so that (see, e.g., [29])

$$ \frac{\partial U}{\partial t} + \frac{1}{2} \frac{\partial |\mathbf{U}|^2}{\partial x} - E = -\frac{1}{\rho} \frac{\partial p}{\partial x}, \\ \frac{\partial W}{\partial t} + \frac{1}{2} \frac{\partial |\mathbf{U}|^2}{\partial z} - e = -\frac{1}{\rho} \frac{\partial (p + \rho g z)}{\partial z}, \qquad (7) $$

where $|\mathbf{U}|^2 = U^2 + W^2$ and $e = U \times \Omega = (E, e)$ is a helicity-type field, with $\Omega$ denoting the vorticity $\Omega = \nabla \times \mathbf{U} = \Omega \hat{j}$, which, in accordance with Equation (1), is split into two terms: $\Omega = \Omega_0 + \omega$. In the present two-dimensional case, the vorticity is fixed in the transverse direction (with $\hat{j}$ denoting the corresponding unit vector), and thus,

$$ E = -W\Omega, \quad e = U\Omega \qquad (8) $$

Subsequently, at each horizontal x position, we may consider the vertical momentum Equation (7) and perform its integration in the vertical direction from a depth level z up to the free surface elevation $\eta(x, t)$. In conjunction with the dynamic free surface boundary condition, $p(x, z = \eta(x, t), t) = 0$, an expression of the local pressure in the column is obtained and reads as follows:

$$ \frac{p(x, z, t)}{\rho} = g(\eta(x, t) - z) - \int_{z}^{\eta(x, t)} e(z) dz + \int_{z}^{\eta(x, t)} \left( \frac{\partial W}{\partial t} + \frac{1}{2} \frac{\partial [U^2 + W^2]^{z=\eta(x,t)}}{\partial z} \right) dz. \qquad (9) $$

This expression of the local pressure might now be differentiated with respect to the horizontal coordinate, and it comes out to

$$ -\frac{1}{\rho} \frac{\partial p}{\partial x}(x, z, t) = -g \frac{\partial \eta}{\partial x}(x, t) + \frac{\partial}{\partial x} \left( \int_{z}^{\eta(x,t)} e(z) dz \right) - \frac{\partial}{\partial x} \left( \int_{z}^{\eta(x,t)} \frac{\partial W}{\partial t} dz \right) - \frac{1}{2} \frac{\partial}{\partial x} \left( \frac{\partial \left[ |\mathbf{U}|^2 \right]^{z=\eta(x,t)}}{\partial z} \right) $$

which, substituted back into Equation (7), finally results in

$$ \frac{\partial U}{\partial t} - E = -g \frac{\partial \eta}{\partial x}(x, t) + \frac{\partial F}{\partial x} - \frac{\partial}{\partial x} \left( \int_{z}^{\eta(x,t)} \frac{\partial W}{\partial t} dz \right) - \left[ U \frac{\partial U}{\partial x} + W \frac{\partial W}{\partial x} \right]^{z=\eta(x,t)}, \qquad (10) $$

for all $x$ and $h(x) < z < \eta(x, t)$, and where $F(x, z, t) = \int_{z}^{\eta(x,t)} e(z) dz$. At this point, we might introduce the flow decomposition (1) in Equation (10), to obtain the horizontal momentum equation satisfied by the only perturbative field. It follows that

$$ \frac{\partial u}{\partial t} - E_1 - E_2 - E_3 = -g \frac{\partial \eta}{\partial x} + \frac{\partial F_1}{\partial x} + \frac{\partial F_2}{\partial x} + \frac{\partial F_3}{\partial x} - \frac{\partial}{\partial x} \left( \int_{z}^{\eta(x,t)} \frac{\partial w}{\partial t} dz \right) - \left[ U_0 \frac{\partial u}{\partial x} + W_0 \frac{\partial w}{\partial x} \right]^{\eta(x,t)} - \left[ u \frac{\partial U_0}{\partial x} + w \frac{\partial W_0}{\partial x} \right]^{\eta(x,t)} - \left[ u \frac{\partial u}{\partial x} + w \frac{\partial w}{\partial x} \right]^{\eta(x,t)} \qquad (11) $$

where the quantities $E_1, E_2, E_3, F_1, F_2, F_3$ are associated with the expansion of the helicity field, and are defined as follows:

$$ \begin{aligned} E_1 &= -W_0\omega, \\ E_2 &= -w\Omega_0, \\ E_3 &= -w\omega, \\ F_1 &= \int_{z}^{\eta(x,t)} e_1(z) dz = \int_{z}^{\eta(x,t)} U_0 \omega dz, \\ F_2 &= \int_{z}^{\eta(x,t)} e_2(z) dz = \int_{z}^{\eta(x,t)} u \Omega_0 dz, \\ F_3 &= \int_{z}^{\eta(x,t)} e_3(z) dz = \int_{z}^{\eta(x,t)} u \omega dz \end{aligned} \qquad (12) $$
In Equation (11), the terms involving only the current have been dropped, since it is understood that the current flow is a solution of this equation by itself.

Finally, Equation (11), in conjunction with the kinematic free surface boundary condition,

$$ \frac{\partial \eta}{\partial t} + U_0 \frac{\partial \eta}{\partial x} + u \frac{\partial \eta}{\partial x} - w = 0, \quad \text{on } z = \eta(x, t), \qquad (13) $$

constitutes the remaining system to be solved.

## 2.3. A Weakly Nonlinear Coupled-Mode System

In this section, the remaining horizontal momentum Equation (11) of an inviscid liquid with a free surface and the kinematic free surface boundary condition, Equation (13), are re-formulated as a non-linear coupled-mode system with respect to the unknown amplitudes of the velocity field $(u_n(x, t))_{n=0,1,2,...}$ and the free surface elevation $\eta(x, t)$. Using the completeness properties of the local vertical functions $\{Z_m^{(1)}(z; h(x), \eta(x, t)), m = 0, 1, 2, ...\}$, both sides of Equation (11) are projected into this basis, for each horizontal $x$ position, and we obtain

$$ \left\langle \frac{\partial u}{\partial t}, Z_m^{(1)} \right\rangle + \left\langle \frac{\partial}{\partial x} \left( \int_z^\eta \frac{\partial w}{\partial t} d\zeta \right), Z_m^{(1)} \right\rangle \\ + \left\langle 1, Z_m^{(1)} \right\rangle \left[ \frac{\partial}{\partial x} [U_0(\eta) \cdot u(\eta) + W_0(\eta)w(\eta)] + \frac{1}{2} \frac{\partial}{\partial x} [u(\eta)^2 + w(\eta)^2] + g \frac{\partial \eta}{\partial x} \right] \\ - \left\langle E_1, Z_m^{(1)} \right\rangle - \left\langle E_2, Z_m^{(1)} \right\rangle - \left\langle E_3, Z_m^{(1)} \right\rangle - \left\langle \frac{\partial F_1}{\partial x}, Z_m^{(1)} \right\rangle - \left\langle \frac{\partial F_2}{\partial x}, Z_m^{(1)} \right\rangle \\ - \left\langle \frac{\partial F_3}{\partial x}, Z_m^{(1)} \right\rangle = 0, \quad \text{for } m = 0, 1, 2, \dots, M, \qquad (14) $$

where $\langle f_1(z), f_2(z) \rangle = \int_{z=-h}^{z=\eta} f_1(z) f_2(z) dz$. As already mentioned, this equation is to be supplemented with the kinematic free surface boundary condition, which, in terms of the wave velocity modal amplitudes, reads as follows:

$$ \frac{\partial \eta}{\partial t} + U_0 \frac{\partial \eta}{\partial x} + \sum_{n=0}^{M} \left\{ \left( Z_n^{(1)}(z=\eta; h, \eta) \frac{\partial \eta}{\partial x} + \frac{\partial Z_n^{(2)}}{\partial x}(z=\eta; h, \eta) \right) u_n \right. \\ \left. + Z_n^{(2)}(z=\eta; h, \eta) \frac{\partial u_n}{\partial x} \right\} = 0, \quad \text{on } z=\eta. \qquad (15) $$

The final form of the above coupled-mode system, Equations (14) and (15), is derived by using the particular forms of the vertical functions defined by Equations (3) and (6), as provided with more details in Appendix A, and the involved coefficients are given in Appendix B. Moreover, a weakly nonlinear system is obtained, as presented with more details in Appendix B, by treating the nonlinear terms as second-order terms of the unknowns $\eta$ and $(u_n(x, t))_{n=0,1,2,...}$, which is studied in present work for the propagation of waves interacting with sheared currents over general bathymetry.

Furthermore, the convergence of the vertical expansion discussed in [24] reveals that, in most cases, the contribution of the first mode $n=0$ is dominant. Thus, it is tempting to produce a two-equation system obtained by considering only this term. By approximating the vertical structure of the wave field, keeping only the mode $n=0$ in the wave velocity expansion Equation (2), and assuming relatively slow currents and dropping approximately
higher-order terms, a simplified version of the above system of Equations (13) and (14) reduces to the following:

$$ 
\begin{aligned}
& A1_{00} \frac{\partial u_0}{\partial t} + B1_{00} \frac{\partial^2 u_0}{\partial t \partial x} + C1_{00} \frac{\partial^3 u_0}{\partial t \partial x^2} + D1_{00} u_0 + E1_{00} \frac{\partial u_0}{\partial x} + F1_{00} \frac{\partial^2 u_0}{\partial x^2} \\
& \qquad + G1_{00} \frac{\partial^3 u_0}{\partial x^3} + H1_{00} \frac{\partial \eta}{\partial x} = NL_0^{(1)}, \text{ and} \\
& \qquad \frac{\partial \eta}{\partial t} + A2_{0} \frac{\partial \eta}{\partial x} + \sum_{n=0}^{M} \left\{ B2_{0} u_{0} + C2_{0} \frac{\partial u_{0}}{\partial x} \right\} = NL^{(2)}
\end{aligned}
\quad (16) $$

where $NL_0^{(1)} = \langle E_3, Z_0^{(1)} \rangle + \langle \frac{\partial F_3}{\partial x}, Z_0^{(1)} \rangle - \langle 1, Z_0^{(1)} \rangle \frac{1}{2} \frac{\partial}{\partial x} [u_0(0)^2]$, $NL_0^{(2)} = -u_0(0) \frac{\partial \eta}{\partial x}$, and the coefficients are defined in Appendix B. Here, the brackets denote $\langle f_1, f_2 \rangle = \int_{z=-h}^{z=0} f_1 f_2 dz$. The above simplified version of the weakly nonlinear CMS could be useful to obtain approximate numerical results at relatively low-cost supporting engineering applications concerning waves propagating in the presence of shear currents over inhomogeneous domains.

### 3. Nonlinear Solutions of Steady Traveling Waves and Sheared Currents in Constant Depth

Following the purpose of initializing our numerical solver of the coupled-mode system (either by means of the boundary condition or the initial condition in the restrictive case of an undisturbed boundary), we need to have access to the modal amplitude corresponding to a traveling wave over unperturbed bathymetry. This purpose can be achieved by solving the fully nonlinear problem and then computing its projection on the CMS modal base. The purpose of this section is to present the method used here.

Restricting ourselves to the case of periodic, two-dimensional, steady-traveling waves at a constant depth $h$, in the presence of a shear current $U_0(z)$, we consider a coordinate system $(\xi, z)$ moving with the wave speed $c$ so that $\xi = x - ct$, with the origin at some point on the undisturbed free surface. Following [26,30], a specific solution is constructed in this section, obtained as solution of the modified Helmholtz equation for the stream function $\Psi(\xi, z)$:

$$ \nabla^2 \Psi - a \Psi = b \quad (17) $$

with $a, b$ constants subjected to Dirichlet boundary conditions at the bottom and at the free surface:

$$ \Psi(\xi, z = -h) = 0, \quad \Psi(\xi, z = \eta(\xi)) = -Q \quad (18) $$

where $\eta(\xi)$ denotes the free surface elevation and $Q$ is a constant controlling the flowrate. The solution is constructed to also satisfy the dynamic free surface boundary condition as it will be described below, in conjunction with horizontal periodicity for a cell of length $\lambda$. Also, the above solution corresponds to a field where the vorticity is proportional to the stream function:

$$ \Omega = a\Psi + b \quad (19) $$

In the sequel, we will consider examples in the case where $a > 0$; however, extension to $a < 0$ is also possible. We consider the following partial solution of Equation (18):

$$ \Psi_0(z) = \frac{B_0}{a^{1/2}} \sinh(a^{1/2}(z+h)) + \frac{b}{a} \left( \cosh(a^{1/2}(z+h)) - 1 \right), \quad (20) $$

where $B_0$ is an unknown coefficient that will be connected with the phase speed of the waves $c = -B_0$. It is easily verified that the function provided by Equation (21) satisfies the partial differential Equation (18) and the bottom boundary condition. By differentiating with respect to $z$ we obtain that it corresponds to a horizontal flow with exponential structure in the vertical direction:

$$ U_0(z) = \frac{\partial \Psi_0}{\partial z} = B_0 \cosh(a^{1/2}(z+h)) + \frac{b}{a^{1/2}} \sinh(a^{1/2}(z+h)), \quad (21) $$
with vorticity described as follows:

$$ \Omega_0(z) = \frac{\partial U_0}{\partial z} = B_0 a^{\frac{1}{2}} \sinh\left(a^{\frac{1}{2}}(z+h)\right) + b \cosh\left(a^{\frac{1}{2}}(z+h)\right). \quad (22) $$

It can be further seen that in the case of small a, Equations (21)–(23) above reduce to the following approximate vertical profile:

$$ \Psi_0(z) \approx B_0(z+h) + \frac{b(z+h)^2}{2} + O(a) \quad (23) $$

$$ U_0(z) = \frac{\partial \Psi_0}{\partial z} \approx B_0 + b(z+h) + O(a) \quad (24) $$

$$ \Omega_0(z) = \frac{\partial U_0}{\partial z} \approx b + O(a). \quad (25) $$

Comparing the above with the result of Ref. [26] we obtain a result that is similar to the corresponding case of waves propagating in the presence of a linearly sheared current with constant vorticity $\Omega_0 \approx b$. On the other hand, for $b = 0$, the vertical current profile is described by the hyperbolic cosine function, and in cases of relatively large a values, it takes the form of an exponential distribution.

### 3.1. Solution Based on Fourier Modal Expansion Method

We proceed now to the formulation of the Fourier series method for the solution of the above problem, based on the following representation:

$$ \Psi(\xi, z) = \frac{B_0}{a^{1/2}} \sinh\left(a^{1/2}(z+h)\right) + \frac{b}{a}\left(\cosh\left(a^{1/2}(z+h)\right) - 1\right) + \sum_{n=1}^{\infty} B_n \Psi_n(\xi, z), \quad (26) $$

where $\Psi_n(\xi, z)$ are eigen solutions of the homogeneous modified Helmholtz equation, Equation (18) with $b = 0$, in the domain fulfilling the bottom boundary condition):

$$ \Psi_n(\xi, z) = \frac{\sinh(k_n(z+h))}{\cosh(k_n h)} \cos(nk\xi), \quad \text{with } k_n^2 = a + (nk)^2, \quad (27) $$

where $k = 2\pi/\lambda$ is the horizontal periodicity parameter. In the present case, where the vorticity is assumed to be proportional to the stream function (20), the Euler equations are integrated after being expressed in the moving frame of reference using $\tilde{u} = u - c$ and $\tilde{\xi} = x - ct$, providing Bernoulli's equation as follows:

$$ \frac{p}{\rho} - R - \frac{a}{2}\Psi^2 - b\Psi + gz + \frac{1}{2}(u^2 + w^2) = 0, \quad (28) $$

where R denotes the Bernoulli constant; see, e.g., [25,29]. In Equation (29), the horizontal and vertical flow velocity components are obtained from the corresponding derivatives of the stream function representation as follows:

$$ \tilde{u}(\xi, z) = \frac{\partial \Psi(\xi, z)}{\partial z}, \quad w(\xi, z) = -\frac{\partial \Psi(\xi, z)}{\partial \xi}. \quad (29) $$

Obviously, the case $a = b = 0$ results to the method of Rienecker and Fenton [25] corresponding to periodic waves without current. Using the above equations, the following system of nonlinear equations is defined:

$$ \frac{B_0}{a^{-1/2}} \sinh\left(a^{1/2}(\eta(\xi) + h)\right) + \frac{b}{a}\left(\cosh\left(a^{1/2}(\eta(\xi) + h)\right) - 1\right) + \sum_{n=1} B_n \Psi_n(\xi, z = \eta(\xi)) + Q = 0, \quad (30) $$
$$g\eta(\xi) + \frac{1}{2}(u^2(\xi, z = \eta(\xi)) + w^2(\xi, z = \eta(\xi))) - R - \frac{a}{2}Q^2 = 0, \quad (31)$$

From which the unknown coefficients {$B_n, n = 0, 1, 2, ...$}, the free surface elevation $\eta(\xi)$, and the unknown constants $Q$ and $R$ are calculated for given values of the parameters $a$ and $b$. For the numerical solution of the above system, the representation of the stream function is truncated, keeping the first $N$ terms {$B_n, n = 1, 2, ... N$}. The free surface elevation is discretized to $N$ equidistant points $\eta_i = \eta(\xi_i)$ with constant step $\Delta\xi = \xi_{i+1} - \xi_i$, and the discrete version of the system (31) and (32) with $2N + 2$ unknowns (including $R$ and $Q$) is numerically solved by the Newton-Raphson method, similarly as presented by [25]. The additional two equations are provided by using the data for the wave height $H$:

$$\max(\eta_i) - \min(\eta_i) = H, \quad (32)$$

and the requirement that the mean water level is at $z = 0$:

$$\text{Mean}(\eta_i) = 0. \quad (33)$$

## 3.2. Examples of Fields and Dispersion Relation Characteristics

In this subsection, indicative numerical solutions are presented as obtained by using the stream function method described above, for periodic waves characterized by the nonlinear parameter $\epsilon = 0.5Hk = 0.15$, where $H$ denotes the waveheight, and the dispersive parameter is chosen to be $kh = 2\pi/10$. The cases of linearly and exponentialy sheared currents are considered here. The selection of the linear case is motivated by the literature review, and the willingness to provide a comparison with already existing results. Water waves propagating in the presence of a vertically linear distribution of the current has been widely studied in previous years for mathematical reasons (see, e.g., Touboul et al. [22]). However, this configuration finds applications in geophysical flows, for instance, when considering rip currents, as demonstrated by Haas and Svendsen [2]. The selection of the exponential case is now allowed by the abilities newly obtained with this approach. It corresponds to the current distribution obtained under the action of wind. By comparing these two solutions, the strong influence of the current distribution is clearly illustrated. Here, the analysis is restricted to opposing currents; however, it could be extended to following currents and will be examined in future work. For normalizing the results, a water depth measured from the mean level of the free surface elevation equal to $h = 1$ m is used in the computations.

(i) *Linearly sheared current*: In this example, the current vorticity coefficients in Equation (21) are chosen, such as $a \approx 0 \text{ m}^{-2}$ and $b = -1 \text{ s}^{-1}$. The current profile and the resulting free surface elevation are presented in Figure 2. It is observed that the wavelength-to-depth ratio is $\lambda/h = 10$ and the wave nonlinearity parameter is $H/h = 0.46$, and thus, the Ursell number takes moderate values $Ur = (H/h)(\lambda/h)^2 = 46$, corresponding to nonlinear Stokes (short) waves.

The free surface presents a strong vertical asymmetry, as it is classically observed in nonlinear periodic wave solutions. The results are obtained by using a number of grid points ($N = 40$) corresponding to a total of 80 grid points per wavelength, which are found to be enough for numerical convergence. In Figure 3, color plots of the calculated stream function and the vorticity corresponding to the total field are presented. In addition, the streamlines are indicated in the subplots using thin black lines. It is observed that the total flow vorticity, in the major part of the domain, presents a value close to $\Omega = -1 \text{ s}^{-1}$, which is due to the constant value of the background current characterized by the linear vertical profile of Figure 2a. Deviations from the latter value are observed close to the free surface and especially near the wave crests as a result of the interaction of the current with the wave flow, which becomes rotational. Along with the solution concerning the free surface elevation and the stream function, the wave celerity is also calculated ($c = -B_0$). Systematic results concerning the dispersion characteristics of the above non-linear solutions will be presented in Section 3.3.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Vertical distribution of the linear background velocity profile. The parameters chosen are $a \approx 0$ and $b = -1 \text{ s}^{-1}$. (b) Free surface elevation obtained for a nonlinear parameter $\epsilon = 0.5Hk = 0.15$ and a dispersive parameter $kh = 2\pi/10$. The dotted line indicates the corresponding mean water level.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Spatial distribution of the total flow stream function, corresponding to the mean flow and the perturbative wave field. The current parameters chosen are $a \approx 0$ and $b = -1 \text{ s}^{-1}$, while the wave parameters are $\epsilon = 0.5Hk = 0.15$ and $kh = 2\pi/10$. (b) Spatial distribution of the total vorticity field in the same conditions. In both subfigures, thin black lines correspond to the streamlines.</figcaption></figure>
(ii) Exponentially sheared current: In this case, we consider an exponentially sheared current field. A particular example is calculated by using the parameters $a = 0.7 \text{ m}^{-2}$ and $b = 0 \text{ s}^{-1}$. The corresponding profile is presented in Figure 4. The free surface obtained numerically by applying the present stream function method is also plotted in this figure. The wave parameters are as before, with $\varepsilon = 0.5Hk = 0.15$ concerning the nonlinear parameter and $kh = 2\pi/10$ for the dispersive parameter, while the numerical configuration is kept identical to the previous case. The free surface elevation obtained with these parameters is found to be close to the one obtained in the presence of a linearly sheared current; this is as expected since the two vertical profiles of Figures 2a and 4a differ mostly at a submergence close to half of the depth, where the wave is relatively attenuated.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Vertical distribution of the background velocity profile. The parameters chosen are $a = 0.7 \text{ m}^{-2}$ and $b = 0 \text{ s}^{-1}$. (b) Free surface elevation obtained for a nonlinear parameter $\varepsilon = 0.5Hk = 0.15$ and a dispersive parameter $kh = 2\pi/10$. The dotted line indicates the corresponding mean water level.</figcaption></figure>

In this case, also, the vertical asymmetry of the wave profile is important, presenting very close values of the wave crests and through as before in Figure 2. This is further verified by the values of the constant controlling the flowrate, which is calculated to be $Q = 2.991 \text{ m}^2/\text{s}$, which is very close to the value in the linearly sheared current case ($Q = 2.998 \text{ m}^2/\text{s}$).

In Figure 5, color plots of the computed stream function and derived vorticity are presented. As before, the streamlines are indicated in both subplots using thin black lines. As far as it concerns the calculated stream function $\Psi$ (Figure 5a), a strong similarity with the linear sheared current is observed, although some fine differences might be noticed. Yet, when comparing the total vorticity field $\Omega$ between the linear sheared current (Figure 3b) and the exponential sheared current (Figure 5b), specific differences are observed. Indeed, it appears clearly in the example of the exponentially sheared current that the vorticity is much smaller (almost zero) in the vicinity of the bottom, which is a direct consequence of
the vertical distribution of the shear current $U_0(z)$ and its derivatives in the neighborhood of $z = -h$, which, for the relatively large $a$ value of this example, is approximated by an exponential distribution, resulting in $\omega_0(z) = \partial U_0/\partial z \simeq 0$ in the vicinity of the bottom. Here, again, the total vorticity including the waves is observed to present a deviation from the corresponding field associated with the current $\Omega_0(z)$ in the vicinity of the free surface and finds its maximum (in magnitude) at the surface.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Spatial distribution of the total stream function, corresponding to the mean flow and the perturbative field. The current parameters chosen are $a = 0.7 \text{ m}^{-2}$ and $b = 0 \text{ s}^{-1}$, while the wave parameters are $\epsilon = 0.5 \text{Hk} = 0.15$ and $kh = 2\pi/10$. (b) Spatial distribution of the total vorticity field in the same conditions. In both subfigures, thin black lines correspond to the streamlines.</figcaption></figure>

### 3.3. Projection of the Fully Nonlinear Solution on a CMS Basis

Having obtained the numerical solution for the fully nonlinear stream function, the components of the wave field are obtained by spatial differentiation of the field using Equation (30). The latter components, in conjunction with the transformation from moving to the laboratory-fixed coordinate system, can be projected accordingly on the vertical basis of Equation (3). Indeed, considering the nonlinear solution of the horizontal component of the velocity, $u$, obtained by the method described in Section 3.1, and selecting a number of modes $M$ for truncation of Equation (2), a classical Galerkin approach can be used. This is easily obtained by considering that

$$ \sum_{n=1}^{M} u_n(x,t) \int_{-h}^{\eta(x,t)} Z_m^{(1)} Z_n^{(1)} dz = \int_{-h}^{\eta(x,t)} u(x,z,t) Z_m^{(1)}(x,z;h,\eta) dz, \quad 1 \le m \le M, \quad (34) $$
which, in matrix form, reads as $u = M^{-1}U$, with the following elements:

$$M_{mn} = \int_{-h}^{\eta(x,t)} Z_m^{(1)} Z_n^{(1)} dz, \quad U_m = \int_{\{-h\}}^{\eta(x,t)} u(x,z,t) Z_m^{(1)}(x,z; h, \eta(x,t)) dz. \quad (35)$$

Here, the functions $Z_n^{(1)}(z; h, \eta)$ and $Z_m^{(1)}(z; h, \eta)$ are provided by Equation (3). The modal amplitudes $(u_n)_{1 \le n \le M}$ are immediately obtained by inverting the linear system of Equation (34). These modal amplitudes $u_n(x, t)$ together with the free surface elevation $\eta(x, t)$ constitute incident wave data for the CMS Equations (14) and (15) or the simplified system Equation (16). Yet, it is emphasized that these integrals differ from the scalar product defined earlier, since the limit of integration is now between $z = -h$ and $z = \eta(x, t)$. The system in (34) is inverted to obtain the expected distribution of $u_n(x, t = 0)$. Whenever needed, boundary data $u_n(x = 0, t)$ at the entrance ($x = 0$) of the CMS domain are obtained by multiplying the above modal amplitude functions by the nonlinear phase velocity. Examples of applications of the latter for various cases will be presented in Section 4.

Yet, it is interesting to consider the converging behavior of the solution represented by the modal expansion of the CMS with the fully nonlinear solution of Section 3.2. In Figure 6, the vertical distribution of the horizontal component of the velocity, obtained in our first example of linearly sheared current ($a \simeq 0, b = -1s^{-1}, \varepsilon = 0.15$, and $kh = 2\pi/10$), is presented. The blue crosses correspond to the fully nonlinear solution, obtained at various phase values. The green curves correspond to the solution obtained for the only mode $n = 0$, at the same phase values. The orange curves present the solution involving the two first modes, $n = 0$ and $n = 1$, while the red curves are the solution obtained with 100 modes. These results are presented in order to provide an idea of the rate of convergence of the solution presented on the basis of the CMS. If only the mode $n_0$ shows some limitations (yet, the reader is reminded that the vorticity and the nonlinearity considered here are really significant), taking into account two modes already provides a very satisfying solution. A few differences can still be observed, located near the free surface. The example involving $M = 100$ modes shows absolutely no difference with the fully nonlinear solution. Another possibility is to consider the modal amplitude and its evolution with the number of modes considered. This is presented in Figure 7. From this figure, it is clear that the modal amplitude provides a rate of convergence which is larger than an order of two.

<figure><img src="image_9.png" /><figcaption>Figure 6. Vertical distribution of the horizontal component of the velocity, for selected phase values. Waves and current conditions correspond to the linearly sheared current case ($a \simeq 0, b = -1s^{-1}, \varepsilon = 0.15$, and $kh = 2\pi/10$).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Normalized amplitude of the modes involved in the CMS decomposition, as a function of the number of modes. Waves and current conditions correspond to the linearly sheared current case ($a \approx 0, b = -1 \text{ s}^{-1}, \varepsilon = 0.15$, and $kh = 2\pi/10$).</figcaption></figure>

### 3.4. Phase Velocity of the Nonlinear Waves

The method presented in this section also allows us to compute the nonlinear phase velocity, since it is related to the constant $B_0 = -c$. Furthermore, the linear phase velocity can be obtained in the case of arbitrary current flow by solving the Rayleigh equation (see [31] for more details).

Thus, it is interesting to compute the phase velocity obtained here for the specific current profiles, Equation (21), with the phase velocity of linearized waves computed by solving the Rayleigh equation. Results are presented in Figure 7, from systematic solutions as the ones examined before in the water strip with depth $h = 1$ m.

In particular, three current configurations are examined and results are presented: (i) water waves propagating in the absence of current ($a, b) = (0,0$), (ii) the current configuration presenting a linear variation in depth ($a, b) \approx (0,-1 \text{ s}^{-1}$), and (iii) the current varying exponentially with depth ($a, b) = (0.7 \text{ m}^{-2}, 0$). For each of these current configurations, three values of the dispersive parameter are considered: $kh = \pi$, $kh = 2\pi/5$, and $kh = \pi/5$, corresponding to water waves propagating in relatively deep water, intermediate water depth, and water waves propagating in relatively shallow water. These configurations constitute the nine curves presented in Figure 8, where the calculated normalized phase speed $c/(gh)^{1/2}$ is plotted. Each curve is obtained by varying the nonlinear parameter $\varepsilon$ in the vertical axis in the range $0.001 \le \varepsilon \le 0.2$.

The above results illustrate the effect of the nonlinear parameter on the phase celerity of traveling waves in the presence of sheared currents. The vertical black lines represent the values obtained for each of the cases considered by solving the Rayleigh equation, which match the calculated results for limiting small nonlinearity $\varepsilon \approx 0$.

Several classical results are confirmed through this figure. First, when considering a given current configuration, the effect of the dispersive parameter $kh$ is observed. The longer the waves, the faster they propagate. Also, for a given current configuration, it is seen that the effect of the nonlinear parameter $\varepsilon$ is more important for more shallow waves. The deviation from the linear theory (the vertical black lines) is faster when $kh$ is small. Finally, when intercomparing the effect of the current on the phase celerity, for a given $kh$, the opposing current is found to slow down propagation. The linearly sheared current presents a stronger effect than the exponentially sheared one, as expected.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Nonlinear phase velocity obtained with the numerical method presented in Section 3, for $kh = \pi$, $kh = 2\pi/5$, and $kh = \pi/5$, for various nonlinear parameters $0.001 \le \varepsilon \le 0.2$. Three current configurations are considered: $(a, b) \simeq (0, 0)$, corresponding to the solution in the absence of current; $(a, b) \simeq (0, -1 \text{ s}^{-1})$, which is the linearly sheared current; and $(a, b) = (0.7, 0)$, which corresponds to the exponentially sheared current.</figcaption></figure>

# 4. Weakly Nonlinear Propagation over Varying Depth and Currents

In this section, numerical results are presented in the case of waves propagating in constant depth in the presence of shear currents, as obtained by the time integration of the weakly nonlinear coupled-mode system (Equations (14) and (15)). Taking into account the observed fast convergence characteristics of the modal expansion, the simplified model based on Equation (16) is finally used to illustrate results of simulations of waves and shear currents in variable-bathymetry regions.

## 4.1. Initialization of the CMS and Comparison with the Fully Nonlinear Method

The purpose of this section is to investigate the ability of the CMS to reproduce the behavior of water waves in the presence of vortical flow in weakly nonlinear conditions. To stay consistent with this aim, three cases are considered, involving weakly nonlinear water waves and various current configurations. First, a no-current case in which $(a, b) \simeq (0, 0)$ is presented. The case of linear shear current in which $(a, b) \simeq (0, -0.5 \text{ s}^{-1})$ is then analyzed. Finally, an exponential case is presented, in which $(a, b) = (0.15 \text{ m}^{-2}, 0)$. For these three cases, the water waves parameters are $kh = 2\pi/10$ and $\varepsilon = 0.05$, characterizing weakly non-linear conditions in intermediate water depth. Results will be compared to those provided by the fully nonlinear method presented in Section 3.

The configuration of the coupled-mode system considered here involves periodical boundary conditions, while time integration of an initial condition is performed. To produce the initial condition, two methods were investigated. The first approach was to use a linear solution and to apply a relaxation term for the nonlinear terms in the first period of propagation (see, e.g., [22,23]). Another possibility is to start from the fully nonlinear solution obtained from the model presented in Section 3, providing the full solution in terms of $\eta$ and $u$ [32]. In this case, it is necessary to obtain its projection on the solution on the basis of the CMS, as described in Section 3.3. For the sake of clarity, it is reminded that the system provided by Equation (34) is solved for each location $x$ of the initial domain for the CMS, providing the data for $u_n(x, t = 0)$.

The discrete version of the present CMS is obtained by using high-order finite differences to approximate horizontal and time derivatives of the wave mode amplitudes $u_n(x, t)$ and the free surface elevation $\eta(x, t)$, on a regular grid ($x_i = x_a + (i - 1)\Delta x$, $t_n = (n - 1)\Delta t$).
The mesh is obtained by using $M_x$ steps of equal length $\Delta x = (x_b - x_a)/M_x$ to subdivide the horizontal interval $x_a \le x \le x_b$. Based on the above analysis concerning the dispersion characteristics of the present model, in most cases, the number of the retained modes required for convergence is very small, and the horizontal subdivision of the domain is based on using 30–40 points per wavelength. In all examples presented and discussed below, the dimension of the discrete coupled-mode system is of the order of several thousands, and for the time integration of the dynamic system, a fourth-order BDF scheme is used, with Courant number $(C\Delta t/\Delta x)_{max}$.

Both initializing procedures are compared in the case involving no current, and results are presented in Figure 9, concerning $\eta(x,t)$ and $u_n(x,t)$, $n = 0,1,2$. In this figure, results indicated by using red lines corresponding to the calculations based on the fully nonlinear method presented in Section 3, obtained by means of the stream function formulation. Green lines present the time evolution computed by means of the CMS using a linear initial condition, together with a relaxation scheme for the nonlinear terms. Finally, blue lines correspond to results obtained with the CMS, initialized with the projection described by Equation (34). It appears that the results of computation are compatible with the steady traveling solutions, by using both initialization approaches, and in the next examples, only the projection method (Section 3.3) is used.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Time evolution of the free surface elevation and the first, second, and third mode of the horizontal velocity. The parameters considered are $(a, b) = (0,0)$, $\epsilon = 0.05$, and $kh = \pi/5$, in the case of no current. Red lines correspond to the results obtained with the fully nonlinear method presented in Section 3. Green presents the results obtained with the CMS, initialized with the relaxation method described in Section 4.1. Blue lines correspond to results obtained with the CMS, initialized with the projection described in Section 3.3.</figcaption></figure>

Next, in Figure 10, the time evolution of the free surface elevation and the first, second, and third mode of the horizontal velocity are shown, as obtained by the present CMS, in the case of an opposing linear shear current characterized by the parameters $(a, b) = (0, -0.5 \text{ s}^{-1})$, $\epsilon = 0.05$, and $kh = \pi/5$.

Here again, the results presented are in good agreement with the fully nonlinear simulations. It is observed that no differences can be seen as far as the free surface elevation and the propagating mode $u_0$ are concerned. When it comes to modal amplitudes $u_1$ and $u_2$, and noting that the vertical scales of the last two subplots of Figure 10 have been increased by 25 and 50 times for comparison purposes, the results obtained by the present method are considered to be very satisfying, both in terms of phase and amplitude. Yet, some small oscillations in the modal amplitude are observed, related to the numerical integration. However, the latter differences are minor concerning the solution for the total velocity field
$u(x,z,t)$, which, in the examined nonlinear case, is approximated by truncating the present expansions defined by Equations (2) and (5) and keeping only the first three terms.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Time evolution of the free surface elevation and the first, second, and third mode of the horizontal velocity. The parameters considered are $(a, b) = (0, -0.5 \text{ s}^{-1})$, corresponding to a linear shear current. Moreover, $\varepsilon = 0.05$ and $kh = \frac{\pi}{5}$. Red lines correspond to the results obtained with the fully nonlinear method presented in Section 3. Blue lines correspond to results obtained with the CMS, initialized with the projection described in Section 4.1.</figcaption></figure>

Finally, in Figure 11, similar results are presented in the case of an opposing current showing an exponential vertical profile, characterized by the parameters $(a, b) = (0.15 \text{ m}^{-2}, 0)$, in the case of nonlinear waves in intermediate water depth $\varepsilon = 0.05$ and $kh = \pi/5$. The time evolution of the free surface elevation and the first, second, and third mode of the horizontal velocity are shown, as obtained by the present CMS, and compared to the fully nonlinear solution of Section 3.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Time evolution of the free surface elevation and the first, second, and third mode of the horizontal velocity. The parameters considered are $(a, b) = (0.15 \text{ m}^{-2}, 0)$ corresponding to the case of exponential current. Red lines (-) correspond to the results obtained with the fully nonlinear method presented in Section 3. Blue lines (-) correspond to results obtained with the CMS, initialized with the projection described in Section 3.3.</figcaption></figure>
In this example, the results presented in Figure 11 indicate a small discrepancy between the stream function solution and the CMS calculation, which is due to a very small difference (of the order of 1.25%) in wave celerity as calculated by the present CMS. Yet, when it comes to the behavior of the amplitudes, for the free surface elevation, but also for the modes $u_0$, $u_1$, and $u_2$, the agreement is still very satisfying.

In both the above examples, we also observe the fast decay of the modal amplitudes (accordingly to Figure 7), indicating that in such cases, $u_n(x,t)$, for $n > 0$, is more than an order of magnitude less than $u_0(x,t)$, confirming the applicability and usefulness of the simplified model provided by Equation (16). Still, some limitations of the numerical integrator of the system are observed, and future work will be devoted to further improving its performance.

## 4.2. Wave Propagation over Rippled Bed in the Presence of Shear Current

Finally, we consider here the application of the present CMS for an example concerning wave propagation over variable bathymetry in the presence of an opposing sheared current. The problem of water waves reflected by periodic undulating topography, known as Bragg resonant reflection, has been studied extensively in the last few decades. When the surface wavelength is about twice the topographic wavelength, the incoming waves are significantly reflected by the periodic bottom, and only a small part of the wave energy is transmitted. Since the problem has strong implications in coastal engineering, it has been widely studied; see, e.g., the study by Mei [33] and the references therein. More recently, the nonlinearity involved in the process has attracted some of the attention of the scientific community [34,35]. In particular, an investigation of the effects of Bragg reflection on harbor oscillations has been presented by Gao et al. [36,37]. The nonlinear problem of wave-bottom interactions involving a piecewise linear depth varying current was investigated recently by Raj and Guha [38]. For this reason, the Bragg reflection problem constitutes an interesting configuration to explore through the CMS.

Here, we consider a sinusoidal bottom topography in the middle of the domain of an otherwise flat horizontal bottom of depth $h = 0.22$ m, i.e., the bottom slope and curvature are zero for a significant length at the ends of the domain. Such a configuration has been tested experimentally in the flume of SeaTech, University of Toulon, France, depicted in Figure 12, and measurements of the reflection coefficient were collected in order to study the Bragg resonance of wave and current systems and its effects on the reflection coefficient. An inclined perforated screen was used to generate the shear of the opposing current at the downwave end of the flume. The position of the wave probes for estimating the reflection and transmission of waves from measured wave data are shown in Figure 13. At the upwave end of the channel, an electromagnetic piston was used to generate regular waves. At the downwave end of the flume, a sloping beach was used to absorb waves. The opposing current was injected into the channel by a hydraulic pump, and a perforated screen was used to control the shear; see Figure 13. The sinusoidal bottom profile consists of 10 ripples of amplitude 0.035 m and length 0.5 m. The opposing shear current vertical profile is approximately linear:

$$U(x,z) = U_0(x) + S(x)z \qquad (36)$$

where $U_0(x)$ is the surface current and $S(x)$ denotes the shear. In this case, the current vorticity field is $\Omega_0(x) = S(x)$.

We consider as an example the reflection of waves over a rippled bed with the presence of an opposing current with a mean speed of 0.17 m/s. More details are provided in Laffite et al. [39], where estimations concerning the spatial variation of shear $S(x)$ are available. In this case, the shear is estimated to be negative near the entrance of the wave and becomes oscillatory above the sinusoidal bottom. Based on the above, the horizontal distribution
of the shear in the domain is calculated by assuming the conservation of the mass of the current flow as follows:

$$S(x) = \frac{2}{h(x)} \left( U_0(x) - \frac{h(0)}{h(x)} U_0(0) \right) + \frac{2h(0)}{h(x)^2} S(0) \quad (37)$$

where $x = a$ and $x = b$ denote the coordinates of the wave entrance and exit in the domain.

<figure><img src="image_15.png" /><figcaption>Figure 12. Periscopic view of the wave-current flume of SeaTech, University of Toulon, showing the sinusoidal bottom.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Sketch of wave flume illustrating the position of wave probes for estimating reflection and transmission of waves from measured data. Also, the screen used at the downwave end of the tank to form the opposing shear current is shown.</figcaption></figure>

Indicative distributions of current parameters used in the sequel for calculations of waves by the present CMS are presented in Figure 14.

<figure><img src="image_17.png" /><figcaption>Figure 14. Current data for waves propagating over the sinusoidal bottom of Figure 4. The surface current $U_0(x)$ in the case of no shear is plotted by using a blue line and the effect of shear $S(x)$ is shown by using a red line, respectively.</figcaption></figure>

As concerns the conditions for Bragg resonance, it is known that the resonant reflection of surface water waves by sinusoidal bars in the presence of a current is manifested when

$$k_p + k_m = K_b, \quad (38)$$
where $K_b$ is the wavenumber of the bottom profile, and $k_p$ and $k_m$ are obtained as the roots of the dispersion relation for opposing and following current, respectively.

$$ \sigma^2(x) = \sigma_0(x)\sigma_2(x) = k_p(x)g\tanh(k_p(x)h(x)), \quad (39) $$

where g is the gravity acceleration. Equation (39) is formulated with respect to an average relative frequency estimated by

$$ \sigma = \sqrt{\sigma_0\sigma_2} = \sqrt{(\omega \pm U_0 \cdot k_{p,m})(\omega \pm U_2 \cdot k_{p,m})}, \quad (40) $$

where $U_0$ is the surface current and $U_2 = U_0 - S(2d)$ is the value of the current at a penetration depth defined as $2d$; see also Touboul et al. [22]. An extended version of the above relation has been derived and discussed in Belibassakis et al. ([23], Section 2.4).

In the experimental tests, the incident wave frequencies ranged from 0.7 Hz to 1.5 Hz with amplitudes between 0.5 – 2 cm, ensuring relatively small steepness, both without a current and with the vertically sheared current. Predictions obtained from the present model are presented in Figures 15 and 16 in the case of sinusoidal incident waves with a frequency of 1.1 Hz, close to wave resonance conditions, and a wave amplitude of 1 cm, without current and with the consideration of the sheared current. In these figures, a space-time plot of the free surface elevation $\eta(x, t)$ is shown, as calculated by the present method. The result corresponds to periodic waves propagating over a sinusoidal bed (also depicted in the figure at an average depth of 22 cm), without any current (Figure 15) and with the effects of a sheared current with a linear vertical profile (Figure 16), as modeled, respectively, by using the data in Figure 14. The bottom profile varies only along the x-axis, and the opposing current is also directed along the x-axis, as indicated in the figures. Since the wave generated at the upwave end of the flume is harmonic, the plot of the free surface history at $t = 0$ (shown along the t-axis) is sinusoidal.

In both cases, a rapid decay of the mode amplitudes is observed. Moreover, we observe in Figures 15 and 16, illustrating the history of the free surface elevation over the sinusoidal bed, that the results indicate strong reflection. In the case of waves opposing the shear current, we clearly observe the shortening of the wavelengths and the steepening of waves in Figure 16, especially over the sinusoidal bottom. Moreover, in this case, the disturbance of the flow due to the combined effect of the rippled bed and the vorticity contained in the vertically sheared current is seen to generate a more complicated pattern, related to the wavenumber asymmetry due to the presence of the current.

<figure><img src="image_18.png" /><figcaption>Figure 15. Numerical results obtained by the present CMS for harmonic incident waves of frequency 1.1 Hz and height $H = 0.02$ m propagating over sinusoidal bed of mean depth $h = 0.22$ m without current. Spatio-temporal evolution of the calculated free surface elevation.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Same as in Figure 15, but in the presence of opposing shear current with data from Equations (37) and (38). The direction of the sheared current is indicated using arrows.</figcaption></figure>

In Figure 17, the reflection coefficient, plotted as a function of the wave frequency (in the laboratory frame of reference), is presented for both cases, considering the current or not. Together with the experimental results, we present results obtained numerically by means of the present CMS model. These results were obtained by considering various incident wave conditions, namely *a* = 0.002 m, *a* = 0.005 m, and *a* = 0.01 m. Thus, the effect of nonlinearity in the reflection process is observed. It appears that both the peak frequency and its amplitude are better described when the considered wave amplitude approaches the one considered in the experiments. This result confirms that although Bragg reflection is a linear process, nonlinear hydrodynamical interactions play a role in the interaction process. Further investigation of the present CMS, including results concerning the calculation of reflection and transmission coefficients for various wave frequencies and comparisons with results from other works and experimental data, will be considered in future work, as well as an extension to three-dimensional problems; see, e.g., [31,40].

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Reflection coefficient obtained experimentally with and without current, compared to the results obtained with the present CMS model.</figcaption></figure>
## 5. Conclusions

A coupled-mode system based on the velocity formulation of the wave field was developed and used for modeling wave scattering by inhomogeneous, sheared currents in variable bathymetry regions. The model is based on a modal series expansion of the wave velocity based on vertical eigenfunctions, which are dependent on local depth and flow parameters, including the propagating and evanescent modes. A new derivation of a simplified weakly nonlinear system was introduced using decomposition to a mean flow and a perturbative wave field, and it was subsequently used to study wave-current-seabed interaction problems in inhomogeneous domains. For defining data describing the incoming waves at the inlet boundary, knowledge of periodic traveling nonlinear water waves over a flat bottom is required. Specific solutions were derived using the semi-analytical method based on the stream function formulation for cases of water waves propagating above sheared currents with linear and exponential profiles. The obtained results clearly illustrate the effect of the current vortical distribution on the nonlinear behavior and phase speed of water waves, in the same order of magnitude as the dispersion parameter, suggesting that sheared currents are expected to have a strong influence on refraction and diffraction phenomena in variable-bathymetry regions. Examples concerning the propagation of waves and currents in cases with depth inhomogeneities were subsequently presented and discussed, illustrating the applicability and performance of the present model. In particular, Bragg reflection by sinusoidal bottom profiles was analyzed, and 5% to 10% differences in the reflection peak were observed in the presence of a current, due to the combined effects of nonlinearity and vorticity on the behavior of water waves. Future work is planned towards the extension of the model to three-dimensional problems, including the case of waves interacting with currents at an angle.

**Author Contributions:** Conceptualization, J.T. and K.B.; methodology, J.T. and K.B.; software, J.T.; writing—original draft preparation, J.T. and V.M.-M.; writing—review and editing, J.T. and K.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work presented here was funded by the French ANR under grant number ANR-21-ASM1-0003. The paper fees have been waived by J. Mar. Sci. Eng.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors would like to thank the French DGA, which funded this research through the ANR grant ANR-21-ASM1-0003. Additionally, J.T would like to thank the French DGA for supporting his visit to Athens. K.B. would like to thank the UTLN for supporting his visit to Seatech University Toulon.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Expansion of the Momentum Equations on the Basis

Starting from Equations (14) and (15), we now introduce the expansions given by Equations (2) and (5). Proceeding term by term, we obtain

$$ 
\begin{aligned}
\frac{\partial u}{\partial t} &= \sum_{n=0}^{\infty} Z_n^{(1)}(z; h, \eta) \frac{\partial u_n(x, t)}{\partial t}, \\
\frac{\partial}{\partial x} \left( \int_z^\eta \frac{\partial w}{\partial t} d\xi \right) &= -\sum_{n=0}^{\infty} \frac{\partial}{\partial x} \left[ Z_n^{(3)} \left( \frac{\partial u_n}{\partial t} \right) - Z_n^{(2)}(\eta) \frac{\partial \eta}{\partial x} \frac{\partial u_n}{\partial t} \right], \\
\frac{\partial}{\partial x} [U_0(\eta) \cdot u(\eta)] &= \sum_{n=0}^{\infty} \frac{\partial}{\partial x} (U_0(\eta) \cdot u_n), \\
\frac{\partial}{\partial x} [W_0(\eta) \cdot w(\eta)] &= -\sum_{n=0}^{\infty} \frac{\partial}{\partial x} \left[ W_0 \frac{\partial}{\partial x} \left( Z_n^{(2)}(\eta) u_n \right) \right].
\end{aligned}
$$
where $Z_n^{(3)}$ is defined by

$$Z_n^{(3)} = \int_z^\eta Z_n^{(2)}(\xi; h, \eta) d\xi.$$

The expressions of $E_1$, $E_2$, $F_1$, and $F_2$, involving the vorticity fields, as derived in [27], are as follows:

$$E_1 = - \sum_{n=0}^{+\infty} W_0 \left(Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2}\right) u_n + \left(2W_0 \frac{\partial Z_n^{(2)}}{\partial x}\right) \frac{\partial u_n}{\partial x} + \left(W_0 Z_n^{(2)}\right) \frac{\partial^2 u_n}{\partial x^2},$$

$$E_2 = \sum_{n=0}^{+\infty} \left( \Omega_{02} \frac{\partial Z_n^{(2)}}{\partial x} \right) u_n + \left( \Omega_{02} Z_n^{(2)} \right) \frac{\partial u_n}{\partial x},$$

$$F_1 = \sum_{n=0}^{+\infty} \left\{ \left[ \int_z^\eta \left( U_0 \left(Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2}\right) \right) d\xi \right] u_n + \left[ \int_z^\eta \left( 2U_0 \frac{\partial Z_n^{(2)}}{\partial x} \right) d\xi \right] \frac{\partial u_n}{\partial x} \right. \\ \left. + \left[ \int_z^\eta \left( U_0 Z_n^{(2)} \right) d\xi \right] \frac{\partial^2 u_n}{\partial x^2} \right\},$$

$$F_2 = \sum_{n=0}^{+\infty} \left[ \int_z^\eta \Omega_{02} Z_n^{(1)} d\xi \right] u_n,$$

where $\Omega_{02}$ is the background vorticity, defined as $\partial U_0 / \partial z - \partial W_0 / \partial x$, and where $Z_n^{(0)}$ is defined by $Z_n^{(0)} = \partial Z_n^{(1)} / \partial z$. For explicit nonlinearities, we obtain

$$\frac{1}{2} \frac{\partial}{\partial x} [u(\eta)^2] = \frac{1}{2} \sum_{n=0}^{\infty} \sum_{l=0}^{\infty} \frac{\partial}{\partial x} (u_n u_l)$$

$$\frac{1}{2} \frac{\partial}{\partial x} [w(\eta)^2] = \frac{1}{2} \sum_{n=0}^{\infty} \sum_{l=0}^{\infty} \frac{\partial}{\partial x} \left\{ Z_n^{(2)}(\eta) Z_l^{(2)}(\eta) \frac{\partial u_n}{\partial x} \frac{\partial u_l}{\partial x} + \frac{\partial Z_n^{(2)}}{\partial x}(\eta) \frac{\partial Z_l^{(2)}}{\partial x}(\eta) u_n u_l \right. \\ \left. + \frac{\partial Z_n^{(2)}}{\partial x}(\eta) Z_l^{(2)}(\eta) u_n \frac{\partial u_l}{\partial x} + Z_n^{(2)}(\eta) \frac{\partial Z_l^{(2)}}{\partial x}(\eta) \frac{\partial u_n}{\partial x} u_l \right\},$$

and

$$E_3 = \sum_{n=0}^{+\infty} \sum_{l=0}^{+\infty} \left\{ \frac{\partial Z_l^{(2)}}{\partial x} \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) u_n u_l + \left( 2 \frac{\partial Z_l^{(2)}}{\partial x} \frac{\partial Z_n^{(2)}}{\partial x} \right) \frac{\partial u_n}{\partial x} u_l \right. \\ \left. + \left( \frac{\partial Z_l^{(2)}}{\partial x} Z_n^{(2)} \right) \frac{\partial^2 u_n}{\partial x^2} u_l + Z_l^{(2)} \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) u_n \frac{\partial u_l}{\partial x} \right. \\ \left. + \left( 2 Z_l^{(2)} \frac{\partial Z_n^{(2)}}{\partial x} \right) \frac{\partial u_n}{\partial x} \frac{\partial u_l}{\partial x} + \left( Z_l^{(2)} Z_n^{(2)} \right) \frac{\partial^2 u_n}{\partial x^2} \frac{\partial u_l}{\partial x} \right\},$$

$$F_3 = \sum_{n=0}^{+\infty} \sum_{l=0}^{+\infty} \int_z^\eta \left( Z_l^{(1)} \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) u_n u_l + \left( 2Z_l^{(1)} \frac{\partial Z_n^{(2)}}{\partial x} \right) \frac{\partial u_n}{\partial x} u_l \right. \\ \left. + \left( Z_l^{(1)} Z_n^{(2)} \right) \frac{\partial^2 u_n}{\partial x^2} u_l \right) d\xi.$$

## Appendix B. Detailed Expression of the Coefficients

Considering the developments provided in Appendix A, the present coupled-mode system reduces to the following form:

$$\sum_{n=0}^{M} \left\{ A1_{mn} \frac{\partial u_n}{\partial t} + B1_{mn} \frac{\partial^2 u_n}{\partial t \partial x} + C1_{mn} \frac{\partial^3 u_n}{\partial t \partial x^2} + D1_{mn} u_n + E1_{mn} \frac{\partial u_n}{\partial x} + F1_{mn} \frac{\partial^2 u_n}{\partial x^2} + G1_{mn} \frac{\partial^3 u_n}{\partial x^3} \right\} \\ + H1_m \frac{\partial \eta}{\partial x} = NL_m^{(1)}, m = 0, 1, \dots, M, \text{ and}$$

$$\frac{\partial \eta}{\partial t} + A2_n \frac{\partial \eta}{\partial x} + \sum_{n=0}^{M} \left\{ B2_n u_n + C2_n \frac{\partial u_n}{\partial x} \right\} = NL^{(2)}$$
In the above system of equations, the right-hand side contains all nonlinear terms, and the coefficients are defined by

$$
\begin{align*}
A1_{mn} &= \left\langle \left( Z_n^{(1)} - \frac{\partial^2 Z_n^{(3)}}{\partial x^2} \right), Z_m^{(1)} \right\rangle \\
B1_{mn} &= - \left\langle 2 \frac{\partial Z_n^{(3)}}{\partial x}, Z_m^{(1)} \right\rangle \\
C1_{mn} &= - \left\langle Z_n^{(3)}, Z_m^{(1)} \right\rangle, \\
D1_{mn} &= \left\langle 1, Z_m^{(1)} \right\rangle \frac{\partial}{\partial x} \left( U_0(\eta) Z_n^{(1)}(\eta) - W_0(\eta) \frac{\partial Z_n^{(2)}}{\partial x}(\eta) \right) + \left\langle \left( W_0 \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) \right), Z_m^{(1)} \right\rangle \\
&\quad - \left\langle \left( \Omega_{02} \frac{\partial Z_n^{(2)}}{\partial x} \right), Z_m^{(1)} \right\rangle - \left\langle \frac{\partial}{\partial x} \left[ \int_z^\eta \left( U_0 \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) \right) d\xi \right], Z_m^{(1)} \right\rangle \\
&\quad - \left\langle \frac{\partial}{\partial x} \left[ \int_z^\eta \left( \Omega_{02} Z_n^{(1)} d\xi \right), Z_m^{(1)} \right\rangle, \\
E1_{mn} &= \left\langle 1, Z_m^{(1)} \right\rangle \left( U_0(\eta) Z_n^{(1)}(\eta) - W_0(\eta) \frac{\partial Z_n^{(2)}}{\partial x}(\eta) - \frac{\partial}{\partial x} \left( W_0(\eta) Z_n^{(2)}(\eta) \right) \right) \\
&\quad + \left\langle \left( 2W_0 \frac{\partial Z_n^{(2)}}{\partial x} \right), Z_m^{(1)} \right\rangle - \left\langle \left( \Omega_{02} Z_n^{(2)} \right), Z_m^{(1)} \right\rangle - \left\langle \left( \int_z^\eta \left( \Omega_{02} Z_n^{(1)} d\xi \right), Z_m^{(1)} \right\rangle \right. \\
&\quad \left. - \left\langle \left( \int_z^\eta \left( U_0 \left( Z_n^{(0)} + \frac{\partial^2 Z_n^{(2)}}{\partial x^2} \right) \right) d\xi + \frac{\partial}{\partial x} \left[ \int_z^\eta \left( 2U_0 \frac{\partial Z_n^{(2)}}{\partial x} \right) d\xi \right] \right), Z_m^{(1)} \right\rangle, \\
F1_{mn} &= - \left\langle 1, Z_m^{(1)} \right\rangle \left( W_0(\eta) Z_n^{(2)}(\eta) \right) + \left\langle \left( W_0 Z_n^{(2)} \right), Z_m^{(1)} \right\rangle \\
&\quad - \left\langle \left( \int_z^\eta \left( 2U_0 \frac{\partial Z_n^{(2)}}{\partial x} \right) d\xi + \frac{\partial}{\partial x} \left[ \int_z^\eta \left( U_0 Z_n^{(2)} \right) d\xi \right] \right), Z_m^{(1)} \right\rangle, \\
G1_{mn} &= - \left\langle \left( \int_z^\eta \left( U_0 Z_n^{(2)} \right) d\xi \right), Z_m^{(1)} \right\rangle, \quad H1_m = g \left\langle 1, Z_m^{(1)} \right\rangle,
\end{align*}
$$

and

$$A2_n = U_0, \quad B2_n = \frac{\partial Z_n^{(2)}}{\partial x}(\eta), \quad C2_n = Z_n^{(2)}(\eta).$$

In the coefficients expressed above, the values of the functions $Z_n^{(0)}$, $Z_n^{(1)}$, $Z_n^{(2)}$, and $Z_n^{(3)}$ at the free surface elevation $\eta$ are often used. Using the definition of these functions, it is easily seen that for all of them,

$$f(z = \eta; h, \eta) = f(z = 0; h, \eta) + \frac{\partial f}{\partial z}(z = 0; h, \eta)\eta + O(\eta^2)$$

and the following approximation is used for vertical integrals as concerns the treatment of implicit nonlinearities.:

$$\int_z^\eta f(\xi)d\xi = \int_z^0 f(\xi)d\xi + f(0)\eta + O(\eta^2)$$

In the present work, a weakly nonlinear system is obtained, and the above second-order terms in $\eta$, which will eventually result in third-order terms in $(u, \eta)$, will be discarded, maintaining consistency with the order of approximation of the explicit nonlinearities. This significantly facilitates the calculation of the system coefficients at the pre-processing stage. On the other hand, as already mentioned, the explicit nonlinearities are retained in the system and the nonlinear terms $NL_m^{(1)}$ and $NL_m^{(2)}$ are as follows:

$$NL_m^{(1)} = \left\langle E_3, Z_m^{(1)} \right\rangle + \left\langle \frac{\partial F_3}{\partial x}, Z_m^{(1)} \right\rangle - \left\langle 1, Z_m^{(1)} \right\rangle \frac{1}{2} \frac{\partial}{\partial x} \left[ u(\eta)^2 + w(\eta)^2 \right], \\ \text{and } NL_m^{(2)} = -u(\eta) \frac{\partial \eta}{\partial x}.$$

## References

1. Rey, V.; Charland, J.; Touboul, J. Wave—Current interaction in the presence of a 3D bathymetry: Deep water wave focusing in opposite current conditions. *Phys. Fluids* **2014**, *26*, 096601. [CrossRef]
2. Haas, K.A.; Svendsen, I.A. Laboratory measurements of the vertical structure of rip currents. *J. Geophys. Res.* **2002**, *107*, 15-1-15-19.
[CrossRef]

3. Furgerot, L.; Sentchev, A.; Bailly-du-bois, P.; Lopez, G.; Morillon, M.; Poizot, E.; Méar, Y.; Bennis, A.C. One year of measurements in Alderney Race: Preliminary results from database analysis. *Phil. Trans. R. Soc. A* **2020**, *378*, 20190625. [CrossRef] [PubMed]

4. Kharif, C.; Abid, M. Miles theory revisited with constant vorticity in water of infinite depth. *J. Mar. Sci. Eng.* **2020**, *8*, 623.
[CrossRef]

5. Zhang, X.; Simons, R.; Zheng, J.; Zhang, C. A review of the state of research on wave-current interaction in nearshore areas. *Ocean Eng.* **2022**, *243*, 110202. [CrossRef]

6. Nwogu, O.G. Interaction of finite-amplitude waves with vertically sheared current fields. *J. Fluid Mech.* **2009**, *627*, 179-213.
[CrossRef]

7. Zhang, J.S.; Zhang, Y.; Jeng, D.S.; Liu, P.L.F.; Zhang, C. Numerical simulation of wave-current interaction using a RANS solver. *Ocean Eng.* **2014**, *75*, 157-164. [CrossRef]

8. Gao, J.; Zhou, X.; Zhou, L.; Zang, J.; Chen, H. Numerical investigation on effects of fringing reefs on low-frequency oscillations within a harbor. *Ocean Eng.* **2019**, *172*, 86-95. [CrossRef]

9. Marino, M.; Faraci, C.; Musumeci, R.E. Shoaling waves interacting with an orthogonal current. *J. Mar. Sci. Eng.* **2020**, *8*, 281.
[CrossRef]

10. Kang, A.; Zhu, B.; Lin, P.; Ju, J.; Zhang, J.; Zhang, D. Experimental and numerical study of wave-current interactions with a dumbbell-shaped bridge cofferdam. *Ocean Eng.* **2020**, *210*, 107433. [CrossRef]

11. Yang, Z.; Huang, B.; Kang, A.; Zhu, B.; Han, J.; Yin, R.; Li, X. Experimental study on the solitary wave-current inter action and the combined forces on a vertical cylinder. *Ocean Eng.* **2021**, *236*, 109569. [CrossRef]

12. Faraci, C.; Musumeci, R.; Marino, M.; Ruggeri, A.; Carlo, L.; Jensen, B.; Foti, E.; Barbaro, G.; Elsaßer, B. Wave- and current-dominated combined orthogonal flows over fixed rough beds. *Cont. Shelf Res.* **2021**, *220*, 104403. [CrossRef]

13. Massel, S. Extended refraction-diffraction equations for surface waves. *Coast. Eng.* **1993**, *19*, 97-126. [CrossRef]

14. Chamberlain, P.G.; Porter, D. The modified mild-slope equation. *J. Fluid Mech.* **1995**, *291*, 393-407. [CrossRef]

15. Raoult, C.; Benoit, M.; Yates, M. Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments. *Coast. Eng.* **2016**, *114*, 194-207. [CrossRef]

16. Belibassakis, K.A.; Athanassoulis, G.A. A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions. *Coast. Eng.* **2011**, *58*, 337-350. [CrossRef]

17. Booij, N. Gravity Waves on Water with Non-Uniform Depth and Current. Ph.D. Thesis, Delft University of Technology, Delft, The Netherlands, 1981.

18. Liu, P.L.F. Wave-current interactions on a slowly varying topography. *J. Geophys. Res.* **1983**, *88*, 4421-4426. [CrossRef]

19. Kirby, J.T. A note on linear surface wave-current interaction over slowly varying topography. *J. Geophys. Res.* **1984**, *89*, 745-747.
[CrossRef]

20. Yang, Z.; Liu, P. Depth-integrated wave-current models. Part 1. Two-dimensional formulation and applications. *J. Fluid Mech.* **2020**, *883*, A4. [CrossRef]

21. Duan, W.Y.; Wang, Z.; Zhao, B.B.; Ertekin, R.C.; Yang, W.Q. Steady solution of solitary wave and linear shear current interaction. *Appl. Math. Model.* **2018**, *60*, 354-369. [CrossRef]

22. Touboul, J.; Charland, J.; Rey, V.; Belibassakis, K. Extended Mild-Slope equation for surface waves interacting with a vertically sheared current. *Coast. Eng.* **2016**, *116*, 77-88. [CrossRef]

23. Belibassakis, K.; Touboul, J.; Laffitte, E.; Rey, V. A Mild-Slope System for Bragg Scattering of Water Waves by Sinusoidal Bathymetry in the Presence of Vertically Sheared Currents. *J. Mar. Sci. Eng.* **2019**, *7*, 9. [CrossRef]

24. Belibassakis, K.; Touboul, J. A Nonlinear Coupled-Mode Model for Waves Propagating in Vertically Sheared Currents in Variable Bathymetry—Collinear Waves and Currents. *Fluids* **2019**, *4*, 61. [CrossRef]

25. Rienecker, M.; Fenton, J. A Fourier approximation method for steady water waves. *J. Fluid Mech.* **1981**, *104*, 119-137. [CrossRef]

26. Francius, M.; Kharif, C. Two-dimensional stability of finite-amplitude gravity waves on water of finite depth with constant vorticity. *J. Fluid Mech.* **2017**, *830*, 631-659. [CrossRef]

27. Coddington, E.A.; Levinson, N. *Theory of Ordinary Differential Equations*; McGraw Hill: New York, NY, USA, 1955.

28. Boyd, J. *Chebyshev and Fourier Spectral Methods*; Dover Publications: New York, NY, USA, 2001.

29. Dingemans, M. *Water Wave Propagation over Uneven Bottoms*; World Scientific: Singapore, 1997.

30. Okamoto, H.; Shōji, M. *The Mathematical Theory of Permanent Progressive Water-Waves*; Advanced Series in Nonlinear Dynamics; World Scientific: Singapore, 2001.

31. Li, Y.; Ellingsen, S. A Framework for Modeling Linear Surface Waves on Shear Currents in Slowly Varying Waters. *J. Geophys. Res.* **2019**, *124*, 2527-2545. [CrossRef]

32. Domermuth, D.G. The initialization of nonlinear waves using an adjustment scheme. *Wave Motion* **2000**, *32*, 307-317. [CrossRef]

33. Mei, C.-C. Resonant reflection of surface water waves by periodic sandbars. *J. Fluid Mech.* **1985**, *152*, 315-335. [CrossRef]

34. Tang, H.-J.; Huang, C.-C. Bragg reflection in a fully nonlinear numerical wave tank based on boundary integral equation method. *Ocean Eng.* **2008**, *35*, 1800-1810. [CrossRef]

35. Peng, J.; Tao, A.-F.; Fan, J.; Zheng, J.-H.; Liu, Y.-M. On the Downshift of Wave Frequency for Bragg Resonance. *China Ocean Eng.* **2022**, *36*, 76-85. [CrossRef]
36. Gao, J.; Ma, X.; Dong, G.; Chen, H.; Liu, Q.; Zang, J. Investigation on the effects of Bragg reflection on harbor oscillations. *Coast. Eng.* **2021**, *170*, 103977. [CrossRef]

37. Gao, J.; Shi, H.; Jun Zang, J.; Liu, Y. Mechanism analysis on the mitigation of harbor resonance by periodic undulating topography. *Ocean Eng.* **2023**, *281*, 114923. [CrossRef]

38. Raj, R.; Guha, A. On Bragg resonances and wave triad interactions in two-layered shear flows. *J. Fluid Mech.* **2019**, *867*, 482–515. [CrossRef]

39. Laffitte, E.; Rey, V.; Touboul, J.; Belibassakis, K. Water wave scattering by sinusoidal bed in the presence of vertically sheared current. *Appl. Ocean Res.* **2021**, *108*, 102549. [CrossRef]

40. Chen, H.; Zou, Q. Effects of following and opposing vertical current shear on nonlinear wave interactions. *Appl. Ocean Res.* **2019**, *89*, 23–35. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.