Analytical Solution for Bichromatic Waves on Linearly Varying Currents

: It is well-known that the bound long waves play an important role on beach morphology. The existence of the current will modify the intensity of the waves. To examine the more realistic problem where the current is non-uniform, a third-order analytical solution for bichromatic waves on currents with constant vorticity is derived using a perturbation method. Earlier derivations of theories for interactions between waves and currents have been limited to cases of bichromatic waves on depth-uniform currents and cases of monochromatic waves on linear shear currents. Unlike the derivations of monochromatic waves, the moving-frame method cannot be used in the case of bichromatic waves because there are multiple waves with different celerities, and the ﬂow can in no way be treated as steady-state. Moreover, for shear currents where the ﬂows are rotational, the velocity potential cannot be simply deﬁned. For these reasons, it is difﬁcult to carry out mathematical operations when the dynamic free surface boundary condition is applied. However, with the consideration of the wave part of the ﬂuid motion remaining irrotational, some of the terms in the expanded boundary conditions can be ignored; thus, the derivations can be further processed. As a result, the third-order explicit expressions for the stream function, the velocity potential, and the surface elevation can be solved. The nonlinear dispersion relation is also derived to account for interacting wave components with different frequencies and amplitudes which can be used to solved for wavenumbers of both wave trains. The obtained solutions are veriﬁed by reducing to those of previous results for monochromatic waves and uniform currents. Furthermore, the inﬂuence of constant vorticity on the wave kinematics is illustrated. A comparison between following/opposing currents is also carried out. Finally, the effects of the shear current on the strength of the bound long waves are also illustrated.


Introduction
Wave-current interaction is a common phenomenon in the ocean. In ocean engineering, it is usually taken into consideration when dealing with the problems about estuaries, rip currents, and longshore currents. In nature, the currents are typically non-uniform, and the waves always consist of multiple wave components, which results in a very complicated flow field. One of the most important phenomena associated with multichromatic waves is the release of bound long waves, which may cause a change in the morphology of the inter-tidal zone and may even impact the stability of beaches [1][2][3]. However, the effects of the vorticity of currents on bound waves have not been widely discussed. Therefore, it is essential to establish theories for related circumstances in order to gain insights into these effects.
There have been several decades during which a large number of studies on wavecurrent interactions have been presented. In the case of depth-uniform currents, a simple concept of the first-order solutions for monochromatic waves was introduced by [4

]. Years
Mathematics 2022, 10, 1666 2 of 21 later, Ref. [5] derived third-order solutions and further described the effects of currents on wave characteristics such as the wavelength and wave height of monochromatic waves. In studies on monochromatic waves and shear currents, Ref. [6] provided insight into concepts used to determine the effects of vertically sheared current on waves. Ref. [7] also provided third-order analytical solutions for monochromatic waves on linearly varying currents and further showed the influence of shear currents on wave characteristics. Solutions for multichromatic waves were first introduced by [8], who gave the general second-order analytical solutions for the superposition of an infinite number of waves. Ref. [9] reduced the wave conditions to bichromatic waves and took uniform currents into consideration. They presented the third-order theory, which was based on three-dimensional fluid motion.
The objective of this study is to provide a fundamental understanding of wave-current interaction by establishing the third-order theory for bichromatic waves on vertically linear shear currents. In the first section of this paper, the mathematical formulation is given. Perturbation technique is used to separate various degree of nonlinearities up to third order. Next, the analytical solutions for wave kinematics are derived and discussed. The derived solutions are verified with the theories from previous research, the solutions of which are the subsets of this study. Additionally, an illustration of the solution is given. Furthermore, the nonlinear dispersion relation is discussed and the effects of vorticity on wave-wave interaction and wave characteristics are examined. Finally, the findings are summarized

Boundary Value Problem
We consider an incompressible, inviscid fluid where the flow is rotational because of the existence of a vertically non-uniform current. The non-uniform current is a linear shear current, where the vorticity is a constant value of Ω 0 . In this study, the flow is considered to be two-dimensional so that the waves and the current are collinear as a first step. The fluid domain is bounded by a horizontal, solid, impermeable bed and a free surface. A sketch of the coordinate system is given in Figure 1. The velocity of the current at the surface, the mid-depth, and the bottom are U s , U m , and U b , respectively. The velocity profile of the current is defined by U c (z) = U s + Ω 0 z = U m + Ω 0 z + 1 2 h = U b + Ω 0 (z + h) (1) In studies on monochromatic waves and shear currents, Ref. [6] provided insight concepts used to determine the effects of vertically sheared current on waves. Ref. [7] provided third-order analytical solutions for monochromatic waves on linearly vary currents and further showed the influence of shear currents on wave characteris Solutions for multichromatic waves were first introduced by [8], who gave the gen second-order analytical solutions for the superposition of an infinite number of wa Ref. [9] reduced the wave conditions to bichromatic waves and took uniform currents consideration. They presented the third-order theory, which was based on th dimensional fluid motion.
The objective of this study is to provide a fundamental understanding of wa current interaction by establishing the third-order theory for bichromatic waves vertically linear shear currents. In the first section of this paper, the mathemat formulation is given. Perturbation technique is used to separate various degree nonlinearities up to third order. Next, the analytical solutions for wave kinematics derived and discussed. The derived solutions are verified with the theories from previ research, the solutions of which are the subsets of this study. Additionally, an illustra of the solution is given. Furthermore, the nonlinear dispersion relation is discussed the effects of vorticity on wave-wave interaction and wave characteristics are examin Finally, the findings are summarized

Boundary Value Problem
We consider an incompressible, inviscid fluid where the flow is rotational becaus the existence of a vertically non-uniform current. The non-uniform current is a linear sh current, where the vorticity is a constant value of . In this study, the flow is conside to be two-dimensional so that the waves and the current are collinear as a first step. fluid domain is bounded by a horizontal, solid, impermeable bed and a free surface sketch of the coordinate system is given in Figure 1. The velocity of the current at surface, the mid-depth, and the bottom are , , and , respectively. The velo profile of the current is defined by Since the fluid is assumed to be incompressible, the fluid motion can be descri with the stream function . Moreover, when the vorticity of the shear current is const the fluid motion can be decomposed using Helmholtz's decomposition [6] into rotatio Since the fluid is assumed to be incompressible, the fluid motion can be described with the stream function ψ. Moreover, when the vorticity of the shear current is constant, the fluid motion can be decomposed using Helmholtz's decomposition [6] into rotational and irrotational fields which represent the current part and the wave part of the motion. Therefore, the wave part of the motion that is irrotational can be described by defining the velocity potential φ and can be written as (2) where u and w are components of the particle velocity in the x and z directions. The governing equations for the problem are The bottom boundary condition is The kinematic free-surface boundary condition is ∂η ∂t The dynamic free-surface boundary condition is ∂φ ∂t where R(t) is the Bernoulli time-dependent term yielded from integration. It should be noted that, typically, it is convenient to omit the time-derivative term using a moving frame to make the waves steady-state. However, in the case of bichromatic waves, owing to existence of two wave celerities, using the same method to obtain steadystate wave field cannot work. The free surface boundary conditions are established on the unknown position η. In order to solve the boundary value problem, they are expanded by a Taylor series at the position z = 0.

Perturbation Method
The stream function ψ, velocity potential φ and the surface elevation η are expanded as perturbation series with the steepness being the parameter for the expansion. Assume that the parameter (steepness) naturally appears in the solutions.
where the subscripts n and m denote the two wave trains, respectively. Note that the zeroth-order velocity potential does not exist because the non-uniform current motion, which is indicated by the zeroth-order problem, is rotational.

Zeroth-Order Problem
The zeroth-order solution is the solution of the steady linear shear current itself with the wave height that is zero. Collecting zeroth-order terms from the governing equation and the expanded boundary conditions yields the following: Governing Equations: BBC: KFSBC: DFSBC: The zeroth-order solutions are Note that Equation (16) is actually the stream function of the linear shear current as defined.

First-Order Problem
The perturbation equations to the first order are Governing Equations: KFSBC: DFSBC: By substituting the zeroth-order solutions into Equations (18)-(21) to solve for the stream function ψ and the velocity potential φ and applying the method of separation of variables, the angular frequency σ is then defined to ensure the wave motion being periodic in time, and the wavenumber k is defined to ensure it being periodic in space.
Therefore, the first-order solutions can be readily obtained as    ψ 1n = Γ 1n sinh[k n (z + h)](a n cos Θ n ) φ 1n = Γ 1n cosh[k n (z + h)](a n sin Θ n ) η 1n ≡ a n cos Θ n where Θ n = k n x − σ n t (23) The angular frequency can be obtained by combining KFSBC and DFSBC as follows: Equation (26) is a quadratic equation of (σ n − U s k n ), which can further be expressed as follows: The plus sign is chosen instead of the minus sign since (σ n − U s k n ) is the modified angular frequency under the Doppler effect and should be positive. Therefore, Equation (27) can be further written as with same expressions for σ m . For later use, the zeroth-order angular frequency is expressed by where σ 0n = gk n tanhk n h σ 0n +Ω 0 tanhk n h = gk n tanhk n h − σ 0n Ω 0 tanhk n h = 1 2 −Ω 0 tanhk n h + Ω 2 0 tanh 2 k n h + 4gk n tanhk n h all with equivalent form for "m" wave train. The dispersion relations of the two wave trains on the linear shear current are thus obtained. The first order Bernoulli term solved with DFSBC is obtained as

Third-Order Problem
Secular behavior is found while solving the third-order problem. In order to eliminate the secular terms, the angular frequencies of both wave trains also have to be expanded with the perturbation method as To eliminate the secular terms, we can use the method proposed by [10], where we let Now, we can substitute all the time-derivative terms with τ-derivative terms in the third-order combined free surface boundary condition, where the derivative terms have the relationship σ 2n and σ 2m are thus solved as The transfer functions Θ and Π are defined as The third-order problem is as follows: Governing Equations: KFSBC: The third-order solutions can be expressed as: Again, "3m" component has the equivalent form to that of "3n". The coefficients for the self-self-self-interaction components are with the equivalent expressions for Γ 3m and Λ 3m . The other coefficients can be expressed by transfer functions as where Finally, the coefficients for the third-order correction to the first-order surface elevation (necessary to eliminate the secular terms) are where Λ 13 (a n , a m , k n , k m , σ 0n , σ 0m , Γ 2n , Γ ± nm , Λ 2n , Λ ± nm ) = 1 gk n σ 2 0n σ 2n coth k n h + 1 2g Γ 2n a 2 n k n σ 0n [sinh 2k n h − coth k n h cosh 2k n h] 8g a 2 n k n σ 2 0n coth k n h − 3 8g a 2 n k n σ 0n Ω 0 + 1 4g a 2 m k n σ 2 0n coth k n h −2k m σ 2 0m coth k m h − 2k n σ 0n σ 0m coth k m h − 2k m σ 0n σ 0m coth k n h − 1 4g a 2 m Ω 0 [k n σ 0n + 2k m σ 0m ] The third-order Bernoulli term is derived as The third-order theory for the effects of bichromatic waves on linear shear currents is thus completed. These analytical solutions provide explicit expressions for the velocity potential, the surface elevation, the frequency dispersion and the amplitude dispersion for bichromatic waves on linear shear currents at a finite depth.
The solutions are verified by reducing to those for bichromatic waves on uniform currents and those for monochromatic waves on linear shear currents. By setting a m . to zero, only the solutions of 1n, 2n, and 3n remain. The solutions thus become those for monochromatic waves on linear shear currents, and the reduced solutions are exactly the same as those obtained [7].
By setting Ω 0 to zero, the solutions are reduced to those for bichromatic waves on uniform currents. The solutions can be verified with those by [9]. The coefficients in stream function and velocity potential thus become The coefficients are exactly the same as those found by [9]. With the solutions having the same form, they are thus successfully verified.

Velocity of Water Particles and Pressure Field
By combining the solutions of each order, the stream function up to the third-order is expressed by The velocity potential to describe the wave motion is expressed by The surface elevation is expressed by η= a n cos θ n + a m cos θ m + a 2 n 2 Λ 2n cos 2θ n + a 2 m 2 Λ 2m cos 2θ m +a n a m Λ + nm cos(Θ n + Θ m ) + a n a m Λ − nm cos(Θ n − Θ m ) + a 3 n 2 Λ 3n cos 3θ n + a 3 n 2 Λ 3m cos 3θ m + a n a 2 m 2 Λ + n2m cos(Θ n + 2θ m ) + a n a 2 Mathematics 2022, 10, 1666 13 of 21 By substituting Equation (82) into Equation (2), the horizontal and vertical velocity of water particles are obtained as The Bernoulli constant up to the third order is The pressure field up to the third order can be obtained by substituting Equation (82) to Equation (85) into the Bernoulli equation, and the pressure field is obtained as with the result above, the essential physical properties of bichromatic waves on linearly varying current up to the third order are thus completely obtained.

Example of Bichromatic Wave on Vertically Linear Shear Current
The solutions are illustrated by considering the following case: h = 5 m, a n = 0.45 m, σ n = 2π × 0.4 s −1 , a m = 0.20 m, σ m = 2π × 0.3 s −1 , U m = 0.35 m/s, Ω 0 = 0.10 s −1 . Figure 2 shows the surface elevation at t = 0 and the horizontal velocity profile at t = 0, x = 0, where the solid lines indicate the third-order solutions, and the dashed lines indicate the first-order solutions. The relevant coefficients corresponding to the example are given in Tables 1 and 2. The difference between the first-and third-order solutions is significant in terms of the surface elevation because of the nonlinear dispersion relation, where the wavelengths are modified by the second-order frequency correction. It is also seen that in the horizontal velocity profile, the second-and third-order components of the horizontal particle velocity weaken the maximum velocity compared to the first-order solution.
Mathematics 2022, 10,1666 15 of 22 the horizontal velocity profile, the second-and third-order components of the horizontal particle velocity weaken the maximum velocity compared to the first-order solution.  In cases where waves propagate in the opposite directions, it is seen that the solutions

Singularity of Solving for Wavenumbers with the Nonlinear Dispersion Relation
In cases where waves propagate in the opposite directions, it is seen that the solutions for the wavenumbers cannot be obtained. In order to obtain wavenumbers, the dispersion relation must be satisfied. The dispersion relation can be rewritten as The solutions for the wavenumbers can be found at where F 1 = 0, and F 2 = 0. For example, in Figures 3 and 4, for the case where Ω 0 = −0.1 s −1 , solutions can be found at U m = −0.5 (m/s) but cannot be found at U m = 0.5 (m/s). Figure 3 shows when U m = −0.5 (m/s), the point of intersection (where the solutions can be found) can be found at k A = −0.4790 m −1 and k B = 0.5029 m −1 , which is a group of solutions for wave A propagating in the -x direction and wave B propagating in the +x direction. When U m = 0.5 (m/s), no solution can be found in the region where k A < 0 and k B > 0. It is noted that, even though in Figure 4, some points of intersection can be found, no solutions can be obtained because either F 1 or F 2 is singular near the points of intersection. Figures 5 and 6 show the region where −2 ≤ F 1 ≤ 2 and −2 ≤ F 2 ≤ 2 when U m = −0.5 (m/s), respectively, the same as in Figures 7 and 8 when U m = 0.5 (m/s). In Figures 5 and 7, the yellow region indicates where −2 ≤ F 1 ≤ 2, which means that no solutions can be found in the white area. Singularity lies near the point where the lines are yellow on one side and white on the other. Therefore, the existence of the points of intersection does not necessarily mean the existence of solutions. For example, even though a point of intersection can be found near k A = −0.40 m −1 and k B = 0.35 m −1 in Figure 4, no solutions can be found because F 1 is near a singularity at that position (see Figure 7). Likewise, in Figure 3, a point of intersection can be found near k A = −0.47 m −1 and k B = 0.11 m −1 , but no solutions for wavenumbers can be found there because F 2 is near a singularity at that point (see Figure 6). athematics 2022, 10, 1666 16 The solutions for the wavenumbers can be found at where = 0, and = 0. For example, in Figures 3 and 4, for the case where = −0.1(s ), solutions ca found at = −0.5 (m/s) but cannot be found at = 0.5 (m/s). Figure 3 shows w = −0.5 (m/s), the point of intersection (where the solutions can be found) can be fo at = −0.4790 (m ) and = 0.5029 (m ), which is a group of solutions for wav propagating in the -x direction and wave B propagating in the +x direction. When 0.5 (m/s), no solution can be found in the region where < 0 and > 0. It is noted even though in Figure 4, some points of intersection can be found, no solutions ca obtained because either or is singular near the points of intersection. Figures 5 and 6 show the region where −2 ≤ ≤ 2 and −2 ≤ ≤ 2 when −0.5 (m/s), respectively, the same as in Figures 7 and 8 Figure 4 solutions can be found because is near a singularity at that position (see Figur Likewise, in Figure 3, a point of intersection can be found near = −0.47 (m ) and 0.11 (m ), but no solutions for wavenumbers can be found there because is ne singularity at that point (see Figure 6).

Influence of the Shear Current on Wavelength
The wavelength can be determined as the solution of the dispersion relation w the wave frequency, the wave amplitude, the water depth, and the current velocity pr are known. For the following results and discussions, the wave properties of wave A B are as follows: (2 = ). In Figures 9 and 10, the wavelength is plotted as a function of vorticity , w different current velocities at the mid-depth for both waves A and B. It can be seen the wavelength increases as the vorticity increases, whether it applies to a follow current or an opposing current. It is also noted that the influence of the vorticity is m significant in the first-order solution than in the third-order solution for the cas opposing current. For wave A, whose kh is larger, the increase is more significant than wave B.

Influence of the Shear Current on Wavelength
The wavelength can be determined as the solution of the dispersion relation where the wave frequency, the wave amplitude, the water depth, and the current velocity profile are known. For the following results and discussions, the wave properties of wave A and B are as follows: a A = 0.45 m, a B = 0.20 m, f A = 0.4 s −1 , f B = 0.3 s −1 , and h = 5.0 m (2Πf = σ). In Figures 9 and 10, the wavelength is plotted as a function of vorticity Ω 0 , with different current velocities at the mid-depth for both waves A and B. It can be seen that the wavelength increases as the vorticity increases, whether it applies to a following current or an opposing current. It is also noted that the influence of the vorticity is more significant in the first-order solution than in the third-order solution for the case of opposing current. For wave A, whose kh is larger, the increase is more significant than for wave B. different current velocities at the mid-depth for both waves A and B. It can be seen that the wavelength increases as the vorticity increases, whether it applies to a following current or an opposing current. It is also noted that the influence of the vorticity is more significant in the first-order solution than in the third-order solution for the case of opposing current. For wave A, whose kh is larger, the increase is more significant than for wave B.   current or an opposing current. It is also noted that the influence of the vorticity is more significant in the first-order solution than in the third-order solution for the case of opposing current. For wave A, whose kh is larger, the increase is more significant than for wave B.

Influence of Shear Current Velocities and Vorticities on Bound Long Waves
For the case of a bichromatic wave group consisting of waves A and B, the frequencies of the bound long waves nm − and n2m − are 0.1 and 0.2 Hz, respectively. Intensity can be described by the coefficients Γ and Λ, which indicate the intensity of the particle velocity and wave height, respectively.
The coefficient Γ and Λ are plotted as functions of current velocity in Figures 11 and 12. In Figure 12, the absolute value of Γ − n2m decreases as the current velocity increases, whereas the absolute value of Λ − n2m also decreases as the current velocity increases. Note that the intensity depends on the absolute value of the coefficient because of the symmetry with respect to the mean water level.
In Figures 11 and 12, it can be seen that for the second-order sub-harmonic bound long wave (nm − ), the intensity increases as the current velocity increases, whereas for the third-order sub-harmonic bound long wave (n2m − ), the opposite occurs.

Influence of Shear Current Velocities and Vorticities on Bound Long Waves
For the case of a bichromatic wave group consisting of waves A and B, the frequencies of the bound long waves and 2 are 0.1 and 0.2 Hz, respectively. Intensity can be described by the coefficients and , which indicate the intensity of the particle velocity and wave height, respectively.
The coefficient and are plotted as functions of current velocity in Figures 11 and  12. In Figure 12, the absolute value of decreases as the current velocity increases, whereas the absolute value of also decreases as the current velocity increases. Note that the intensity depends on the absolute value of the coefficient because of the symmetry with respect to the mean water level.
In Figures 11 and 12, it can be seen that for the second-order sub-harmonic bound long wave ( ), the intensity increases as the current velocity increases, whereas for the third-order sub-harmonic bound long wave ( 2 ), the opposite occurs.
However, it is also noticed in the illustration of in Figure 12 that the nonmonotonicity of with vorticity is found as a stronger opposing current (i.e. − 0.5 m/s) occurs. Figure 13 shows that is a function of vorticity when The effects of the shear current vorticities on bound long waves were also examined. Particularly, the coefficients as functions of vorticity with the same wave properties are shown in Figures 15 and 16. In this case, the intensity increases as the vorticity increases regardless of whether the current is following or opposing both the second-order and the third-order bound long waves.     However, it is also noticed in the illustration of Λ − n2m in Figure 12 that the nonmonotonicity of Λ − n2m with vorticity is found as a stronger opposing current (i.e. U m − 0.5 m/s) occurs. Figure 13 shows that Λ − n2m is a function of vorticity when U m = 0.5 (m/s) is monotonic (Λ − n2m increases as the vorticity increases). However, in Figure      The effects of the shear current vorticities on bound long waves were also examined. Particularly, the coefficients as functions of vorticity with the same wave properties are shown in Figures 15 and 16. In this case, the intensity increases as the vorticity increases regardless of whether the current is following or opposing both the second-order and the third-order bound long waves.

Concluding Remarks
In this work, new third-order analytical stream function and surface elevation solutions for bichromatic waves on linearly varying currents are presented. The effects of the shear current on the wave kinematics are discussed.
Based on the present results, some concluding remarks can be made: 1.
The nonlinear dispersion relation, which is derived while solving for the third-order components, shows the amplitude dispersion (see Equation (53)). For bichromatic waves on linear shear currents, it can be seen that the amplitude dispersion of the two wave trains involves the effects of the amplitude on each other.

2.
A bichromatic wave example is illustrated, where the first-order and the third-order solutions are compared in terms of both surface elevation and their velocity profiles.
In the example, a significant difference is found in the wavelength. Moreover, it is seen that in the discussed case, the first-order solution overestimates the maximum velocity, especially near the free surface.

3.
Wavenumbers may not always be found in the framework of proposed approach. They can only be found at the points of intersection of the two dispersion relations, and if the points are near the singularity, the solution cannot be obtained, either. For the discussed case, a group of wavenumbers with opposite directions can be obtained when U m = −0.5 (m/s) and Ω 0 = −0.1 s −1 , but no solutions are found when U m = 0.5 (m/s) and Ω 0 = −0.1 s −1 .

4.
From the discussed cases, it is found that as the vorticity of the current increases, both the wavelength and the maximum horizontal particle velocity increase. It is also found that the increase is more prominent for a shorter wave component. 5.
The intensities of the bound long waves are also discussed. It is found that for the second-order sub-harmonic bound long wave (nm − ), the intensity increases as the current velocity increases, whereas for the third-order sub-harmonic bound long wave (n2m − ) the intensity decreases as the current velocity increases. When the current velocity at the mid-depth U m is between −0.35 and 0.35 m/s, for both the second-and the third-order bound long waves, the intensity is found to increase as the vorticity increases, regardless of whether the current is following or opposing. However, when a stronger opposing current occurs (i.e. U m = −0.5 m/s), the intensity of the thirdorder bound long wave as a function of vorticity is found to be non-monotonic. A minimum Λ − n2m is found at about Ω 0 = −0.3 s −1 , with a value of 0.1938 m −2 . 6.
For primary wave trains in the opposite directions, the bound long wave can still be illustrated as functions of current velocity or vorticity. In the discussed cases, as the current velocity increases, Γ + nm increases, but Λ + nm decreases for all of the discussed vorticities. However, Γ + n2m and Λ + n2m are shown to be non-monotonic as the current velocity increases. Furthermore, for both positive and negative current velocities, as the vorticity increases, Γ + nm increases, but Λ + nm decreases. For the positive current velocity case, as the vorticity increases, Γ + n2m increases, and Λ + n2m decreases, and for the negative current velocity case, the result is the opposite. 7.
The derived solutions have only been verified with the theories from previous researches which are the special cases of this study. Therefore, it is suggested that a set of numerical simulations or laboratory experiments of bichromatic waves on shear currents are desired to further compare the findings in the present study.