The Fourth-Order Nonlinear Schrödinger Equation and Stability Analysis for Stokes Waves on Slowly Varying Topography

The surface gravity wave equation is expanded to the fourth-order wave steepness on slowly varying topography, obtaining a topographic modified nonlinear Schrödinger (TMNLS) equation. When the time scale is longer than ε -3 times of the dominant wave period or the space scale is larger than ε -3 times the dominant wavelength, the second water depth derivative and the square of the first water depth derivative affect the first-order wave amplitude. The instability area for a uniform Stokes wave train by small perturbations is the entire wavenumber space, except for a specific stability curve on infinite and slowly varying depth. The depth variation terms affect the growth rate of uniform Stokes wave train on the order of 0.01. The stability curve shows more sensitive to the depth variation in x direction than that in y direction. The increment of the value for depth variation in x direction contributes the stable wave number of perturbation to approach or parallel to y axis. The increment of the value for depth variation in y direction helps the stable wave number of perturbation to approach or parallel to x axis.


INTRODUCTION
The interactions among wave packets with narrowband range frequencies and wavelengths received considerable attention. Benjamin and Feir (1967) theoretically proved that wave packets were unstable when kh is larger than 1.363, where k is the dominant wavenumber, and h is the water depth. Whitham (1967) identified and explained the Benjamin-Feir instability by theoretical analysis. Benney and Roskes (1969) complemented Whitham's theory. With a pair of nonlinear conservation equations introduced by Whitham, Lighthill (1967) analyzed the nonlinear wave evolution process after the initial stable stage of a single wave packet. Chu and Mei (1971) added the modulation rate term to Whitham's equation for long-term wave packet evolution processes to interpret the high-order dispersion effect.
The wave packet evolution process can be studied by the cubic nonlinear Schrödinger (NLS) equation [Zakharov (1968); Benney and Roskes (1969)], which is equal to the conservation equation proposed by Chu and Mei (Hasimoto and Ono, 1972). Zakharov and Shabat (Davey, 1972) solved the analytical NLS solution. Hasimoto and Ono (Zakharov and Shabat, 1972) obtained a onedimensional NLS for the wave packet envelope from multiple scale expansions for finite depth. Davey and Stewartson (1974) investigated the transformation of slowly varying wave packets in a three-dimensional finite depth, concluding that the wave packet envelope is confined to two nonlinear partial differential equations, similar to NLS in form. Lo and Mei (1985) highlighted the simulation value with NLS's preferably approached measured value if e < 0.1, while the departure between the simulation and measured values was larger when e > 0.1. Martin (Martin H. Yuen, 1980) found that the simulated wave energy with NLS does not satisfy the conservation law due to energy attenuation. Besides the confined condition of small wave steepness, NLS showed Benjamin-Feir instability in an unbounded region by twodimensional sideband perturbation, resulting in leaking energy from low wavenumber components to high wavenumber components (Dysthe, 1979).
The NLS has been modified to overcome these defects above. By expanding the equation to the fourth order on finite depth, Dysthe (1979) established the modified nonlinear Schrödinger equation (MNLS) to improve wave packets' instability properties. Lo and Mei (1987) numerically solved and transformed MNLS in moving coordinates, showing poor longtime wave packet evolution reproducibility on infinite depth and asymmetry of sideband perturbation evolution. By considering that the wave spectrum in a realistic ocean is not narrowband, Trulsen and Dysthe (1996) developed a broader band modified nonlinear Schrödinger equation (BMNLS) by extending the bandwidth to the e 1 2 order. It showed the same precision as MNLS in nonlinear terms but higher precision in linear dispersive terms than MNLS. In deep water, to compute the dispersive relation efficiently using the pseudo-spectrum method and keep the simple structure of the Dysthe equation, Trulsen et al. (2000) used cubic nonlinear terms to modify linear dispersive relations, improving wave packets' instability property by comparing the result with that from the Stokes wave analytic solution of Mclean (McLean, 1982), ensuring the boundness of Benjamin-Feir instability and stopping the wave energy leaking to high wavenumber components. Craig et al. (2012) proposed that NLS is not a Hamilton partial differential equation but an approximation to the Euler equation. Craig et al. (Craig et al., 2010;Craig et al., 2011) adopted Hamilton's method to solve the nonlinear wave modulation process and provided a Hamilton structure of the Dysthe equation (Dysthe, 1979) to describe the gravity wave evolution process on finite and infinite depth. Craig et al. (2012) developed the Hamilton method by introducing the Hamilton pair to the equations proposed by Trulsen et al. (Trulsen and Dysthe, 1996;Trulsen et al., 2000). The equations developed using this method were compatible with water wave equations. Zhang and Li (2012) modified the pseudo-spectrum method by splitting the technique to make the MNLS equation suitable for nonperiodic boundary conditions. The method can efficiently solve nonlinear wave equations through numerical examples by nonlinear parabolic and MNLS equations.
The topography under nature waters is complicated and a crucial factor affecting the propagation process of surface gravity waves in coastal areas. Mei (2005) solved the weakly nonlinear narrowband wave packet equation on a finite flat bottom, indicating the effect of a three-order nonlinearity on the first-order wave height when the time scale is longer than e -2 of the dominant wave period or space scale is larger than e -2 of the dominant wave length. Brinch-Nielsen and Jonsson (1986) extended the nonlinear Schrödinger equation to the fourth order in three dimensions on an arbitrary constant depth, concluding that water depth can affect the applicability of wave instability expressions in deep water. Mild slope equations (Berkhoff, 1972;Lozano and Meyer, 1976) and their extensions (Kirby, 1986;Chamberlain and Porter, 1995;Miles and Chamberlain, 1998;Agnon and Pelinovsky, 2001) are powerful tools, aiming either at steeper slopes on large length scales or shorter irregularities, primarily used for calculating wave fields on the background of ocean engineering. According to Yue and Mei (1980), restricting ∇ h h to O (e 2 ), Kirby (Kirby and Dalrymple, 1983) introduced two variable x scales and one variable y scale. A parabolic equation with time independence was developed, avoiding the caustics and irregular focusing on the ray approximation while precluding wave instability analysis (Kirby and Dalrymple, 1983). Xiao and Lo (Xiao and Lo 2004) introduced the first-order depth variation terms to NLS by expanding the equation to the third-order to allow . No stable region exists for a uniform Stokes wave on varying bottoms, and a higher order instability beyond the Benjamin-Feir type is introduced by depth variation (Xiao and Lo 2004). Combined with experimental results and numerical analysis, Li et al.  found additional wave packets propagating freely and arising at first and second orders in wave steepness in a Stokes expansion as the wave packet travels over a sudden depth transition area. Free and bound waves coexisting with different phases at the second-order wave steepness indicated that the combination of the local transient peak and the magnitude of the linear free waves explained the rogue waves observed after a sudden depth transition. Zhang and Benoit (2021) proposed that the wavebottom interaction in coastal areas forms rogue waves and increases the possibility of big waves occurring.
Neither MNLS nor BMNLS can describe the wave packet evolution process on varying bottoms. Considering the extensive application of NLS, the wave evolution process and its instability features in realistic situations can be evaluated by improving NLS to the fourth order for variable depths. Based on a mathematical technique introduced by Mei (Chu and Mei, 1970) and a boundary condition adopted by Kirby (Kirby and Dalrymple, 1983), the narrowband wave packet evolution equation was expanded to the fourth-order wave steepness. A topographic modified nonlinear Schrödinger (TMNLS) equation is obtained and investigated for the instability of a uniform Stokes wave train.

EVOLUTION EQUATION FOR NARROWBAND WAVE PACKETS ON THE FINITE DEPTH
We assign (x, y, z) as spatial coordinates with z pointing vertically upward and assume that the water depth h (x, y) is slowly varying and finite at p 10 < kh < p. In the irrotational current field of inviscid and incompressible fluids, velocity potential Ф (x, y, z, t) and free surface displacement ς (x, y, t) describe surface wave propagation. P a is the local atmospheric pressure, and the equations to describe waves are as follows.
Laplace equation: Kinematic boundary condition on the bottom: Dynamic boundary condition on the surface: Kinematic boundary condition on the surface: We act the operator ∂ ∂ t + u ! · ∇ on two sides of Equation (3). p a is a constant. Equation (3) can be given as The variable of Equation (5) is expanded into the Taylor series about (z = 0) to the fourth order, yielding Thus Equations (6) and (8) are the (ka) 4 order. k is the dominant wavenumber, and a is the dominant wave's amplitude. It is supposed that ka=e≪1.
It is supposed that the dominant wave direction is along the xaxis, and the wave packet is slowly modulated. The multiscale variables are x, x 1 = ex, x 2 = e 2 x, …::, y 1 = ey, y 2 = e 2 y, …::, t, t 1 = et, t 2 = e 2 t, …::, We expand the velocity potential and wave displacement into a perturbation series Where f n = f n x, x 1 , x 2 , …; y 1 , y 2 , …; z; t, t 1 , t 2 , … ð Þ ς n = ς n x, x 1 , x 2 , …; y 1 , y 2 , …; t, t 1 , t 2 , … ð Þ (10 À 1) The variable of Laplace Equation (1) is expanded into a perturbation series to the fourth order, yielding The variable of Equation (6) is expanded into a perturbation series to the fourth order, yielding The variable of boundary condition (8) is expanded into a perturbation series to the fourth order, yielding f n-m = (f n, m )*, ()*, and c.c. are complex conjugate numbers. According to the boundary condition introduced by Kirby (Kirby and Dalrymple, 1983), depth h is modulated at the x and y directions as The variable of the bottom boundary condition (2) is expanded into a perturbation series to the fourth order, yielding Equations (11), (12), (13), (14), and the bottom boundary condition (16) constitute definition conditions.
The free surface's leading-order displacement is A is the free surface leading-order displacement amplitude. Under third and fourth-order definition conditions, the firstorder wave height's dimensionless equation are The coefficients in Equations (17) and (18) are f 10 is the velocity potential of wave-induced current in equation (17) and (18). Equation (17) indicates that the water depth's second derivative and the square of the water depth's first derivative affect the first-order wave amplitude when the time scale is longer than e -2 times of the dominant wave period or the space scale is larger than e -2 of the dominant wavelength. Equation (18) is more complicated than Equation (17). Equation (18) improves the coefficient of Equation (17) by incorporating depth variation and wave-induced current terms. In Equation (18), higher-order dispersive and nonlinear terms are added. Equation (18) indicates that when the time scale is longer than e -3 times of the dominant wave period or the space scale is larger than e -3 times the dominant wavelength, the second water depth derivative and the square of the water depth's first derivative affect the first-order wave amplitude.
Equation (19) is the same as Equation (2.18) of reference (Kirby and Dalrymple, 1983). It cannot analyze wave instability properties because the equation is steady.
When kh limits to infinite, the coefficients of Equation (18) Omitting the terms about f 10 , except for ∂ f 10 ∂ x before A, Equation (18) In contrast with the MNLS (Trulsen and Dysthe, 1996), Equation (20) is added by topography variation terms of( − 1 The first and second-order water depth derivatives affect the first-order wave amplitude on infinite depth. Equation (20) can be called a topographic modified nonlinear Schrödinger equation (TMNLS).

INSTABILITY OF A UNIFORM STOKES WAVE TRAIN
It is supposed that the Stoke wave solution is.A = a 0 e − i 2 a 2 0 t Its instability can be evaluated by assuming small perturbations in amplitude and phase. m and l are the wavenumbers of small perturbations in x and y direction.
It is supposed that small perturbations have the plane wave solution According to Equation (20), the dispersion relation for perturbation is Thus, As shown in Equation (26), ∂ 2 h ∂ y 2 influences the perturbation phase without affecting the perturbation amplitude. Equation (27) demonstrates that ∂ h ∂ y and ∂ h ∂ x affect the small perturbation's amplitude. ImW is defined as growth rate of Stokes wave disturbed by perturbation by Lo and Mei (1987). To ensure the ungrowth of perturbations, ImW = 0, meaning A uniform Stokes wave disturbed by small perturbations is stable only when Equation (28) is satisfied. Therefore, the small perturbation's instability area is the entire perturbation wavenumber space, except for the curve satisfying Equation (28). It is shown that there are solutions for ImW = 0 when Equation (28) is satisfied. When ∂ h ∂ x is a higher order of magnitude than ∂ h ∂ y , ( ∂ h ∂ x ) 2 is neglected. Then, The first depth derivative perpendicular to the dominant wave packet imposes on the wave instability disturbed by small perturbations when ∂ h ∂ x is a higher order of magnitude than ∂ h ∂ y . Without regard to the bottom slope, ∂ h ∂ x = ∂ h ∂ y = 0, then ImW must have a real root, then The left-hand of Inequality (31) is the same as the right-hand of Equation (18) in reference (Trulsen and Dysthe, 1996) when kh limits to infinite. Inequality (31) stands for the Stokes wave instability area disturbed by the MNLS perturbation on a flat bottom. A uniform Stokes wave train is unstable disturbed by small perturbations on infinite and slowly varying depth, except for the curve satisfying Equation (28). Compared with the MNLS perturbation analysis on the flat bottom when h x = h y = 0, the small perturbation's instability area is the entire wavenumber space, except for the curve satisfying Equation (28) on a slowly varying bottom.

DISCUSSIONS
According to the results of Mclean (McLean, 1982) and Trulsen (Trulsen and Dysthe, 1996), we choose a 0 = 0.0995 and a 0 = 0.196 to plot the stability curves for |Im W| and Im W = 0, corresponding to e = 0.1 and e = 0.2. |Im W| is the growth rate of Stokes wave disturbed by perturbation and the curve of Im W = 0 is the stable curve for Stokes wave. Figures 1-12 and Figures 1-24 in the Supplementary Material show the curves of |Im W| and Im W = 0 for a 0 = 0196, the value of ∂ h ∂ x and ∂ h ∂ y ranging from 0 to 0.3. Figure 25 to Figure 60 in the Supplementary Material show the curves of |Im W| and Im W = 0 for a 0 = 0.0995, the value of ∂ h ∂ x and ∂ h ∂ y varying from 0 to 0.3. To distinguish the influence of the orders of bottom variation, the selected orders of ∂ h ∂ x and ∂ h ∂ y are 0, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 0.2, and 0.3. In Figures 1-15, h x and h y are ∂ h ∂ x and ∂ h ∂ y respectively. In Figures 1,  13-15, h x = h y = 0 stands for the MNLS perturbation wave stability curve, neglecting the bottom slope and to compare with other curves. Figures 1-12 and Figures 1-6 in the Supplementary Material indicate curves of |Im W| and the corresponding scatter diagrams of Im W = 0 with the invariable value of ∂ h ∂ x and the variable value of ∂ h ∂ y from 0 to 0.3. The discussions are as follows: 1. As Figures 1-12 shown, for a 0 = 0.196 and the invariable value of ∂ h ∂ x , the value of |Im W| increases as the value of ∂ h ∂ y rises. As the value of ∂ h ∂ y varying from 0 to 0.01, the value of |Im W| is maximum around the original point in wave number plane. As the value of ∂ h ∂ y is 0.1, 0.2 and 0.3, there is no maximum value of |Im W| around the original point in wave number plane. As the value of ∂ h ∂ x increasing from 0 to 0.001, the contours for |Im W| are similar. As ∂ h ∂ x = 0:01, the contours show obvious change by comparing with that As   10. In conclusion, as the value of ∂ h ∂ x increases, the curve of Im W = 0 approximates to m axis as ∂ h ∂ y < ∂ h ∂ x . The increment of depth variation in x direction contributes the the Stokes wave to be stable in or to parallel to y axis disturbed by small perturbation for ∂ h ∂ y < ∂ h ∂ x . The curve of Im W = 0 is the broken line composed by m = 1 and m axis for the value of ∂ h ∂ y = ∂ h ∂ x from 0.01 to 0.1. Figure 7 to Figure   n summary, as the value of ∂ h ∂ y increases, the curve of Im W = 0 approximates to l axis for ∂ h ∂ x < ∂ h ∂ y . The increment of depth variation in y direction helps the Stokes wave to be stable in or to   Figure 25 to Figure 60 in supplementary material indicate curves of |Im W| and corresponding scatter diagrams of Im W = 0 for the value of ∂ h ∂ x and the value of ∂ h ∂ y both varying from 0 to 0.3. It is indicated that ∂ h ∂ x begins to influence the shape of curve for Im W = 0 on the order of 0.000001 and to affect the value of |Im W| on the order of 0.01. It is shown ∂ h ∂ y begins to influence the shape of curve for Im W = 0 on the order of 0.00001 and affect the value of |Im W| on the order of 0.01. They show similar properties to the curves and corresponding scatter diagrams for a 0 = 0.196, with little value of |Im W| than that for a 0 = 0.196.
It is concluded the curve for Im W = 0 is more sensitive to depth variation terms than the curve of |Im W|. The curve for Im W = 0 is more sensitive to depth variation terms as a 0 = 0.0995 than that as a 0 = 0.196. The curve for Im W = 0 is more sensitive to ∂ h ∂ x than that to ∂ h ∂ y . Scatter maps of Im W = 0 for a 0 = 0.0995, a 0 = 0.196 and ∂ h ∂ x = 0 are indicated in Figure 13. The values of a 0 and ∂ h ∂ y determine the characteristics of curves of Im W = 0. The value of a 0 determines the intercept in l axis, with larger intercept for a 0 = 0.196 than for a 0 = 0.0995. The value of ∂ h ∂ y changes the amplitude of curve, with little amplitude for larger value of ∂ h ∂ y . Figures 14, 15 show the curves of Im W = 0 as a 0 = 0.196 and a 0 = 0.0995 for the value of ∂ h ∂ x = ∂ h ∂ y varying from 0 to 0.05. The scatters of Im W = 0 form smooth curves for ∂ h ∂ x = ∂ h ∂ y = 0. The scatters gather to be groups for ∂ h ∂ x = ∂ h ∂ y ≠ 0. A broken line is formed as ∂ h ∂ x = ∂ h ∂ y = 0:05, in the line of m = 1 and m axis.

CONCLUSIONS
On a finite slowly varying depth, the surface gravity wave equation is expanded to the fourth order by multiscale expansion in the narrowband range, and the TMNLS equation is obtained. When the time scale is longer than e -3 times of the dominant wave period or the space scale is larger than e -3 times of the dominant wavelength, the second depth derivative and square of the first depth derivative influence on the first-order wave height. Compared with the MNLS perturbation analysis results on the flat bottom when h x = h y = 0, the small perturbation's instability area by TMNLS is the entire wavenumber space, except for the curve satisfying Equation (28), which means that TMNLS increases the small perturbation's instability area by including depth variation terms to MNLS.
For a 0 = 0.196, ∂ h ∂ x starts to influence the shape of curve for Im W = 0 on the order of 0.00001 and to affect the curve of |Im W| on the order of 0.01. ∂ h ∂ y starts to influence the shape of curve for Im = 0 on the order of 0.0001 and affect the curve of |Im W| on the order of 0.01.
For a 0 = 0.0995, ∂ h ∂ x begins to influence the shape of curve for Im W = 0 on the order of 0.000001 and to affect the curve of |Im W| on the order of 0.01. ∂ h ∂ y begins to influence the shape of curve for Im W = 0 on the order of 0.00001 and affect the curve of |Im W| on the order of 0.01.
The curve for Im W = 0 is more sensitive to depth variation terms than the curve of |Im W|. The curve for Im W = 0 is more sensitive to depth variation terms as a 0 = 0.0995 than that as a 0 = 0.196. The curve for Im W = 0 is more sensitive to ∂ h ∂ x than that to ∂ h ∂ y . As the value of ∂ h ∂ x increases, the curve for Im W = 0 approximates to m axis as ∂ h ∂ y < ∂ h ∂ x . The increment of the value for depth variation in x direction contributes the Stokes wave to be stable in or paralleling m axis disturbed by small perturbation for ∂ h ∂ y < ∂ h ∂ x . The curve of Im W = 0 is the broken line composed by m = 1 and m axis for ∂ h ∂ y = ∂ h ∂ x ≥ 0:05. As the value of ∂ h ∂ y increases, the curve approximates to l axis for ∂ h ∂ x < ∂ h ∂ y . The increment of the value for depth variation in y direction contributes the Stokes wave to be stable in or paralleling l axis for ∂ h ∂ x < ∂ h ∂ y . The curve of Im W = 0 is a broken line combined by m = 1 and m axis for ∂ h ∂ x = ∂ h ∂ y ≥ 0:05.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.