# Numerical Method for Calculating Fourier Coefficients and Properties of Water Waves with Shear Current and Vorticity in Finite Depth

## Article information

## Abstract

Many numerical methods have been developed since 1961, but unresolved issues remain. This study developed a numerical method to address these issues and determine the coefficients and properties of rotational waves with a shear current in a finite water depth. The number of unknown constants was reduced significantly by introducing a wavelength-independent coordinate system. The reference depth was calculated independently using the shooting method. Therefore, there was no need for partial derivatives with respect to the wavelength and the reference depth, which simplified the numerical formulation. This method had less than half of the unknown constants of the other method because Newton’s method only determines the coefficients. The breaking limit was calculated for verification, and the result agreed with the Miche formula. The water particle velocities were calculated, and the results were consistent with the experimental data. Dispersion relations were calculated, and the results are consistent with other numerical findings. The convergence of this method was examined. Although the required series order was reduced significantly, the total error was smaller, with a faster convergence speed.

**Keywords:**Fourier approximation; Fenton’s method; Newton’s method; Shooting method; Incompressible Euler equations

## 1. Introduction

Since Chappelear (1961) developed an early numerical method, many other numerical methods have been derived (Dean, 1965; Rienecker and Fenton, 1981; Chaplin, 1979; Cokelet, 1977; Fenton, 1988). Rienecker and Fenton (1981) solved the coefficients and some wave properties directly using Newton’s method. This method was further simplified by Fenton (1988). The major modification was that all the required partial differentials were calculated numerically (Tao et al., 2007). Although neither of the above numerical methods could be applied to waves in deep water (Fenton, 1990), Shin (2023) calculated deep water breaking waves with an error of less than 2.04×10^{−3} percent, demonstrating that Fourier approximation is effective even for deep water waves. The objective of this study was to extend Shin’s method (Shin, 2023) to include rotational water waves with a shear current in a finite depth (Shin, 2022). The basic concept of Shin’s method (Shin, 2023) for establishing a set of equations was similar to Fendon’s method (Rienecker and Fenton, 1981; Fenton, 1988), and some issues with Fenton’s method have been addressed. The first issue is inherent in Newton’s method, which necessitates an initial solution in close proximity to the final solution. In common with other versions of the Fourier approximation method (Chappelear, 1961; Dean, 1965; Chaplin, 1979), it is sometimes necessary to solve a sequence of lower waves, extrapolating forward in height steps until the desired height is reached (Fenton, 1988). These methods can occasionally converge to the wrong solution with very long waves (Fenton, 1988). This paper shows that regression analysis can easily avoid this issue (Shin, 2019).

The second issue arises from the coordinate system. Because Fenton’s method adopts the moving coordinate system proposed by Dean (1965), the abscissa is the position from the crest in the range, −*L*≤*x*≤*L*. Therefore, the relative position with respect to the wavelength, L, is unknown because the wavelength is unknown. As a result, the number of unknown variables to be determined in Fenton’s method is more than double the number in this study. Furthermore, a nonlinear equation for calculating the wavelength was coupled to the set of equations for calculating the coefficients. This coupling complicates the problem because Newton’s method requires partial derivatives with respect to the wavelength. The issue can be avoided easily using the dimensionless coordinate system in Fig. 1 because it is independent of the wavelength. The abscissa is the phase in the range −*π*≤*β*≤*π*, and is known. After determining the steepness, *S*, the wave-number is calculated as *k*=*S*/*H* because wave height, *H* is known. The steepness is calculated using the wave height condition, which is integrated automatically into the set of equations because it is linear. Therefore, it is not necessary to calculate the partial derivatives with regard to the wavelength in this study.

The third issue arises from the water depth condition used to calculate the reference depth, *D*, which is unknown. Fenton’s method requires partial derivatives with respect to the reference depth because the water depth condition is coupled to the set of nonlinear equations used to calculate the coefficients. Most errors in Fenton’s method occur when the water depth condition is calculated numerically. The required series order must be increased to reduce the error, and Simpson’s rule was used for the numerical integration. This issue was avoided by calculating the reference depth independently using the shooting method in the study. In numerical analysis, the shooting method is used to solve a boundary value problem by transforming it into an initial value problem. A set of “*n*” equations with “*n*” unknowns can be converted to a set of “*n*−1” equations with “*n−*1” unknowns because one variable can be calculated independently, assuming it is a known value in advance. The assumed value was determined using the shooting method. Although it has the disadvantage of requiring the calculation of one procedure as two procedures, the numerical formulation simplifies the process and significantly reduces the total number of steps needed to converge to a complete solution. The reference depth is not coupled to the set of equations for calculating the Fourier coefficients because the water depth condition is calculated independently.

The fourth issue is a method for calculating the partial derivatives required in Newton’s method. While all the necessary partial derivatives were obtained numerically in Fenton (1988), they were calculated analytically by tensor analysis in the present study. Therefore, there are no errors in the partial derivatives in this study.

Although the above-mentioned studies are irrational flows, this study is rotational flow based on Shin (2022). The numerical method used to calculate the Fourier coefficients and the related wave properties was not presented by Shin (2022) but is presented in the present study. Shin (2018) presented the numerical method for irrotational waves. The method was extended to rotational waves in this study. Shin (2022) reported the complete solutions to the Euler equations and calculated the waves by Le Mehaute et al. (1968). The fluid velocities calculated by irrotational wave theories (Fenton, 1990; Tao et al., 2007) and by rotational wave theory (Shin, 2022) were compared with the experimental data reported by Le Mehaute et al. (1968) on the same graphs. Unlike the irrotational waves, the results of the rotational wave theory agreed well with the experimental data. The calculated wave-breaking limit was in good accordance with the Miche formula. While Rienecker and Fenton (1981), Fenton (1988), Tao et al. (2007), Cokelet (1977), and Vadden-Broeck and Schwarts (1979) calculated the dispersion relations for some waves in a relative depth of 0.7, this study calculated these relations for many waves with different heights at various depths. The results were then extrapolated to all waves using regression analysis (Shin, 2019). The dispersion relations for a relative depth of 0.7 were compared and showed good accordance with the other results. The convergence of this study was tested according to the series orders. Even for relatively small orders, the error was significantly lower than Fenton’s method and Tao et al. (2007).

## 2. The Complete Solution to the Euler Equations

This chapter summarizes the solution reported by Shin (2022). The two coordinate systems were adopted. The first was a conventional coordinate system (*t*, *x*, *y*), as reported by Dean et al. (1984), in which the origin is located at a fixed point on the still water line (SWL). The *x*-axis is in the direction of wave propagation; the *y*-axis points vertically upward, and *t* is the time. The fluid domain is bounded by a flat bed, *y* = −*h* (water depth), and by a free surface *y* = *η* (*t*, *x*). The second is a dimensionless coordinate system, (*β α*), as shown in Fig. 1, which is a stationary frame. The origin is located under the crest of the bed. Therefore, the wave profile is fixed with a periodic, even function in the system. The abscissa is phase *β* = *kx* − *ωt* in −*π* ≤ *β* ≤ *π*, where *k* = 2*π*/*L* is the wave number and *ω* = 2*π*/*T* is the angular frequency. Here, *T* is the wave period; *L* is the wavelength, and the ordinate is the dimensionless elevation from the bed, *α* = *k*(*y* + *h*) in 0 ≤ *α* ≤ *γ*, where *γ* = *k*(*η* + *h*) is the dimensionless free surface elevation from the bed. The term *ζ* = *γ* − *D* is a dimensionless free surface elevation measured from the reference line (the horizontal line passing through two points on the free surface at two phases *β* = ±*π*/2. The reference depth *D* = *k*(*η _{o}* +

*h*) is the dimensionless distance from the bed to the reference line where

*η*=

_{o}*η*(±

*π*/2). Because the coordinate system is independent of the wavelength, it has several advantages, as presented by Shin (2023). The dimensionless quantity of a flow field

*f*is denoted by

*f̄*except for those defined separately, such as

*S*and

*θ*. Using the notation, the stream function is defined by

*ψ*=

*ψ̄ω/k*

^{2}.The horizontal and vertical velocities are defined by

*u*=

*Cū*and

*v*=

*Cv̄*, respectively, where

*C*=

*ω*/

*k*is the celerity. The vorticity is defined by

*Ω*=

*ωΩ̄*. The pressure is defined by

*p*=

*ρC*

^{2}

*p̄*, where

*ρ*is the water density. The dimensionless wave height (steepness) is defined by

*S*=

*kH*. The linear steepness

*θ*=

*ω*

^{2}

*H*/

*g*(where

*g*is the gravity) is a constant for a particular wave. When

*h*→ ∞ and

*H*→ 0,

*S*→

*θ*, which gives the dispersion relation (

*ω*

^{2}=

*gk*) of deep water linear waves. The stream function to satisfy the incompressible Euler equations is as follows:

*N*is the series order; and

*c*are the Fourier coefficients. The first and second terms describe a current. On the other hand, the two terms can be deleted for an irrotational wave (ε = 0). The vorticity is presented as follows:

_{n}This solution also contains irrotational waves, where the constant, *ε* is referred to as the strength of vorticity because *Ω̄* = 0 when *ε* = 0. When *ε* is non-positive, the vorticity has the same direction as the motion of the particles. The iso-stream and iso-vorticity lines can be calculated using Newton’s method in the Appendix. The free surface is an iso-vorticity line given by *Ω̄*(*β*,*γ*) = *C*_{1}*ε*, where *C*_{1} is a constant. Substituting *Ω̄*(*β*,*γ*) = *C*_{1}*ε* into Eq. (2), the wave profile is determined to satisfy the KFSBC (Kinematic free surface boundary condition) as follows:

The constant *C*_{1} is determined from the definition of the reference line in Fig. 1; that is, *γ*(±*π*/2) = *D*. The profile is an implicit function that can be calculated using Newton’s method in the Appendix. When *ε* → 0, Eq. (3) yields the profile defined by Shin (2018). Moreover, the sea bed is an iso-vorticity line given by *Ω̄*(*β*,*0*) = 0. Differentiating Eq. (1) with respect to *α* and *β*, the horizontal velocity becomes the following.

The vertical velocity is as follows:

Substituting Eqs. (1)–(5) to the Euler momentum equations, Bernoulli’s principle for rotational flow can be calculated as follows.

*U*

*̄*

^{2}=

*V*

*̄*

^{2}−

*ε*(

*α*−

*ψ̄*)

^{2}, and

*Q*is a constant. A mechanical analog to irrotational and rotational flows was depicted by considering a carnival Ferris wheel in Fig. 2.11 in Dean et al. (1984). While the irrotational motion of chairs involves rotating around the center of the Ferris wheel, the rotational motion of the chairs involves both revolving around the center of the Ferris wheel and rotating around their center. Therefore, the total kinetic energy of the rotational motion is represented by the sum of the revolving kinetic energy around the center of the Ferris wheel and the rotating kinetic energy around their center. Analogous to the Ferris wheel, the total kinetic energy of rotational flow is represented by

*U*

*̄*

^{2}=

*V*

*̄*

^{2}−

*ε*(

*α*−

*ψ̄*)

^{2}, where the first term corresponds to the revolving kinetic energy around the center of the Ferris wheel, and the second term corresponds to the rotating kinetic energy around its center. When

*ε*= 0,

*U*

*̄*

^{2}=

*V*

*̄*

^{2}= (1 −

*u*

*̄*)

^{2}+

*v*

*̄*

^{2}, which is the kinetic energy for irrotational flow. In contrast to Baddour et al. (1998),

*ε*is not positive because kinetic energy functions should be positive-definite functions. Therefore, let

*σ*is a non-negative real number. Hence, 〈0〉 =

*σi*. Substituting Eqs. (1), (4), and (5) into the relation,

*U*

*̄*

^{2}=

*V*

*̄*

^{2}−

*ε*(

*α*−

*ψ̄*)

^{2}

*,*the dynamic pressure (kinetic energy)

*U*

*̄*

^{2}is as follows:

*Z*

_{1}= (〈

*n*〉〈

*m*〉 +

*nm*− ε)/4

*W; Z*

_{2}= (〈

*n*〉〈

*m*〉 −

*nm*+ ε)/4

*W*;

*Z*

_{3}= (〈

*n*〉〈

*m*〉 −

*nm*− ε)/4

*W*;

*W*= cosh〈

*n*〉

*D*cosh〈

*m*〉

*D*. If Λ

_{0}

*is the component of column vector {Λ*

_{n}_{0}

*} and*

_{n}*c*is the component of the column vector {

_{n}*c*}, the second term on the right side in Eq. (7) becomes an inner product of the two vectors. Accordingly,

_{n}*Λ*can be considered a component of the square matrix [

_{nm}*Λ*] in

_{nm}*N*-dimensional space. Therefore, the third term on the right side of Eq. (7) is presented as {

*c*}

_{n}*[*

^{t}*Λ*]{

_{nm}*c*}, where {

_{m}*c*}

_{n}*is a row vector. This concept helps formulate numerical calculations in this study. Considering the average speed*

^{t}*Ū*over a wave period on the sea bed, the constant

_{b}*c*

_{0}is determined as follows:

For the average speed *Ū _{s}* over a wave period on still water, the strength of vorticity is as follows:

From Stokes’ breaking criterion, |*Ū _{b}*| ≤ 1 and |

*Ū*| ≤ 1. The Bernoulli’s constant

_{s}*Q*is determined from the definition of the reference line (i.e.,

*γ*(π/2) =

*D*) and the DFSBC (Dynamic free surface boundary condition), i.e.,

*p̄*(

*π*/2,

*D*) = 0 because the point (

*π*/2,

*D*) is on the free surface. Therefore, the pressure field is

By applying the DFSBC, i.e., *p̄*(*β*, *γ*) = 0 on *α* = *γ* to Eq. (12), the other wave profile can be determined as follows:

The wave height condition is defined as

The water depth condition is defined as

## 3. Newton’s Method for Coefficients and Steepness

The solution contains *N*+2 unknown constants: *N* Fourier coefficients, *c _{n}*, steepness,

*S*, and reference depth,

*D*. Eqs. (3) and (13) present two wave profiles with implicit functions. The two wave profiles are even functions considering the coordinate system. The Fourier series of a periodic even function is presented as

*π*≤

*x*≤

*π*. The coefficients are generally determined from the orthogonality of trigonometric functions. The method is inapplicable when

*f*(

*x*) is an implicit function, as expressed in Eqs. (3) and (13). The problem can be solved by converting the series to a set of algebraic equations, which are obtained by calculating the series at some phases instead of all phases. This is done by replacing the infinite series with a truncated series, i.e.,

*m*= 1, …,

*N*. There are

*N*algebraic equations for calculating the coefficients,

*c*, because cos

_{n}*nx*and

_{m}*f*(

*x*) are known. Hence, there are a set of

_{m}*N*equations, and Eqs. (3) and (13) are equal to each other at

*N*+1 phases (Note that Eqs. (3) and (13) are already equal to each other at phase ±

*π*/2 because the integral constants were determined to satisfy the definition of the reference line) and two equations, i.e., Eqs. (14) and (15). Therefore, there are

*N*+2 equations to determine the unknown constants.

Referring to Eqs (3)–(5), the denominators are expressed as transcendental functions of the reference depth. When determining the reference depth using Newton’s method in conjunction with the coefficients, it is necessary to calculate the partial derivative with respect to the reference depth. Furthermore, Eq. (15) only permits a numerical approach. This difficulty can be overcome by calculating the reference depth independently using the shooting method, assuming that the reference depth is known. Based on the assumption, Eq. (15) is independent of the other equations and will be calculated in the next chapter. In addition to this assumption, *N* parameters, *X _{m}*, are introduced in this chapter to simplify the numerical formulation. Unlike Fenton’s method,

*X*are merely parameters for calculating the coefficients and the steepness in this chapter. After determining the coefficients,

_{m}*X*are no longer used. Hence, the wave elevation is denoted with

_{m}*X*instead of

_{m}*ζ*in this chapter.

_{m}When the wave profile *ζ* and the reference depth *D* are prescribed, it is possible to convert Eq. (3) into a set of linear equations for the coefficients. Because the wave profile is an even function, the phases *β _{m}*(

*m*= 1,2, …,

*N*) are considered in the range 0 ≤

*β*≤

_{m}*π*and

*β*≠

_{m}*π*/2 because the two wave profiles are already equal at phase

*β*=

*π*/2

*.*In addition,

*β*

_{1}= 0 and

*β*=

_{N}*π*. Denoting

*ζ*(

*β*) as

_{m}*X*, i.e.,

_{m}*X*=

_{m}*ζ*(

*β*), Eq. (3) is presented.

_{m}The summation convention is considered in Eq. (16). The repeated subscript “*n*” is a dummy subscript, whereas *K _{mn}* is a second-order tensor, and the terms

*c*and

_{n}*X*are vectors in N-dimensional space. From Eq. (3), the component of the tensor

_{m}*K*is presented as

_{mn}*X̄*=

_{m}*X*+

_{m}*D.*Because 〈0〉=

*σi*, the component of the first-order tensor,

*Y*is presented as

_{m}The summation convention is not applied when the component of a tensor is presented as in Eqs. (17) and (18). Using the inverse tensor *G _{mn}* of the tensor

*K*, the solution to Eq. (16) is determined easily as follows.

_{mn}Eq. (13) is also evaluated at the same phase. A set of *N* nonlinear equations for calculating the *N* parameters *X _{n}* is obtained by substituting Eq. (19) into Eq. (13) at the same phase because the wave height condition in Eq. (14) is expressed as

*S*=

*X*

_{1}−

*X*in this approach. For Newton’s method, the error vector

_{N}*E*is defined from Eq. (13) as follows:

_{m}From Eqs. (8) and (9), the second-order tensor *A _{pn}* and the third-order tensor

*B*are defined as

_{pnm}Substituting Eq. (19) into Eq. (20)*E _{p}* = 0 gives a set of nonlinear equations for calculating the parameters,

*X*. Therefore, unlike Fenton’s method (Rienecker and Fenton, 1981; Fenton, 1988), the set of N nonlinear equations is solved using Newton’s method in this study. Therefore, when using the same series order, the number of unknown constants is less than half of that in Fenton’s method.

_{q}Denoting partial derivatives of a tensor by making use of commas and indices as
*X _{q}* points in terms of

*ΔX*and ignoring higher-order terms, the following set of linear equations is obtained:

_{q}Because

The superscript (*r*) means the step of Newton’s method. It is clear that all the steps are *r ^{th}* in all equations except Eq. (24). Thus, for the simplification of equations, the superscript (

*r*) was omitted in all the equations except Eq. (24). Differentiating Eq. (20) with respect to

*X*

_{q}*,*the partial derivative

*E*of the error vector is

_{p,q}Note that *c _{n,q}*

*c*=

_{m}*c*

_{n}*c*Because

_{m,q}*S*=

*X*

_{1}−

*X*,

_{N}Differentiating Eq. (16) with respect to *X _{p}* results in

Multiplying Eq. (27) by the tensor *G _{im}*, the partial derivative of the coefficient is determined as follows.

*G*

_{im}*K*=

_{mn}*δ*and

_{in}*δ*

_{in}*c*=

_{n,p}*c*and

_{i,p}*δ*is the second-order isotropic tensor. Differentiating Eq. (17) with respect to

_{in}*X*, the components of the third-order tensor

_{p}*K*are presented as

_{mn,p}Differentiating Eq. (18) with respect to *X _{p}*, the component of the second-order tensor

*Y*is presented as

_{m,p}Differentiating Eq. (21) with respect to *X _{q}*, the component of the third-order tensor

*A*is expressed as follows:

_{pn,q}Differentiating Eq. (22) with respect to *X _{q}*, the component of the forth-order tensor

*B*is

_{pnm,q}Based on the assumption that the reference depth is known, Eq. (16) is linear and decoupled from Eq. (20), unlike in Fenton’s method (Rienecker and Fenton, 1981; Fenton, 1988). Dividing the inner product of the error vector, *E _{m}* by N, the mean squared error (MSE) in the DFSBC is defined as

The MSE is used as a criterion for determining the convergence of the Newton method.

Shin (2019) was used as the first step,
*f*(*β*) with a common ratio, *λe ^{iβ}* (where “i” is the imaginary unit), closely resembles the shape of a water wave. A reliable approximation of the profile of a water wave can be achieved by shifting the function down by “

*f*(

*π*/2)” to align with the reference line in Fig. 1, normalizing the result by its height and multiplying the result by the actual wave amplitude “

*S*/2” to meet the wave height condition, as reported by Shin (2019).

The profile satisfies the definition of the reference depth, i.e., *ζ*(±*π*/2) = 0, the wave height condition in Eq. (14), and the water depth condition in Eq. (15). Replacing *ζ* with
*β _{m}* is presented as follows:

*S*

^{(0)}and the shape factor

*λ*are presented in Shin (2019), and (0) means the first step in Newton’s method. The function is much simpler than the other nonlinear waves because it contains only two parameters, i.e., wave steepness,

*S*and the shape factor,

*λ*. When the shape factor approaches 0, it becomes an Airy wave. When the shape factor approaches 0.5, it becomes a fifth-order Stokes wave. When the shape factor approaches 1, it becomes a Solitary wave. The shape factor is a positive real number less than 1. Because the function was derived independently from the Euler equations, the two parameters were determined to minimize the error in two boundary conditions on the free surface using a variational method. Because Eq. (34) is very close to the complete solution, Newton’s method rapidly and converges for all waves.

## 4. Shooting Method for Reference Depth

In contrast to the assumption considered in Ch. 3, the reference depth is closely linked to the coefficients and is unknown. Therefore, while maintaining the validity of the method presented in Ch. 3, the shooting method is suitable for determining the reference depth. The reference depth is adjusted until the calculated water depth from the water depth condition matches the actual value using the shooting method. The secant method is used as a root-finding algorithm. In each step of the secant method, substituting the calculated coefficients and the prescribed reference depth into Eqs. (A1)–(A4), the wave elevations

*q*) stands for the q’

^{th}step of the secant method, and

*ζ*(

_{i}= ζ*β*).

_{i}*ζ*is an even function. Hence,

*β*

_{1}= 0;

*β*=(

_{i}*i*− 1)

*π*/

*M*; and

*β*

_{M}_{+1}=

*π*. In Fenton (1988), Eq. (35) must be calculated using only

*N*data of

*X*because Eq. (35) is coupled to the other equations. Therefore,

_{i}*M*=

*N*in Fenton (1988). On the other hand, M is independent of

*N*in this study. Eq. (35) can be calculated accurately compared to Fenton (1988) because M can be freely increased in this study. The water depth is calculated as follows:

The reference depth is calculated using the secant method as follows.

*D*

^{(0)}is determined using the regression result in Shin (2019), and

*h*is the actual water depth. The actual value is in the vicinity of

*D*

^{(0)}. Hence,

*D*

^{(1)}is selected with a value close to

*D*

^{(0)}, e.g.,

*D*

^{(1)}= 1.001

*D*

^{(0)}. For each step, the error in the water depth is as follows:

## 5. Calculation Procedure and Results

### 5.1 Calculation Procedure

Fig. 2 presents a flowchart. While there is only one do-loop in Fenton’s method, the chart comprises three do-loops. The outer do-loop implements the secant method for calculating the reference depth, as described in Ch. 4. The first inner do-loop implements the Newton method for calculating the coefficients and steepness, as described in Ch. 3. The other inner do-loop implements Newton’s method for calculating the wave profile, as described in the Appendix.

Initial approximations of the shape factor, steepness, and reference depth are first calculated by substituting the period, height, and depth into the regression analysis results presented by Shin (2019). The initial approximations of the wave profile are calculated using Eq. (34).

The first inner do-loop is performed. The coefficients,
*q*” represents the step of the secant method in Ch. 4. Substituting
*D*^{(}^{q}^{)} into Eqs. (20) and (25),

The second inner do-loop is then performed to calculate the free surface elevations at more phases than those considered in the first inner do-loop. The elevations,
*r*” represents the step of Newton’s method.

By substituting
*h*^{(}^{q}^{)} is calculated using Eq. (36). The reference depth *D*^{(}^{q}^{+1)} is calculated using Eq. (38). Replacing *D*^{(}^{q}^{)} with *D*^{(}^{q}^{+1)}, the outer do-loop should be repeated until the calculated depth converges to the actual depth. The first step solution of Newton’s method in the *q*+1 step of the secant method,
*D*^{(0)} is accurate, the secant method converges within three steps to the acceptable value whose error is less than 10^{−4} percent. The outputs are the coefficients, the reference depth, and the steepness.

### 5.2 Application Range of this Study

The coefficients and wave properties are functions of two variables: the reference depth and the linear steepness. The functions are presented with closed-form solutions in Shin (2019), which are the regression results determined using numerically calculated data for 2091 waves, as shown below. The domain of the functions is {(*θ, D*) | 0 ≤ *θ* < 1 and 0 < *D*}. A subset of the domain, {(*θ, D*) | 0 ≤ *θ* < 1 and 0 < *D* < 5}, was considered instead of the entire domain because waves are similar to deep water waves when *D* > 5. Seventeen cases of the reference depth were considered by dividing the domain by *D* = 0.1, 0.15, 0.2, ..., 1.0, 1.2, 1.5, 2, 3, 4, and 5. For each reference depth, the linear steepness was considered from 0 to 0.98 at 0.01 intervals, from 0.98 to 0.989 at 0.001 intervals, from 0.989 to 0.9899 at 0.0001 intervals, and from 0.9899 to 0.98997 at 0.00001 intervals, resulting in 123 cases of linear steepness. Therefore, 17 × 123 = 2091 irrotational waves were calculated with *N* = 35 and *M* = 180. The limits can be determined easily because the Newton method presented in the Appendix diverges at the crest when the steepness reaches its limit. These limits were represented by 17 red solid circle symbols in Fig. 3, which are calculated with *D* = 0.1, 0.15, 0.2, ..., 1.0, 1.2, 1.5, 2, 3, 4, and 5 in sequence from left to right. The solid curve is the regression analysis result determined using the 17 data. As shown in Fig. 3, the breaking limit is in good accordance with the Miche formula. Because the 2091 waves are dimensionless waves, they represent all waves below the breaking limit. Therefore, this study confirms stable convergence characteristics for all water waves.

### 5.3 Water Particle Velocities Under the Crest

The eight waves considered in Le Méhauté et al. (1968) were calculated. Fig. 3 and Table 1 present the waves with solid diamond symbols. The results were compared with experimental data in Fig. 2 by Shin (2022). Fig. 4 shows the wave (h) during the eight waves. Referring to Fig. 2 of Shin (2022), the other waves give the same pattern, as shown in Fig. 4. The irrotational results were calculated with *ε* = 0. The regression analysis results were calculated by Shin (2019).

Rienecker and Fenton (1981) also calculated the waves, as shown in Fig. 2 of that work. The waves were also calculated by Fenton (1990), and the results are shown in Fig. 3(c) of that work. The cnoidal theory is the new cnoidal theory by Fenton (1990), which closely agrees with the Fourier approximation by Fenton (1988). Tao et al. (2001) also calculated the waves; the results are shown in Fig. 2 of that work. HAM (Homotopy analysis method) by Tao et al. (2001) results closely agree with the Fourier approximation reported by Rienceker et al. (1981). Therefore, the Fourier approximation and HAM results are not presented in Fig. 4. This study agrees with the experimental data, unlike those obtained using HAM or Fenton’s method.

### 5.4 Dispersion Relation

Table 3 of Tao et al. (2007) provides a detailed comparison of the HAM solutions of waves in finite water depth and the results of Cokelet (1977) and Rienecker and Fenton (1981). The results are given in terms of the non-dimensional phase speed squared *kC*^{2}/*g* at various values of *H*/*h*, for a constant value of exp(−*kR*/*C*) = 0.5 (*R* is the unit span denoting the total volume rate of flow underneath the steady wave per unit length in a direction normal to the *x, z* plane), corresponding to a wavelength to water depth ratio of *L*/*h* ≈ 9 (Tao et al., 2007). Table 1 of Rienecker and Fenton (1981) presents the same comparison between the solutions of Rienecker and Fenton (1981) for *N* = 16, 32, and 64, the results of Cokelet (1977), and the results of Vanden-Broeck & Schwarts.

The non-dimensional phase *kC*^{2}/*g* is presented with physical quantities defined in this study as follows:

*H*/*h* is also presented with physical quantities defined in this study.

*kh*= 2

*π*/9 in Eq. (41) because

*L*/

*h*≈ 9. Therefore, the comparison shows the relationship between steepness, S, and linear steepness,

*θ*for the relative water depth,

*kh*= 2

*π*/9. Hence, the comparison shows the dispersion relation of waves in a finite water depth,

*kh*= 2

*π*/9. Table 2 lists the relation. The data in the first and second columns were cited from Rienecker and Fenton (1981), which are the results calculated using

*N*= 64. The data in the fourth column was obtained by multiplying the data in the second column by 2

*π*/9, as stated in Eq. (41). The data in the third column was obtained by multiplying the data in the first column by the data in the fourth column, as stated in Eq. (40). Shin (2018, 2019) also calculated the relation for various water depths,

*D*= 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞. For a given reference depth, the linear steepness

*θ*was considered with a highly dense interval from 0 to the application limit of this study. Fig. 5 presents the calculation results. The 15 curves are the results calculated with

*D*= 0.1, 0.15, 0.2, ... , 1.0, 1.2, 1.5, 2, and ∞ in sequence from left to right. The peaks of each curve represent the limit of wave breaking. The nine data in Table 2 are also presented with solid circles in Fig. 5, which agree with the curve for

*D*= 0.7. The small difference results from the difference between

*kh*= 2

*π*/9 and

*D*= 0.7. The relative depth for

*D*= 0.7, is presented as

*kh*=

*D*−

*kη*

_{0}= 0.7 −

*kη*

_{0}> 2

*π*/9. because

*D*=

*kh*+

*kη*

_{0}and kη

_{0}is negative. Therefore, the

*D*= 0.7 curve lies slightly to the right of the eight data because

*kη*

_{0}is very small.

### 5.5 Error Check of this Study

The numerical method focuses on minimizing the errors in the DFSBC defined in Eq. (20) and the water depth conditions defined in Eq. (15). The convergence of this study according to the series order was analyzed by calculating the eight waves considered in Le Méhauté et al. (1968) using *N* = 8, 10, 16, and 35. Fig. 6 presents the calculated wave profiles for the waves. The abscissa represents the phase *β* = *kx* − *ωt* in −180° ≤ *β* ≤ 180°, while the ordinate represents the dimensionless elevation, *kη*. The errors are summarized in Tables 3 and 4. The first step is the result of Shin (2019). Table 3 lists the error of Shin (2019). There are three steps for calculating the water depth condition using the secant method and three steps for the Newton method at each secant step. Therefore, there are 10 steps; Table 4 presents the errors for step 10. The blanks indicate that they are unsuitable. The series order, *N* = 8, is unsuitable for waves (c) and (d), and the series order, *N* = 10, is unsuitable for wave (d) because they do not satisfy the condition for regular waves. Hence, there is no local extrema between the crest and the trough in a regular wave. The deeper the water depth and the lower the wave height, the smaller the number of series orders required. If the series order is increased beyond the optimal number, the convergence of this study decreases, as presented by Shin (2023). Even for very long waves, the convergence is very stable, unlike other methods.

## 6. Conclusions

This paper reported a numerical method to calculate the Fourier coefficient, steepness, and reference depth parameter. The method is much simpler than Fenton’s method or HAM. The partial derivatives with respect to wavelength and reference depth, as well as certain parameters, were eliminated utilizing a dimensionless coordinate system. The water depth condition was calculated independently using the shooting method in this study.

## Notes

The authors declare that they have no conflict of interest.

## References

## Appendices

#### Newton’s Method for the Wave Profile

Note that the coefficient *c _{n}* and the reference depth

*D*are known in this Appendix. An implicit function

*f*(

*β*,

*ζ*) = 0 can be considered as an equation with respect to the dependent variable

*ζ*because the independent variable β is known. Thus, using Newton’s method, the explicit representation of the function is

*r*) stands for the step of Newton’s method in this Appendix and

*f*(

*β*,

*ζ*) is defined as follows using Eq. (3):

Because 0 ≤ |ζ| < 1 for all waves, the first step solution *ζ*^{(1)} is calculated with the power series expansion of *f* in *ζ* and ignoring the higher order terms than the second order. As *f*(*β*, *ζ*)= 0, the power series expansion is a quadratic equation with respect to *ζ*. Then we have

The other root of the quadratic equation is trivial. Eq. (A4) has the limit that is determined with the discriminant of the quadratic equation. On the other hand, the limit is greater than the wave-breaking limit. Therefore, Eq. (A4) is valid for all waves. Iso-streamlines and iso-vorticity lines are calculated by applying the method presented in this Appendix to Eq. (1) and Eq. (2), respectively. An iso-streamline for the specified value of *ψ̄* = constant can be obtained by applying Eq. (A5) to the Newton method instead of Eq. (A3) and replacing *ζ* with *α*.

The iso-vorticity line for the specified value of *Ω̄ =* constant can also be obtained by applying Eq. (A6) to the Newton method instead of Eq. (A3) and replacing *ζ* with *α*.