Kostas Belibassakis

Julien Touboul
email: julien.touboul@mio.osupytheas.fr

A Nonlinear Coupled-Mode Model for Waves Propagating in Vertically Sheared Currents in Variable Bathymetry-Collinear Waves and Currents

Keywords: coupled modes system, nonlinear water waves, wave-bathymetry interaction, wave-current interaction

A novel coupled-mode model is developed for the wave-current-seabed interaction problem, with application in wave scattering by non-homogeneous, sheared currents over general bottom topography. The formulation is based on a velocity representation defined by a series of local vertical modes containing the propagating and evanescent modes, able to accurately treat the continuity condition and the bottom boundary condition on sloping parts of the seabed. Using the above representation in Euler equations, a coupled system of differential equations on the horizontal plane is derived, with respect to the unknown horizontal velocity modal amplitudes. In the case of small-amplitude waves, a linearized version of the above coupled-mode system is obtained, and the dispersion characteristics are studied for various interesting cases of wave-seabed-current interaction. Keeping only the propagating mode in the vertical expansion of the wave potential, the present system is reduced to a one-equation, non-linear model, generalizing Boussinesq models. The analytical structure of the present coupled-mode system facilitates extensions to treat non-linear effects and further applications concerning wave scattering by inhomogeneous currents in coastal regions with general 3D bottom topography.

Introduction

In coastal areas, steep bathymetries and strong currents are often observed. Among several causes, the presence of cliffs, rocky beds, or human structures may cause strong variations of the sea bed, while oceanic circulation, tides, wind action, or wave breaking can be responsible for the generation of strong currents. For both coastal safety and engineering purposes, there is strong interest in providing efficient models predicting the nonlinear, phase-resolved behavior of water waves in such areas. This leads to computational difficulties, since large computational areas often need to be considered, and the evolution of water waves propagating in such inhomogeneous media involves multiple-scale phenomena. Moreover, many physical processes influence the dynamics of water waves in such cases, due to reflection, refraction, and diffraction of water waves, in conjunction with nonlinear wave-bottom, wave-current, and wave-vorticity interactions, and the description of these processes requires careful attention.

Various attempts to partially treat the above problems are found in the literature. Starting with depth-averaged models, the mild-slope equation [START_REF] Dingemans | Water Wave Propagation over Uneven Bottoms[END_REF][START_REF] Berkhoff | Computation of combined refraction-diffraction[END_REF] was derived to describe refraction, diffraction, and reflection of dispersive water waves over varying bathymetry. The derivation is based on the assumption that the vertical structure of the wave field is represented by hyperbolic cosine function, and thus, the mild-slope equation is restrained to the description of linear water waves evolving on slowly varying bathymetries with respect to the wavelength. Similarly, Boussinesq models [START_REF] Dingemans | Water Wave Propagation over Uneven Bottoms[END_REF][START_REF] Massel | Hydrodynamics of Coastal Zones[END_REF] have been derived by using a polynomial function. Later, Massel [START_REF] Massel | Extended refraction-diffraction equations for surface waves[END_REF] and Chamberlain and Porter [START_REF] Chamberlain | The modified mild-slope equation[END_REF] provided an extension of mild-slope models by considering higher order approximations with respect to the bed-slope parameter. Furthermore, Massel [START_REF] Massel | Extended refraction-diffraction equations for surface waves[END_REF] suggested a multi-modal expansion to improve the representation of the vertical distribution of the wave field deriving a set of coupled equations, involving both propagating and evanescent modes. The latter approach has been further improved by Athanassoulis and Belibassakis [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF] by including extra terms to consistently treat the boundary conditions. This approach has been extended to nonlinear equations [START_REF] Athanassoulis | New evolution equations for non-linear water waves in general bathymetry with application to steady travelling solutions in constant, but arbitrary depth[END_REF][START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF], and many solvers are now based on similar techniques, e.g., Raoult et al. [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF]. The above techniques have demonstrated their efficiency in modeling wave-bottom interactions in inhomogeneous environments. Similar developments for modelling wave-current interactions have been presented by various authors, e.g., Booij [START_REF] Booij | Gravity Waves on Water with Non-Uniform Depth and Current[END_REF], Liu [START_REF] Liu | Wave-current interactions on a slowly varying topography[END_REF], and Kirby [START_REF] Kirby | A note on linear surface wave-current interaction over slowly varying topography[END_REF]. The current considered in the above works is assumed to be slowly varying in horizontal directions and uniform in depth. In the same direction, coupled-mode models have also been derived [START_REF] Belibassakis | A coupled-mode model for water wave scattering by horizontal, non-homogeneous current in general bottom topography[END_REF].

More recently, it was established that the vorticity involved in vertically non-uniform currents could have significant effects on the propagation of water waves [START_REF] Laffitte | Wave-bottom-current interaction, effects of the wave vorticity on the Bragg resonance[END_REF]. This gave rise to new derivations of models aiming to describe this interaction. In the case of currents with locally constant shear, extended mild slope models, similar to those derived by Berkhoff and Kirby, have been obtained [START_REF] Rey | Wave-Current interaction in the presence of a 3d bathymetry: Deep water wave focusing in opposite current conditions[END_REF][START_REF] Touboul | Extended Mild-Slope equation for surface waves interacting with a vertically sheared current[END_REF], and at the next stage, a coupled mode model was derived [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF]. The importance of this approach was recently validated experimentally in a previous study [START_REF] Belibassakis | A Mild-Slope System for Bragg Scattering of Water Waves by Sinusoidal Bathymetry in the Presence of Vertically Sheared Currents[END_REF]. In these works, the currents were assumed to vary linearly with depth. If this case is realistic in some specific configurations [START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF][START_REF] Nwogu | Interaction of finite-amplitude waves with vertically sheared current fields[END_REF], it might be suggested that when the current contains strong vorticity it will have more complicated vertical distribution. In the case when the second vertical derivative of the mean current flow is zero, the vorticity conservation equation for water waves involves no source term [START_REF] Nwogu | Interaction of finite-amplitude waves with vertically sheared current fields[END_REF]. Furthermore, in the case of mild horizontal variations of the shear, its effect can be approximately neglected, and water waves can be approximately described as irrotational [START_REF] Belibassakis | A Mild-Slope System for Bragg Scattering of Water Waves by Sinusoidal Bathymetry in the Presence of Vertically Sheared Currents[END_REF].

The present work aims to extend the validity of the above models to more general current configurations. Indeed, if the vorticity conservation equation involves a source term, due to the interaction of the flow with the sheared current, water waves cannot be treated as irrotational anymore. For this purpose, the present novel coupled-mode system is obtained in the framework of Euler equations by using a velocity-based formulation, allowing one to describe wave-vorticity interactions in more generic conditions.

The present paper is structured as follows. After the introduction, the modal expansion of the velocity field is introduced in Section 2 using a Fourier-type basis generated by local vertical eigenproblems. In particular, the present formulation is based on a velocity representation defined by a series of local vertical modes containing the propagating and evanescent modes, able to accurately treat the continuity condition and the bottom boundary condition on sloping parts of the seabed. Using the above representation in Euler equations, in Section 3 a coupled system of differential equations on the horizontal plane is derived, with respect to the unknown horizontal velocity modal amplitudes. Subsequently, Section 5 is dedicated to the analysis of the dispersive behavior of this system under linearity assumptions, and finally, numerical examples are presented and discussed in Section 6 in order to illustrate the ability of the present system to reproduce classical propagation configurations. Finally, the main conclusions are presented, including directions for further research.

Vertical Expansion of the Velocity Field

In the present work we will restrict ourselves to 2D problems corresponding to waves propagating in a vertical strip under the additional effects by collinear currents. However, the present analysis can be extended in two horizontal dimensions, and this is left to be considered along with 3D applications in future work.

Let D(t) denote the flow domain bounded below by the bottom boundary defined by the depth function, z = -h(x), and above by the free-surface defined by the corresponding elevation, z = η(x; t); see Figure 1. We assume that η(x; t) + h(x) > 0. In this work, we exploit a local-mode series expansion of the wave velocity u(x, z, t) = (U(x, z, t), W(x, z, t)) in variable bathymetry regions. It is based on extension of similar developed models, under the assumption of irrotationality, for the wave potential by Athanassoulis and Belibassakis [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF] for the linearized problem and subsequently for non-linear wave problems [START_REF] Athanassoulis | New evolution equations for non-linear water waves in general bathymetry with application to steady travelling solutions in constant, but arbitrary depth[END_REF][START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF]. This expansion has the general form of summation of local modes (each one numbered by an integer index). In particular, the local-mode representation of the horizontal velocity component reads as follows:

U(x, z, t) = ∑ n=0 U n (x, t)Z (1) n (z; h(x), η(x, t)), (1)
where Z

n (z; h(x), η(x, t)) denote vertical functions that are parametrically dependent on the bathymetry and the free-surface elevation, having the properties to satisfy the continuity equation, the bottom boundary condition, and to be complete in local vertical intervals z ∈ [-h(x), η(x, t)]. More specifically, the functions Z (1) n (z; h(x), η(x, t)) are defined as follows

Z (1) n (z; h, η) = cos[k n (z + h)] cos[k n (η + h)] , n = 0, 1, 2, 3, . . . ( 2
)
where the quantities k n = k n (h, η), n = 0, 1, 2 . . . , are obtained as the roots of the (dispersion-like) equations

µ 0 + k n tan[k n (h + η)] = 0 (3)
The first root is imaginary k 0 = i|k 0 | and the rest k n , n = 1, 2 . . . , are real, with k n ≈ nπ/(h + η), n → ∞ . The completeness property of the vertical basis stems from the fact that Z

(1) n (z) are generated from regular vertical Sturm-Liouville eigenproblems in the local vertical z ∈ [-h(x), η(x, t)], controlled by the frequency-type parameter µ 0 ; see also Belibassakis and Athanassoulis [START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF]. The latter parameter (or in non-dimensional form µ 0 h 0 , with h 0 denoting a characteristic depth) is not subjected to a priori restrictions. However, the specific values of these parameters influence the accuracy of truncated versions of the local-mode series, as will be illustrated in the sequel concerning the convergence of the expansion and the dispersion characteristics of the present model.

The vertical component W(x, z, t) of the velocity field is represented by a corresponding modal series based on a new set of vertical functions Z

n (z; h(x), η(x, t)), as follows

W(x, z, t) = -∑ n=0 ∂ x U n (x; t)Z (2) n (z; h; η) , (4)
and thus, it is expressed as follows

W(x, z, t) =-∑ n=0 ∂ x U n (x; t)Z (2) n (z; h; η)+U n (x; t)∂ x Z (2) n (z; h; η), ( 5
)
where the new vertical functions Z

n (z; h, η) are defined to automatically satisfy the wave kinematics, which are defined by continuity condition in the water column

∂ x U(x, z, t) + ∂ z W(x, z, t) = 0, (6)
and the bottom boundary condition

W(x, z = -h, t) + ∂ x hU(x, z = -h, t) = 0. (7)
It is easily verified that the above constraints are fulfilled by selecting the functions Z

n (z; h, η) to satisfy the requirement

∂ z Z (2) n (z) = Z (1) n (z), (8)
leading to the following analytic expression of the new set of vertical functions

Z (2) n (z; h, η) = 1 k n sin[k n (z + h)] cos[k n (η + h)] , n = 0, 1, 2, . . . . (9)
In fact, from Equations ( 5) and [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF] we see that the bottom vertical velocity is given by

W(x, z = -h, t) = -∑ n=0 U n (x; t)∂ x Z (2)
n (z = -h; h; η) [START_REF] Booij | Gravity Waves on Water with Non-Uniform Depth and Current[END_REF] The latter used in the bottom boundary condition, Equation ( 6), it finally leads to

W(x, z = -h, t) + h x U(x, z = -h, t) = -∑ n=0 U n (x; t)∂ x Z (2)
n (z = -h) + ∂ x h ∑ n=0 U n (x, t)Z (1) n (z = -h) = -∑ n=0 U n (x; t) ∂ x Z (2)
n (z = -h) -∂ x hZ (1)
n (z = -h) = 0, (11)
since, from Equation ( 8), it is obtained that

∂ x Z (2)
n (z = -h) = ∂ x hZ (1)
n (z = -h). Details concerning the latter results are provided in Appendix A.

Moreover, the continuity condition, for all x, and -h(x) < z < η(x, t), is identically satisfied

∂ x U(x, z, t) + ∂ z W(x, z, t) = ∑ n=0 ∂ x U n (x, t)Z (1) n (z; h, η) -∂ z ∑ n=0 ∂ x U n (x; t)Z (2) n (z; h, η) = ∑ n=0 ∂ x U n (x, t)Z (1) n (z; h, η) -∑ n=0 ∂ x U n (x; t)∂ z Z (2) n (z; h, η) = 0. ( 12
)
The usefulness of the above local-mode representations, Equations ( 1) and [START_REF] Chamberlain | The modified mild-slope equation[END_REF], is that used in the momentum equations will lead to a non-linear coupled-mode system of differential equations on the horizontal plane, with respect to unknown modal amplitudes U n (x, t) and the unknown free-surface elevation η(x, t). The coupled-mode system, in conjunction with the fast convergence properties of the local mode series, greatly facilitates the numerical solution of the present problem and will be presented in the next section. Similar modal-type series expansion has been previously introduced by Nadaoka et al. [START_REF] Nadaoka | A fully dispersive weakly nonlinear model for water waves[END_REF] for the development of a fully dispersive, weakly nonlinear, multiterm-coupling model for water waves, with application to slowly varying bottom topography. In that work, the vertical modes have been selected to have the form cosh(k n (z + h)) cosh -1 (k n h), where the parameters k n > 0 are independent from the free surface elevation η(x, t). This approach can be considered as an extension of Fourier methods above a flat bottom [START_REF] Baldock | Numerical calculations of large transient water waves[END_REF][START_REF] Tao | Nonlinear progressive waves in water of finite depth-An analytic approximation[END_REF] to variable bathymetry regions characterized by a mildly sloped bottom.

surface elevation   , x t  . This approach can be considered as an extension of Fourier methods above a flat bottom [START_REF] Baldock | Numerical calculations of large transient water waves[END_REF][START_REF] Tao | Nonlinear progressive waves in water of finite depth-An analytic approximation[END_REF] to variable bathymetry regions characterized by a mildly sloped bottom. As it is presented and discussed in more detail in Appendix B, the expected rate of convergence of the modal series (1) is

    2 , n U x t O n  
, where n is the mode number. Under the assumption of smooth upper   , x t  and lower   h x boundaries, the above rate ensures that

      1 , ; U x z t H D t 
. The latter is sufficient for the solution of the wave flow problem by a spectral Galerkin-type method that will be presented in the next section, as obtained by projecting the momentum flow equations on the vertical basis The modal representation Equation (1) can be further improved by introducing additional modes in the series expansion, as presented by Athanassoulis and Belibassakis [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF][START_REF] Athanassoulis | New evolution equations for non-linear water waves in general bathymetry with application to steady travelling solutions in constant, but arbitrary depth[END_REF] for the wave potential [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with non-planar boundaries[END_REF]. In the latter works the additional slopping-bottom and free-surface modes treat the inconsistency between the boundary conditions of the potential at Two examples illustrating the above results in the case of irrotational wave-like field under a sinusoidal upper surface of large amplitude and a vertically sheared Couette flow between horizontal planes are presented and discussed in Appendix B. As it is presented and discussed in more detail in Appendix B, the expected rate of convergence of the modal series (1) is U n (x, t) = O n -2 , where n is the mode number. Under the assumption of smooth upper η(x, t) and lower h(x) boundaries, the above rate ensures that U(x, z; t) ∈ H 1 (D(t)). The latter is sufficient for the solution of the wave flow problem by a spectral Galerkin-type method that will be presented in the next section, as obtained by projecting the momentum flow equations on the vertical basis Z n (z; h(x), η(x, t)), n = 0, 1, 2 . . . , for each instant and at each horizontal x-position. The modal representation Equation (1) can be further improved by introducing additional modes in the series expansion, as presented by Athanassoulis and Belibassakis [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF][START_REF] Athanassoulis | New evolution equations for non-linear water waves in general bathymetry with application to steady travelling solutions in constant, but arbitrary depth[END_REF] for the wave potential [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with non-planar boundaries[END_REF]. In the latter works the additional slopping-bottom and free-surface modes treat the inconsistency between the boundary conditions of the potential at z = -h(x), z = η(x, t), and the boundary conditions satisfied by the Sturm-Liouville problem generating the vertical basis. This accelerates the convergence of the modal series by increasing the rate of decay of the potential amplitudes from ϕ n (x, t) = O n -2 to ϕ n (x, t) = O n -4 . In the present case, there is also an inconsistency for the velocity field U(x, z; t) but only at the free surface z = η(x, t), and could be similarly treated accelerating the convergence. However, it should be mentioned that application of the present model to the special case of irrotational flows, without any additional modes for the expansions of U(x, z; t) and W(x, z; t), leads to a decay rate U n (x, t) = O n -2 , which would correspond to ϕ n (x, t) = O n -3 concerning the underlying potential modes, as obtained by direct integration.

          1 ; , , , 0,1, 2.

Momentum Equations

Two examples illustrating the above results in the case of irrotational wave-like field under a sinusoidal upper surface of large amplitude and a vertically sheared Couette flow between horizontal planes are presented and discussed in Appendix B.

Momentum Equations

In the previous section it is shown that through the appropriate selection of vertical expansions of the velocity field, the kinematics of the wave problem concerning continuity and bottom surface boundary condition are analytically satisfied. Thus, we are left with the requirement of the momentum equations and the free surface boundary conditions, including the kinematic and the dynamic ones.

Neglecting the effect of viscous dissipation in the fluid volume, the momentum equations in the horizontal and vertical components are, respectively,

∂ t U + 1 2 ∂ x |u| 2 -E = - ∂ x p ρ and ∂ t W + 1 2 ∂ z |u| 2 -e = - ∂ z (p + ρgz) ρ , ( 13
)
where |u| 2 = U 2 + W 2 and e = u × ω = (E, e), an helicity-type field, with ω denoting the vorticity ω = ∇ × u = Ω ĵ, which in the present two-dimensional case is fixed in the transverse direction (with ĵ denoting the corresponding unit vector). Thus, in the present 1DH case,

E = -WΩ, e = UΩ. (14)
Subsequently, we may consider at each horizontal x-position, the vertical momentum Equation ( 13) integrated in the vertical direction from a depth level z up to the free surface elevation η(x, t), from which an expression for the pressure in the local column is obtained

- p + ρgz ρ z=η(x,t) z = - η(x,t) z e(z)dz + η(x,t) z ∂ t Wdz + 1 2 U 2 + W 2 z=η(x,t) z . ( 15
)
Employing in Equation ( 15) the dynamic free-surface boundary condition

p(x, z = η(x, t), t) = 0, (16)
we obtain

p(x, z, t) ρ = g(η(x, t) -z) - η(x,t) z e(z)dz + η(x,t) z ∂ t Wdz + 1 2 U 2 + W 2 z=η(x,t) z . ( 17
)
The above equation is then differentiated with respect to the horizontal coordinate to obtain

- ∂ x p(x, z, t) ρ = -g∂ x η(x, t) + ∂ x η(x,t) z e(z)dz -∂ x η(x,t) z ∂ t Wdz - 1 2 ∂ x |u| 2 z=η(x,t) z , (18)
which substituted back to Equation ( 13) finally results in

∂ t U + 1 2 ∂ x |u| 2 -E = -g∂ x η(x, t) + ∂ x η(x,t) z e(z)dz -∂ x η(x,t) z ∂ t Wdz -1 2 ∂ x |u| 2 z=η(x,t) -|u| 2 or ∂ t U -E = -g∂ x η + ∂ x F -∂ x η(x,t) z ∂ t Wdz -[U∂ x U + W∂ x W] z=η(x,t) , (19a)
for all x and h(x) < z < η(x, t), where F(x, z, t) = η(x,t) z e(z)dz. The last Equation (19a), in conjunction with the kinematic free-surface boundary condition,

∂ t η + U(x, z = η, t) ∂ x η -W(x, z = η, t) = 0, (19b)
constitute the remaining system to be solved. In the next section, Equations (19a) and (19b) in conjunction with the representations defined by Equations ( 1) and (4) will be put in the form of a coupled-mode system with respect to the wave velocity modal amplitudes U n (x; t) and the free surface elevation η(x, t), constituting the unknown fields in the horizontal domain.

Irrotational Waves Over Variable Bathymetry

We first consider the case of irrotational waves without a current. In order to reduce it to a coupled system of equations we exploit the completeness properties of the local eigenfunctions

Z (1)
m (z; h(x), η(x, t)), m = 0, 1, 2, . . . and project both sides of Equation (19a) into the latter vertical basis. Thus, for each horizontal x-position, we obtain

∂ t U + ∂ x η(x,t) z ∂ t Wdz, Z (1) m η -h + g∂ x η + [U∂ x U + W∂ x W] z=η(x,t) 1, Z (1) m η -h = 0, m = 0, 1, 2, . . . , where f 1 (z), f 2 (z) η -h = z=η z=-h f 1 (z) f 2 (z)dz.
Furthermore, using the local mode representations defined by Equations ( 1) and ( 4), and carrying out the algebra, the following set of equations is derived

Z (1) m 2 ∂ t U n + ∑ n W (1) n , Z (1) m U n ∂ t η + 1, Z (1) m g∂ x η + [U∂ x U + W∂ x W] z=η(x,t) + ∂ x η z ∂ t Wdz, Z (1) m η -h = 0, m = 0, 1, 2, . . . , (20a)
where the orthogonality properties of the vertical basis

Z (1) n , Z (1) m η -h = δ nm Z
(1) m 2 with δ nm denoting Kronecker's delta, and

W (1) n (z; h, η) = ∂Z (1)
n (z; h, η)/∂η have been also used. Moreover, using the fact that the vertical local eigenfunctions Z [START_REF] Dingemans | Water Wave Propagation over Uneven Bottoms[END_REF] n are normalized by their free surface values (i.e., Z

n (z = η) = 1), the kinematic free-surface boundary condition, Equation (19b), takes the following form:

∂ t η + ∂ x η ∑ n=0 U n (x; t) + ∑ n=0 ∂ x U n (x; t)Z (2) n (z = η) = 0. ( 20b
)
Recalling the fact that the expected rate of decay of modal amplitudes U n (x; t) = O n -2 , specific attention is needed concerning the convergence of the above term-wise differentiated series at z = η representing the horizontal flow velocity derivatives. In particular, under the assumption that ∂ x U n (x; t) are expected to decay as ∂ x U n (x; t) = O n -2 based on evidence that all modes present similar wavelike behavior U n (x; t) ≈ αk -2 n k 0 2 U 0 (x; t), n > 1, which will be illustrated in the examples below, and using the fact that ∂ x Z

(1)

n (z = η) = O(1)
, which is obtained by straightforward algebra (see also [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with non-planar boundaries[END_REF]), the convergence of the series

∑ n=0 ∂ x U n (x; t)Z (1) n (z = η) + U n (x; t)∂ x Z (1) n (z = η) to the horizontal derivative ∂ x U(x, z = η; t) is obtained. Furthermore, using the fact that Z (2) n (x, z = η) = k -1 n tan(k n (η + h)) = µ 0 k -2 n and thus, Z (2)
n (x, z = η) = O n -2 , n → ∞ , the series appearing as the last term in Equation (20b) converges to W(x, z = η, t) = µ 0 ∑ n=0 ∂ x k -2 n U n (x; t) = µ 0 ∑ n=0 k -2 n ∂ x U n (x; t) -2k -3 n (∂ x k n )U n (x; t)
. Finally, using similar arguments, it can be shown that term-wise differentiation of the latter series converges to ∂ x W(x, z = η; t). Substitution of the above modal expansions for ∂ x U(x, z = η; t) and ∂ x W(x, z = η; t) in Equation ( 20), in conjunction with the representations Equations ( 1) and ( 4), eventually leads to the present nonlinear coupled-mode system (nCMS) formulated with respect to the horizontal flow velocity amplitudes U n (x; t) and the free surface elevation η(x, t) as the unknowns. The detailed analysis of the nCMS will be presented in future works.

Waves and Current Over Bathymetry

In the case of waves propagating under the additional effects of steady currents, the flow variables are split in two components: (i) the wave part (U, η), and (ii) the steady background current (U c , η c ). We furthermore assume that the variation of the current flow is smaller than the corresponding wave flow and that |∂ x η| >> |∂ x η c | ≈ 0. Using the above decomposition in Equation [START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF] and retaining terms to the lowest order [START_REF] Quinn | Explicit wave action conservation for water waves on vertically sheared flows[END_REF], these equations are put in the form

D t U -E = -g∂ x η + ∂ x F -∂ x η c (x,t) z ∂ t Wdz -[U∂ x U c + W∂ x W] z=η c (x,t) , ( 21a
)
D t η + U(x, z = η c , t) ∂ x η -W(x, z = η c , t) = 0, (21b)
where D t = ∂ t + U c ∂ x denotes differentiation moving with the current.

To reduce to a coupled system of equations, we again exploit the completeness properties of the local eigenfunctions Z

n (z; h(x), η c (x, t)), n = 0, 1, 2, . . . and project both sides of Equation (21a) into the latter vertical basis. Thus, for each horizontal x-position, we obtain

D t U -E + g∂ x η + [U∂ x U c + W∂ x W] z=η c (x,t) -∂ x F + ∂ x η c (x,t) z ∂ t Wdz, Z (1) m η c -h = 0, m = 0, 1, 2, . . . , where now f 1 (z), f 2 (z) η c -h = z=η c z=-h f 1 (z) f 2 (z)
dz. Furthermore, using in Equation ( 21) the local mode representations defined by Equations ( 1) and ( 4), and carrying out the algebra, the following set of equations is derived, assuming also that the vorticity is essentially contained in the current

Ω ≈ Ω c , Z (1) m 2 D t U m + 1, Z (1) m η ψ -h g∂ x η + [U∂ x U c + W∂ x W] z=ηc(x,t) + ΩW -∂ x ηc(x,t) z (UΩ -D t W)dz, Z (1) m ηc -h + ∑ n=0 U n D t Z (1) n , Z (1) m ηc -h = 0, m = 0, 1, 2, . . . . (22)
Moreover, using the fact that the vertical local eigenfunctions Z

n (z = η c ) = 1), the kinematic free-surface boundary condition takes the form:

D t η + ∂ x η ∑ n=0 U n (x; t) + ∑ n=0 ∂ x U n (x; t)Z (2) n (z = η c ) = 0. ( 23
)
The above Equations ( 22) and ( 23) constitute the present CMS in the case of waves and currents, again formulated with respect to the horizontal flow amplitudes U n (x; t) and the free surface elevation η(x, t) as unknowns.

Dispersion Characteristics of the Linearized CMS

The set of Equations ( 22) and ( 23) is linearized by first dropping the terms involving products of the unknowns. Furthermore, assuming that η c ≈ 0, the vertical integrals are restricted up to mean water level z = 0. In addition, assuming that the vorticity is essentially contained in the current Ω ≈ Ω c , we obtain the following expressions of the corresponding terms in Equation ( 22)

Ω c W, Z (1) m 0 -h = -∑ n U n Ω c ∂ x Z (2) n , Z (1) m 0 -h + ∂ x U n Ω c Z (2) n , Z (1) m 0 -h , -∂ x z=0 z (UΩ c -D t W)dz = -∑ n ((∂ x U n )F n + U n ∂ x F n ) -∑ n Z (3) n ∂ 2 x (D t U n ) + 2∂ x (D t U n )∂ x Z (3) n + (D t U n )∂ 2 x Z (3) n
, where

F n = z=0 z Ω c Z (1)
n dz, and the functions Z

(1)

n (z), Z (2)
n (z) are obtained from Equations ( 1) and (4) using η = 0. Moreover, in the present case Z

(3)

n (z) = 0 z Z (2)
n (z)dz and the brackets denote

f 1 (z), f 2 (z) 0 -h = z=0 z=-h f 1 (z) f 2 (z)dz.
Thus, Equations ( 22) are put in the following form:

Z (1) m 2 D t U m + g 1, Z (1) m 0 -h ∂ x η -∑ n=0 Z (3) n , Z (1) m ∂ 2 x (D t U n ) + 2 ∂ x Z (3) n , Z (1)
m 0 -h ∂ x (D t U n ) + ∂ 2 x Z (3) n , Z (1)
m 0 -h (D t U n ) -∑ n=0 ∂ x U n F n , Z (1)
m 0 -h + Ω c Z (2) n , Z (1)
m 0 -h + U n ∂ x F n , Z (1)
m 0 -h + Ω c ∂ x Z (2) n , Z (1)
m 0 -h -U c ∂ x Z (1)
n (0) 1, Z (1) m 0 -h = 0, m = 0, 1, 2, . . . . (24a)
Also, the kinematic free-surface condition, Equation [START_REF] Tao | Nonlinear progressive waves in water of finite depth-An analytic approximation[END_REF], is linearized by dropping the second term and transferring the condition to the level z = 0

D t η + ∑ n=0 ∂ x (c n U n ) = 0, ( 24b
)
where

c 0 = Z (2) 0 (z = 0) = tanh(|k 0 |h)/|k 0 | and c n = Z (2) 0 (z = 0) = tan(k n h)/k n .
We next consider an environment characterized by constant parameters (constant depth, constant background flow) and use again the local mode representations. Assuming now harmonic time dependence and seeking periodic solutions U n (x; t) = U n (x)e i(kx-ωt) , η(x; t) = η(x)e i(kx-ωt) in a constant depth strip with range independent steady current U(x), we obtain from Equation (24b)

η = k (ω -U 0 k) ∑ n=0 c n U n (25)
Using the latter result in Equation (24a) we obtain the following homogeneous system

∑ n δ mn - k 2 g σ 2 α m c n + k 2 A mn + k σ B mn U n = 0 m = 0, 1, 2, . . . ( 26
)
where σ = ω -U c0 k, with U c0 = U c (z = 0) the value of the current at the surface. Moreover, δ mn is Kronecker's delta, and the rest coefficients are given by

α m = 1, Z (1)
m 0 -h Z (1) m -2 , A mn = Z (3) n , Z (1)
m 0 -h Z (1) m -2
and

B mn = F n + Ω c Z (2) n , Z (1)
m 0 -h Z (1) m -2 . ( 27
)
For a given wave frequency ω, the periodicity parameter k is found by requiring the determinant of Equation ( 26) to be zero, which enables us to investigate the dispersion properties of the linearized system and compare against analytical solutions. For this purpose, we truncate the local-mode expansion, Equation (1), keeping a finite number of terms, and investigate the linearized truncated CMS dispersion curve against known analytical solutions.

(i) No current: In the case of waves propagating in constant depth without current (σ = ω and B mn = 0), numerical results obtained as non-trivial solutions of Equation ( 26) concerning the normalized phase speed Ĉ(kh), which depends only on the non-dimensional wavenumber kh, are compared against the analytical solution, given by

Ĉ(kh) = C gh = tanh(kh)/(kh) (28)
Results are presented in Figure 2 for values of the nondimensional wavenumber kh from 0 to 2π, ranging from very shallow to deep wave conditions, and using two values of the parameter controlling the vertical basis: µ 0 h = 0.1π, 0.5π. A rapid convergence of results to the analytical solution is observed, shown by using thick lines, as the modes included in the present series increases. We also observe in Figure 2 that the present representation, keeping only the first term (n = 0) in the expansion, provides excellent results over a bandwidth of non-dimensional wavenumbers around the selected value of the parameter µ 0 h, indicating that a simplified model based on a single-term approximation of the flow velocity U(x, z; t) = U 0 (x; t)Z 0 (z; η(x; t), h(x)) is able to provide good results with appropriate selection of parameter µ 0 h.

   .     / ˆ, h h k C C S F gh gh     (29a)
More specifically, for each frequency ω, the corresponding wavenumber   k  is derived from the roots of the following dispersion relation

    2 tanh h h S F kh kh kh   (29b)
where the case of opposing waves to current correspond to (ii) Parallel flow: In the case of a uniform current in depth, the flow vorticity is zero (B mn = 0). Numerical results are obtained as non-trivial solutions of Equation ( 26) concerning the normalized phase speed Ĉ(S h , F h ), which now depends both on the frequency parameter S h = ω 2 h/g and the bathymetric Froude number F h = U c / gh, and compared against the analytical solution

Ĉ(S h , F h ) = C gh = ω/k(ω) gh (29a)
More specifically, for each frequency ω, the corresponding wavenumber k(ω) is derived from the roots of the following dispersion relation

S h -F h kh 2 = khtanh(kh) (29b)
where the case of opposing waves to current correspond to F h < 0 and following waves to F h > 0, respectively. Results concerning a relatively strong current characterized by F h = 0.1 for following waves are presented in Figure 3, again for all values of the non-dimensional wavenumber kh, ranging from 0 to 2π, corresponding from very shallow to deep wave conditions. The same as before, two values of the parameter controlling the vertical basis (µ 0 h = 0.1π, 0.5π) are used. Again, rapid convergence of present results to the analytical solution, which is shown by using thick lines, is observed as the modes included in the present series increases. The corresponding results in the case of opposing wave and current are presented in Figure 4. In this case the convergence to the analytical results is slower, especially when the propagation speed becomes smaller approaching the blocking condition. However, in both cases, the present representation, keeping only the first term (n = 0) in the expansion, provides good results over a bandwidth of non-dimensional wavenumbers around the selected value of the parameter µ 0 h, rendering the single-term model a good approximation. (iii) Current with constant vorticity: In the case of a linear vertical profile U c (z) = U c0 + Sz, the vorticity is constant in the domain Ω c = S and the coefficient B mn in Equation ( 26) becomes B mn = Sc n α m . Nontrivial solutions of Equation ( 26) are compared against the dispersion relation. Numerical results obtained as non-trivial solutions of Equation ( 26) concern the normalized phase speed Ĉ(S h , F h ). Following a previous study [START_REF] Touboul | Extended Mild-Slope equation for surface waves interacting with a vertically sheared current[END_REF], in this case the dispersion relation is formulated in terms of the frequency parameter S h = ω 2 h/g and two bathymetric Froude numbers: one corresponding to the current flow at the surface F h = U c0 / gh and a second one F * = (U c0 -2Sd)/ gh associated with the current speed at a depth z = -2d = -tanh(kh)/k below the free surface. The latter parameter is also expressed in terms of Froude number F h and the nondimensional shear coefficient S c = S 2 h/g, as follows: F * = F h -√ S c tanh(kh)/kh. Again, for each frequency ω, the corresponding wavenumber k(ω) is derived from the roots of the following dispersion relation

S h -F h kh S h -F * kh = khtanh(kh) (30a)
and are compared against the analytical solution

Ĉ(S h , F h , S c ) = C gh = ω/k(ω) gh . ( 30b
)
Results for following waves and a relatively strong current characterized by F h = 0.1 and positive shear S c = 0.03 are presented in Figure 5a,b, while the case of negative shear is presented in Figure 5c,d. Corresponding results for opposing waves and vertically sheared currents characterized by the same parameters are presented in Figure 6. The behavior of the present CMS concerning the dispersion is again verified. In all cases the single-term model offers a good approximation for a band of frequencies around the selected value of the parameter µ 0 h defining the vertical basis. Further investigation of the dispersive properties of the present CMS in cases corresponding to a current with a more general vertical profile and the simulation of waves will be the subject of future work.

Numerical Results

In this section, first numerical results concerning the simulation by the present CMS for waves propagating in variable bathymetry regions will be presented. We will first investigate in Sections 5.1 and 5.2 the behavior of the linearized CMS in cases of waves propagating in variable bathymetry regions without and with current effects, respectively, and provide some evidence concerning the convergence of the modal series. This will reveal the characteristics and rate of decay of the modal amplitudes   ; n U x t and will illustrate the significance of the first term n = 0 in the present modal expansion. Subsequently, the effect of nonlinearity will be demonstrated in Section 5.3, in the case of wave propagation over bathymetry without current, by using a simplified model based on the truncated local mode series, and keeping only the first n = 0 term.

Numerical Solution

The modal series Equations ( 1) and ( 4) are truncated, keeping only a finite number of terms, and the present CMS is discretized by using second-order, central finite differences to approximate horizontal and time derivatives of the unknowns

  ; n U x t and   ; x t  based on a uniform grid       1 , 1 i a n x x i x t n t        ,

Numerical Results

In this section, first numerical results concerning the simulation by the present CMS for waves propagating in variable bathymetry regions will be presented. We will first investigate in Sections 5.1 and 5.2 the behavior of the linearized CMS in cases of waves propagating in variable bathymetry regions without and with current effects, respectively, and provide some evidence concerning the convergence of the modal series. This will reveal the characteristics and rate of decay of the modal amplitudes U n (x; t) and will illustrate the significance of the first term n = 0 in the present modal expansion. Subsequently, the effect of nonlinearity will be demonstrated in Section 5.3, in the case of wave propagation over bathymetry without current, by using a simplified model based on the truncated local mode series, and keeping only the first n = 0 term.

Numerical Solution

The modal series Equations ( 1) and ( 4) are truncated, keeping only a finite number of terms, and the present CMS is discretized by using second-order, central finite differences to approximate horizontal and time derivatives of the unknowns U n (x; t) and η(x; t) based on a uniform grid (x i = x a + (i -1)∆x, t n = (n -1)∆t), defined by subdividing the horizontal extent of the domain ranging in x a ≤ x ≤ x b by using M x steps of equal length ∆x = (x bx a )/M x .

Based on the above analysis presented in Section 4 concerning the dispersion characteristics of the present model, in conjunction with extensive experience from similar applications of the CMS, in most cases the number of the retained modes required for convergence is 4-5, and the horizontal subdivision of the domain is based on using 30-40 points per wavelength. In the examples that will be presented and discussed below, the dimension of the discrete coupled-mode system is of the order of 1500, which is expected to be significantly (orders of magnitude) less than the analysis required for a solution of Euler or Navier-Stokes equations based on a full spatial grid (see Koutrouveli and Dimas [26]). Moreover, the time interval is 0 ≤ t ≤ t max by using constant time steps ∆t = t max /M t . An implicit Crank-Nicolson method for time-integration is considered, satisfying Courant numbers (C∆t/∆x) max < 1.

In the examples presented below the depth function h(x) changes smoothly from h a = h(x a ) to h b = h(x b ), such as the bottom slope and curvature are zero at the ends of the domain, the flow motion starts from rest, and boundary conditions are imposed on the wave inlet boundary x = x a corresponding to a simple periodic wave corresponding to the depth h a . Finally, the reflected back propagating waves in the vicinity of the wave entrance x = x a and the radiated waves in the wave exit x = x b are absorbed by using a relaxation scheme extending over one characteristic wave length in these regions.

An example concerning harmonic incident waves of period T = 2 s and a height of 5 cm propagating over a linear upslope with decreasing depth from h a = 0.4 m to h b = 0.1 m, presenting a bottom slope of 4.5% is presented in Figure 7. The calculated modes U n (x; t), n = 0, 1, 2, 3 at a given time instant are plotted along with the velocity field U(x, z = 0; t) and the free surface elevation η(x; t) at the same instant, as obtained by the present linearized CMS, using µ 0 h = 0.4. Using the results presented in Appendix B concerning the rate decay of the mode amplitude, Equations (A4) and (A5), in conjunction with the fact that all modes (n > 0) follow the phase of the propagating mode (n = 0) as it can be also observed in Figure 7, the following trend is expected to be valid:

U n (x; t) ≈ αk -2 n k 0 2 U 0 (x; t), n > 1.
The latter result has been used in the discussion concerning the convergence of the present modal expansion in Section 2. Results obtained by the present linearized CMS for harmonic incident waves of period T = 2 s and a height of 5 cm propagating over a linear upslope characterized by bottom slope 4.5%. The first 4 subplots present the calculated modes U n (x;t), n = 0, 1, 2, 3 at a given time instant. The vertical dashed lines indicate the absorbing (relaxation) zones. The last three subplots show (using thick lines) the velocity field U(x,z = 0;t) on the mean free surface (z = 0) and the free surface elevation η(x;t) at the same instant, and finally the sloping bottom topography.

Wave Propagation Over an Underwater Bar in the Presence of Opposing Shearing Current

As a second example, we consider the case of normal incident waves propagating over a trapezoidal bar, seated on a flat horizontal bottom in the presence of an opposing vertically sheared current, which is presented and discussed in Belibassakis et al. [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF]. In this case, experiments have also been carried out in the wave flume at SeaTech, University of Toulon, France. At the upwave end of the channel, an electromagnetic piston generates regular waves by horizontal motion. Also, at the downwave end, a sloping beach is used to absorb waves. An opposing current is injected in the channel by a hydraulic pump and a perforated screen is used to control the shear. The opposing current is adjusted to generate flow profiles that are linear in depth U c (x, z) = U 0 (x) + S(x)z, where U 0 (x) is the surface current and S(x) the shear, and thus, in this specific case the vorticity is Ω(x) = S(x).

In the experiments, the water depth was 0.305 m, and the incident wave frequencies were ranging from 0.65 Hz to 1.3 Hz with amplitudes between 1-2 cm ensuring small wave steepness, both without and with the vertically sheared current. The distributions of the bathymetry, as well as the surface current and shear data, are presented in Figure 8. [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF]. In the lower subplot, the surface current U 0 (x) is plotted by using a solid line and the shear S(x) = Ω(x) is shown by using a dashed line, respectively. Furthermore, a coupled-mode model has been developed and discussed in Belibassakis et al. [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF] for the above problem formulated in the frequency domain. The latter model is based on the assumption that the wave field is irrotational, which is a plausible approximation for cases where the current flow and the vorticity structure are simplified and the corresponding parameters are slowly varying horizontal functions. Results from the above coupled mode system have been compared in a previous study [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF] against experimental measurements concerning the reflection coefficient, illustrating that the model is able to provide good predictions.

In order to illustrate the performance of the present system, calculated results are presented in Figures 9 and 10 concerning harmonic incident waves of period T = 1 s and an amplitude of 2 cm propagating over the trapezoidal bar topography, without and with the consideration of the sheared current, using the data of Figure 8. In particular, in Figure 9 the calculated modes U n (x; t), n = 0, 1, at a given time instant are plotted along with the surface wave velocity field U(x, z = 0; t) and the free surface elevation η(x; t) at the same instant, as obtained by the present linearized CMS given by Equation ( 24), first without the consideration of any current. Again we observe the very fast decay of the mode amplitudes. Moreover, the present solution concerning the free surface elevation calculated at an instant where the wave has reached the absorbing beach is also compared with the time-harmonic solution indicated by using dashed lines, illustrating the compatibility of the present model with previously established ones. Furthermore, a coupled-mode model has been developed and discussed in Belibassakis et al. [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF] for the above problem formulated in the frequency domain. The latter model is based on the assumption that the wave field is irrotational, which is a plausible approximation for cases where the current flow and the vorticity structure are simplified and the corresponding parameters are slowly varying horizontal functions. Results from the above coupled mode system have been compared in a previous study [START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF] against experimental measurements concerning the reflection coefficient, illustrating that the model is able to provide good predictions.

In order to illustrate the performance of the present system, calculated results are presented in Figures 9 and 10 concerning harmonic incident waves of period T = 1 s and an amplitude of 2 cm propagating over the trapezoidal bar topography, without and with the consideration of the sheared current, using the data of Figure 8. In particular, in Figure 9 24), first without the consideration of any current. Again we observe the very fast decay of the mode amplitudes. Moreover, the present solution concerning the free surface elevation calculated at an instant where the wave has reached the absorbing beach is also compared with the time-harmonic solution indicated by using dashed lines, illustrating the compatibility of the present model with previously established ones.

The same behavior is also confirmed in the case of waves propagating against linear vertically sheared current, as shown in Figure 10, where the corresponding result is presented with the The same behavior is also confirmed in the case of waves propagating against linear vertically sheared current, as shown in Figure 10, where the corresponding result is presented with the modifications by the presence of the opposing sheared current and compared with the frequency-domain solution. We clearly observe the shortening of the wavelengths due to the effects of the current, and the steepening of waves, especially over the submerged structure.

Moreover, in this case, the disturbance of the flow due to combined effect of variable bathymetry and vorticity contained in the vertically sheared current is seen to generate a more complicated pattern concerning the higher-order modes, especially at the face and rear side of the trapezoidal bar. Further investigation of this behavior, including optimization of the numerical scheme to reduce possible instabilities, with application of the present model to more complex situations and detailed analysis, including comparison against experimental measurements, will be the subject of future work.

Figure 10. The same as in Figure 9, but in the presence of opposing shear current with data as shown in Figure 8.

Moreover, in this case, the disturbance of the flow due to combined effect of variable bathymetry and vorticity contained in the vertically sheared current is seen to generate a more complicated pattern concerning the higher-order modes, especially at the face and rear side of the trapezoidal bar. Further investigation of this behavior, including optimization of the numerical scheme to reduce possible instabilities, with application of the present model to more complex situations and detailed analysis, including comparison against experimental measurements, will be the subject of future work.

A Single-Mode Weakly Nonlinear Model

Based on the above remarks concerning the convergence of the present modal series, it is useful to consider here a simplified nonlinear model for waves without current effects, derived from Equations ( 20) by keeping only the n = 0 term in the expansions of the horizontal U velocity component and dropping the term x W W  as higher-order quantity in the momentum equation.

The single-mode, simplified non-linear model reads as follows:

      2 0 0 0 0 0 0 0 0 0 0 t t x x x t t t U U g U U A U B U C U                    (11a)   0 0 0 0, t x x U cU        (31b)
where the coefficients are given by

    2 1 1 0 0 0 1, Z Z    ,       2 1 1 1 0 0 0 0 , W Z Z    ,   0 0 0 tanh / c k h k  (12a)       2 3 1 1 0 0 0 , A Z Z Z    ,       2 3 1 1 0 0 0 2 , x B Z Z Z     ,       2 3 1 1 2 0 0 0 , x C Z Z Z     . ( 32b
)
Figure 10. The same as in Figure 9, but in the presence of opposing shear current with data as shown in Figure 8.

A Single-Mode Weakly Nonlinear Model

Based on the above remarks concerning the convergence of the present modal series, it is useful to consider here a simplified nonlinear model for waves without current effects, derived from Equation ( 20) by keeping only the n = 0 term in the expansions of the horizontal U velocity component and dropping the term W∂ x W as higher-order quantity in the momentum equation. The single-mode, simplified non-linear model reads as follows:

∂ t U 0 + β 0 U 0 ∂ t η + α 0 (g∂ x η + U 0 ∂ x U 0 ) + A∂ 2 x (∂ t U 0 ) + B(∂ t U 0 ) + C∂ t U 0 = 0 (31a) ∂ t η + U 0 ∂ x η + ∂ x (c 0 U 0 ) = 0, (31b)
where the coefficients are given by

α 0 = 1, Z (1) 0 Z
(1) 0 -2 , β 0 = W (1) 0 , Z (1) 0 Z
(1) 0 -2 , c 0 = tanh(|k 0 |h)/|k 0 | (32a) A = -Z (3) 0 , Z (1) 0 Z
(1) 0 -2 , B = -2 ∂ x Z (3) 0 , Z (1) 0 Z
(1) 0 -2 , C = -∂ 2 x Z (3) 0 , Z (1) 0 Z
(1) 0 -2 . ( 32b
)
The linear dispersion characteristics of the above simplified model are the same as the ones of the present CMS with n = 1 terms shown in Figure 2. Numerical results in the case of monochromatic waves of period T = 2 s and a height of 2cm propagating over a trapezoidal bar are shown in Figures 11 and12, as obtained by a simplified version of the present single-mode model defined by Equations (31) and (32), using µ 0 h = 0.1π. In particular, comparison of the latter model results against the experimental data by Beji and Battjes [START_REF] Beji | Numerical simulation of wave propagation over a bar[END_REF] concerning a time series of free surface elevation measured at stations 2-7 of Figure 11 are shown in Figure 12 by using dotted lines. It is observed that the nonlinear simplified model is able to provide quite good predictions, which is evidence of the good properties of the nonlinear system developed in the present work. Further applications and additional investigation in the case of wave-current interaction over bathymetry will be presented in future work. It is observed that the nonlinear simplified model is able to provide quite good predictions, which is evidence of the good properties of the nonlinear system developed in the present work. Further applications and additional investigation in the case of wave-current interaction over bathymetry will be presented in future work.

Conclusions

In this work, a novel coupled-mode system is derived based on a velocity formulation and Euler equations, overcoming the irrotational assumption of standard models of waves propagating in variable bathymetry regions, focusing on the study of additional effects due to inhomogeneous currents. The present work is limited to collinear waves and currents, and the new model is based on a modal expansion of the horizontal component of the velocity field by using a vertical basis defined by local Sturm-Liouville eigenproblems. A corresponding expansion of the vertical velocity component is also derived fulfilling the continuity and the bottom boundary condition. Subsequently, projection of momentum equations in the local vertical basis leads to the new coupled system of equations with respect to the horizontal velocity mode amplitudes and the free-surface elevation, closed by the kinematic free-surface boundary condition. The dispersive behavior of the present system is analyzed in detail, for: (i) water waves propagating without current, (ii) water waves propagating in the presence of a depth-uniform current (thus involving no vorticity), and (iii) water waves propagating in the presence of a linearly sheared current, thus presenting constant vorticity. In every case, the convergence of the system is rapidly achieved, and good results are obtained by keeping only a few modes in the modal series, for all values of the shallowness parameter ranging from shallow to deep water conditions. First

Conclusions

In this work, a novel coupled-mode system is derived based on a velocity formulation and Euler equations, overcoming the irrotational assumption of standard models of waves propagating in variable bathymetry regions, focusing on the study of additional effects due to inhomogeneous currents. The present work is limited to collinear waves and currents, and the new model is based on a modal expansion of the horizontal component of the velocity field by using a vertical basis defined by local Sturm-Liouville eigenproblems. A corresponding expansion of the vertical velocity component is also derived fulfilling the continuity and the bottom boundary condition. Subsequently, projection of momentum equations in the local vertical basis leads to the new coupled system of equations with respect to the horizontal velocity mode amplitudes and the free-surface elevation, closed by the kinematic free-surface boundary condition.

The dispersive behavior of the present system is analyzed in detail, for: (i) water waves propagating without current, (ii) water waves propagating in the presence of a depth-uniform current (thus involving no vorticity), and (iii) water waves propagating in the presence of a linearly sheared current, thus presenting constant vorticity. In every case, the convergence of the system is rapidly achieved, and good results are obtained by keeping only a few modes in the modal series, for all values of the shallowness parameter ranging from shallow to deep water conditions. First results concerning the numerical performance of the present system are also presented, and a demonstration of the nonlinear behavior of the system is shown through the application of a simplified single-mode model to the case of waves propagating over a trapezoidal bar and comparison against experimental data. After further elaboration, the present system is expected to provide a useful tool for modeling wave-current interactions under the effects of vorticity. Finally, the analytical structure of the present coupled-mode system facilitates extensions to treat non-linear effects and applications concerning 3D wave scattering by inhomogeneous currents in coastal regions with general bottom topography.

Figure 1 .

 1 Figure 1. Waves propagating in variable bathymetry in the presence of current with velocity c U .

Figure 1 .

 1 Figure 1. Waves propagating in variable bathymetry in the presence of current with velocity U c .

  are normalized (i.e., Z

Fluids 2019, 4 , 24 Figure 2 .

 4242 Figure 2. Dispersion characteristics of the linearized Coupled-Mode System (CMS) in the case of waves only in constant depth strip, using one term (n = 0, indicated by 1), two terms (n = 0,1, indicated by 2), and three terms n = 0,1,2, indicated by 3) in the local mode expansions; (a) 0 0.1 h   

0 hF  and following waves to 0 hF

 00  , respectively.Results concerning a relatively strong current characterized by h F =0.1 for following waves are presented in Figure3, again for all values of the non-dimensional wavenumber kh , ranging from 0 to 2π, corresponding from very shallow to deep wave conditions. The same as before, two values of the parameter controlling the vertical basis   results to the analytical solution, which is shown by using thick lines, is observed as the modes included in the present series increases. The corresponding results in the case of opposing wave and current are presented in Figure4. In this case the convergence to the analytical results is slower, especially when the propagation speed becomes smaller approaching the blocking condition. However, in both cases, the present representation, keeping only the first term (n = 0) in the expansion, provides good results over a bandwidth of non-dimensional wavenumbers around the selected value of the parameter 0 h  , rendering the single-term model a good approximation.(iii) Current with constant vorticity: In the case of a linear vertical profile the coefficient mn B in Equation (26) becomes

Figure 2 .

 2 Figure 2. Dispersion characteristics of the linearized Coupled-Mode System (CMS) in the case of waves only in constant depth strip, using one term (n = 0, indicated by 1), two terms (n = 0,1, indicated by 2), and three terms n = 0,1,2, indicated by 3) in the local mode expansions; (a) µ 0 h = 0.1π, (b) µ 0 h = 0.5π.

  numbers: one corresponding to the current flow at the surface below the free surface. The latter parameter is also expressed in terms of Froude number h

Figure 3 .

 3 Figure 3. Dispersion characteristics of the linearized CMS in the case of waves and following parallel current h F =0.1 in constant depth strip, using one term (n = 0, indicated by 1), and two terms (n = 0, 1, indicated by 2) in the local mode expansions; (a) 0 0.1 h   

Figure 4 .F

 4 Figure 4. Dispersion characteristics of the linearized CMS in the case of waves and opposing parallel current h F =0.1 in constant depth strip, using one term (n = 0, indicated by 1), two terms (n = 0,1, indicated by 2), and three terms (n = 0, 1, 2, indicated by 3) in the local mode expansions; (a) 0 0.1 h   

Figure 3 .

 3 Figure 3. Dispersion characteristics of the linearized CMS in the case of waves and following parallel current F h = 0.1 in constant depth strip, using one term (n = 0, indicated by 1), and two terms (n = 0, 1, indicated by 2) in the local mode expansions; (a) µ 0 h = 0.1π, (b) µ 0 h = 0.5π.

Figure 3 .

 3 Figure 3. Dispersion characteristics of the linearized CMS in the case of waves and following parallel current h F =0.1 in constant depth strip, using one term (n = 0, indicated by 1), and two terms (n = 0, 1, indicated by 2) in the local mode expansions; (a) 0 0.1 h   

Figure 4 .F

 4 Figure 4. Dispersion characteristics of the linearized CMS in the case of waves and opposing parallel current h F =0.1 in constant depth strip, using one term (n = 0, indicated by 1), two terms (n = 0,1, indicated by 2), and three terms (n = 0, 1, 2, indicated by 3) in the local mode expansions; (a) 0 0.1 h   

Figure 4 .

 4 Figure 4. Dispersion characteristics of the linearized CMS in the case of waves and opposing parallel current F h = 0.1 in constant depth strip, using one term (n = 0, indicated by 1), two terms (n = 0, 1, indicated by 2), and three terms (n = 0, 1, 2, indicated by 3) in the local mode expansions; (a) µ 0 h = 0.1π, (b) µ 0 h = 0.5π.

Fluids 2019, 4 , 24 Figure 5 .

 4245 Figure 5. Dispersion characteristics of the linearized CMS in the case of waves and following shear current with constant vorticity in constant depth strip, using one term (n = 0, indicated by 1), and two terms (n = 0,1, indicated by 2) in the local mode expansions; (a,c) 0 0.1 h   

Figure 5 .

 5 Figure 5. Dispersion characteristics of the linearized CMS in the case of waves and following shear current with constant vorticity in constant depth strip, using one term (n = 0, indicated by 1), and two terms (n = 0,1, indicated by 2) in the local mode expansions; (a,c) µ 0 h = 0.1π, (b,d): µ 0 h = 0.5π. The specific values of F h and S c are included in the subplots.

Figure 6 .

 6 Figure 6. Similar to Figure 5, but for opposing waves and sheared current.

.

  defined by subdividing the horizontal extent of the domain ranging in Based on the above analysis presented in Section 4 concerning the dispersion characteristics of the present model, in conjunction with extensive experience from similar applications of the CMS, in most cases the number of the retained modes required for convergence is 4-5, and the horizontal subdivision of the domain is based on using 30-40 points per wavelength. In the examples that will

Figure 6 .

 6 Figure 6. Similar to Figure 5, but for opposing waves and sheared current.

Fluids 2019, 4 ,

 4 x FOR PEER REVIEW 16 of 24

Figure 7 .

 7 Figure 7. Results obtained by the present linearized CMS for harmonic incident waves of period T = 2 s and a height of 5 cm propagating over a linear upslope characterized by bottom slope 4.5%. The first 4 subplots present the calculated modes

Figure 7 .

 7 Figure 7. Results obtained by the present linearized CMS for harmonic incident waves of period T = 2 s and a height of 5 cm propagating over a linear upslope characterized by bottom slope 4.5%. The first 4 subplots present the calculated modes U n (x;t), n = 0, 1, 2, 3 at a given time instant. The vertical dashed lines indicate the absorbing (relaxation) zones. The last three subplots show (using thick lines) the velocity field U(x,z = 0;t) on the mean free surface (z = 0) and the free surface elevation η(x;t) at the same instant, and finally the sloping bottom topography.

Figure 7 .

 7 Figure 7. Results obtained by the present linearized CMS for harmonic incident waves of period T = 2 s and a height of 5 cm propagating over a linear upslope characterized by bottom slope 4.5%. The first 4 subplots present the calculated modes

Figure 8 .

 8 Figure 8. Bathymetry and current data for waves over the trapezoidal bar presented and discussed in Belibassakis et al. [17]. In the lower subplot, the surface current

Figure 8 .

 8 Figure 8. Bathymetry and current data for waves over the trapezoidal bar presented and discussed in Belibassakis et al.[START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF]. In the lower subplot, the surface current U 0 (x) is plotted by using a solid line and the shear S(x) = Ω(x) is shown by using a dashed line, respectively.

Fluids 2019, 4 , 24 Figure 9 .

 4249 Figure 9. Results obtained by the present CMS for harmonic incident waves of period T = 1 s propagating over the trapezoidal bar of Figure 8 without current. The last three subplots show (using thick lines) the free-surface velocity field

  n  at a given time instant are plotted along with the surface wave velocity field   , 0; U x z t  and the free surface elevation   ; x t  at the same instant, as obtained by the present linearized CMS given by Equation (

Figure 9 .

 9 Figure 9. Results obtained by the present CMS for harmonic incident waves of period T = 1 s propagating over the trapezoidal bar of Figure8without current. The last three subplots show (using thick lines) the free-surface velocity field U(x,z = 0;t) and the free surface elevation η(x;t) at the same instant, and finally the sloping bottom topography. The time-harmonic solution[START_REF] Belibassakis | A coupled-mode model for water wave scattering by vertically sheared currents in variable bathymetry regions[END_REF] is indicated by using dashed lines.

Figure 11 ..

 11 Figure 11. Horizontal and vertical flow components in the case of monochromatic incident waves of period T = 2 s and a height of 2 cm propagating over a trapezoidal bar, as calculated by the present nonlinear simplified model. The bathymetry is shown by using a thick solid line. The linear dispersion characteristics of the above simplified model are the same as the ones of the present CMS with n = 1 terms shown in Figure 2. Numerical results in the case of monochromatic waves of period T = 2 s and a height of 2cm propagating over a trapezoidal bar are shown in Figures 11 and 12, as obtained by a simplified version of the present single-mode model defined by Equations (31) and (32), using 0 0.1 h   

Figure 11 .

 11 Figure 11. Horizontal and vertical flow components in the case of monochromatic incident waves of period T = 2 s and a height of 2 cm propagating over a trapezoidal bar, as calculated by the present nonlinear simplified model. The bathymetry is shown by using a thick solid line.

Fluids 2019, 4 , 24 Figure 12 .

 42412 Figure 12. Comparison of present simplified nonlinear model results in the case of monochromatic waves of period T = 2 s propagating over the trapezoidal bar of Figure 9 against the experimental data by Beji and Battjes [27], shown by using dotted lines. Time series of free surface elevation measured at stations 2-7 of Figure 9.

Figure 12 .

 12 Figure 12. Comparison of present simplified nonlinear model results in the case of monochromatic waves of period T = 2 s propagating over the trapezoidal bar of Figure 9 against the experimental data by Beji and Battjes [27], shown by using dotted lines. Time series of free surface elevation measured at stations 2-7 of Figure 9.

Figure A1 ..

 A1 Figure A1. Modal representation of potential flow      , co sh sin U x z z h x         

Figure A1 .

 A1 Figure A1. Modal representation of potential flow z) = κ cosh[κ(z + h)] sin(κx) at x = 0 between horizontal bottom z/h = -1 and an upper surface η(x) = εh sin(κx) with ε = 0.7 and κh = 1. (a) Comparison of the present modal series defined by Equation (1) shown by using dashed lines against the potential velocity field U(x = 0, z), shown by using solid lines. (b) The calculated rate of decay of the mode amplitudes |U n | is shown by using solid lines and the improvement obtained by including additional free-surface mode by using dashed lines. The asymptotic trends are indicated in the figure by using blue dotted lines.

Fluids 2019, 4 , 24 Figure A2 .,

 424A2 Figure A2. Modal representation of Couette flow between parallel planes at / 1 z h   and / 0 z h  . (a) Comparison of the present modal series defined by Equation (1), shown by using dashed lines, against the shear flow     0 1 / U z U z h   , shown by using solid lines. (b)

Figure A2 .

 A2 Figure A2. Modal representation of Couette flow between parallel planes at z/h = -1 and z/h = 0. (a) Comparison of the present modal series defined by Equation (1), shown by using dashed lines, against the shear flow U(z) = U 0 (1 + z/h), shown by using solid lines. (b) Calculated rate of decay of the mode amplitudes |U n |. The asymptotic trend is indicated in the figure by using blue dotted lines.

  ....

	n Z z h x		x t n		, for
	each instant and at each horizontal x-position.

Acknowledgments:

The support of Kostas Belibassakis to visit Université de Toulon is greatly acknowledged. The invitation of Julien Touboul by the National Technical University of Athens, partially supported by the French Direction Générale de l'Armement (DGA), is also deeply acknowledged.

Author Contributions:

The main ideas of this work, as well as the first draft of the text, belong to K.B. and J.T.

Funding: This research was funded by the French Direction Générale de l'Armement (DGA), partially supporting the venue of Julien Touboul at the National Technical University of Athens.

Conflicts of Interest:

The authors declare no conflict of interest.

Appendix A

Using the analytical form of the vertical functions Z

(2) n (z; h, η), as defined by Equation (8), we obtain by direct differentiation with respect to the horizontal x-coordinate

From the above equation we obtain at the seabed, z = -h, the following result is obtained:

n (z = -h; h, η), n = 0, 1, 2, . . . . (A2)

Appendix B

Assuming that the field U(x, z; t) is bounded and its derivatives are integrable in the domain, and taking into account the orthogonality of the system Z

n (z; x), n = 0, 1, 2, . . . , for each horizontal x-position, we have

n (z; x)U(x, z; t)dz = -

where

x) is used. Moreover, and using the fact that ∂ z Z

(1)

and using standard asymptotic results for the regular Sturm-Liouville problem, k n ≈ nπ/(η + h), for n → ∞ , and the fact that Z

The above result ensures that the horizontal velocity field U(x, z; t), as represented by Equation (1) U(x, z; t) ∈ H 1 (D(t)).

In order to illustrate the above results, we consider below two examples. The first concerns an irrotational wavelike flow, which is also considered by Nicholls and Reitich [START_REF] Nicholls | Stability of High-Order Perturbative Methods for the Computation of Dirichlet-Neumann Operators[END_REF], and Athanassoulis and Papoutselis [START_REF] Athanassoulis | Exact semi-separation of variables in waveguides with non-planar boundaries[END_REF], and for which analytical solutions are known. This flow is defined by Φ(x, z) = cosh[κ(z + h)] cos(κx) in the domain 0 < κx < 2π, -h < z < η(x) for sinusoidal shapes of the upper boundary upper η(x) = εh sin(κx) and the corresponding horizontal flow component is U(x, z) = κ cosh[κ(z + h)] sin(κx). The result of the representation given by Equation (1), using as control frequency-type parameter µ 0 = κ tan[κh], is plotted in Figure A1 by using dashed lines and is compared against the given field shown by using solid lines. In this case the differences are mainly observed near the top boundary, as observed in the detail of Figure A1a.

The rate of convergence of the horizontal velocity series can be improved by including additional modes increasing the rate of decay of the mode amplitudes |U n |, as shown in Figure A1b by using dashed lines.

As a second example, we consider a Couette flow characterized by a linear distribution U(z) = U 0 (1 + z/h) between two parallel planes at z/h = -1 and z/h = 0. The result of the representation Equation (1) using the same value as before for the frequency-type parameter µ 0 controlling the vertical basis is plotted in Figure A2. In this case, the differences are mainly observed near the top and the bottom boundaries, as observed in the details of Figure A2a, and are due to the incompatibility between the boundary conditions of the field U(z) at z/h = -1 and z/h = 0, and the corresponding boundary conditions satisfied by the Sturm-Liouville problem, generating the vertical basis Z