A numerical model for edge waves on a compound slope

An edge wave is a kind of surface gravity wave basically travelling along a shoaling beach. Based on the periodic assumption in the longshore direction, a second order ordinary differential equation is obtained for numerical simulation of the cross-shore surface elevation. Given parameters at the shoreline, a cross-shore elevation profile is obtained through integration with fourth-order Runge–Kutta technique. For a compound slope, a longshore wavenumber is obtained by following a geometrical approach and solving a transcendental equation with an asymptotic method. Numerical results on uniform and compound sloping beaches with different wave periods, slope angles, modes and turning point positions are presented. Some special scenarios, which cannot be predicted by analytical models are also discussed.


Introduction
Edge wave can be treated as a unique kind of surface gravity waves, trapped in the nearshore region, propagating along the shoreline (Mei et al., 2005). It has long been considered significant in the formation of beach cusps under reflective wave conditions (Guza and Inman, 1975). Partly due to difficulties in field observations and laboratory reproduction, most studies on edge-wave are of theoretical interest; whereas the practical use is also suggested by recent studies on applying data-driven techniques to predicting coastal morphological changes (Alvarez and Pan, 2016;Horrillo-Caraballo et al., 2016). In classical wave theory, analytical solution for edge waves exists for only two types of beach topography: one is for the uniform sloping beach, where the exact analytical solution of zero-mode and fullmode are given by Stokes (1846) and Ursell (1952) respectively in the context of full water and that of the infinitemode is given by Eckart (1951) in the context of shallow water. The other is for the exponential beach, where the approximate solution is given by Ball (1967) in the context of shallow water. Though explicit, Ball's dispersion relationship is thought to be too complicated.
As for more complex topography, Munk et al. (1964) used a numerical scheme to solve for wave motions on the California shelf. His approach is to approximate the shelf lim y→∞ (y) = 0 profile with a number of horizontal steps and to match elevation and mass flux between steps. Darbyshire (1977) used another numerical model to find edge wave lengths on a real beach. Holman and Bowen (1979) used a scheme similar to Darbyshire's which solved the equations for any form of bottom profile from an approximating analytical function to a specified array of depth data with linear interpolation between data points. Given the velocity potential at the shoreline φ(0) and its partial derivative φ'(0), φ on each grid point can be obtained through integration using the Runge-Kutta technique. Values of k x , the longshore wavenumber, are picked out when the condition is satisfied in numerical experiments. However, this method of defining k x by trials is somewhat rough, and the governing equation is shallow water equation (SWE), which fails in simulating high mode edge waves accurately. In this paper, the model is set up in the context of full water, while Holman and Bowen's method of integration is employed. Governing equations and numerical model are introduced in Sections 2 and 3, respectively. Shen et al. (1968) have used geometrical optics to study edge waves for the first time. However, it is not until Schäffer and Jonsson (1992) have presented this approach in a simple way as opposed to the complex mathematical treatment by Shen et al. (1968) that it is known to more people.
In Schäffer and Jonsson's paper, the equation is constructed with a simple physical meaning, that is, the edge wave should be trapped in the nearshore region. With the application of Snell's law, explicit dispersion relationships of edge waves on uniform sloping and exponentially sloping beach are derived respectively in shallow water theory. Dispersion relationship in full water theory is also derived for uniform sloping beach. Most of them match well with the corresponding analytical solutions. Since this approach has a clear physics-based background, here it is followed to construct a transcendental equation of k x , where the definite integral for uniform sloping beach is converted to integral of piecewise functions for a compound slope. While k x is known, it drives the numerical model established in Section 3. The detailed derivation of k x is given in Section 4. Results and analyses of the numerical experiment are given in Section 5. Section 6 is for conclusions.

Governing equations
It is known that the dynamics of the surface wave propagating over a seabed of mild slope is well described by the mild-slope equation. For monochromatic waves, Smith and Sprinks (1975) derived the equation below for numerical simulation: where p=cc g , , and ω, c, c g , η are respectively the angular frequency, the phase velocity, the group velocity, and the surface elevation. h and φ are respectively the water depth and velocity potential (hereinafter). c and c g can be determined by the linear dispersion relationship below: . (2) For edge waves, η may be assumed to be in a sinusoidal longshore form with a wave number k x in the x direction (x denotes the long-shore direction; y denotes the cross-shore direction; z is the vertical coordinate, hereinafter) then we can rewrite Eq. (1) as:

A numerical model
Employing the shallow water equation, Holman and Bowen (1979) developed a numerical model to simulate edge waves on complex beach profiles. In their paper, boundary conditions at the shoreline can be determined by expanding the depth h and velocity potential φ in series form where σ denotes the angular frequency in that paper. By substituting Eqs. (5) and (6) into Eq. (7) and equating coefficients of powers of y, the boundary condition can be determined in the form: Now we try to extend this approach to the full water theory for numerical simulation. Similarly, on slowly varying topography, the variables in Eq. (10) may be expanded in powers of y, where the subscript s denotes variables at the shoreline, and α 1 , α 2 , β 1 , β 2 , γ 1 , and γ 2 are parameters. Substituting Eqs.
In fact, variable values such as (k i ) s , c s , (c g ) s , (cc g ) s at the shoreline, where the theoretical depth is zero, are determined by introducing a virtual depth h s =ε at the shoreline. Since c g ≈c and , So is generally determined by (h y ) s rather than h s . Besides, Eq. (4) can be equally well expressed as two coupled first-order differential equations: Now the problem lies in the determination of the longshore wavenumber k x . With definite values of k x , Eq. (18) can be solved numerically by starting from the boundary conditions in Eq. (16) to find for each discrete node, using the Runge-Kutta technique.

Determination of k x on a compound slope
For waves of very little steepness on shallow beaches, dissipation can be of minor importance, and full reflection may be assumed (Schäffer and Jonsson, 1992). In classical wave theory, if the longshore component of the wave number k x (which is independent of the depth from Snell's Law) is larger than the deep water wave number (or equivalently ω 2 ≤gk x ), the wave cannot exist in deep water. A caustic will develop at some depth h=h c for y=y c , and the wave will be trapped to the coast. For edge waves, the whole interference pattern is confined to the near-shore, and then at the caustic a crest from the incoming wave train must be cancelled by a trough in the reflected one. (c.f. Figs. 5 and 6 in the paper of Schäffer and Jonsson, 1992; Note that here we use y rather than x to denote the cross-shore coordinate, which is adopted in the original paper). In this case, the equation below is established: Let identify non-dimensional quantities and let the length scale be (k ∞ is the wavenumber at infinity, hereinafter), and then the dimensionless wavenumber and water depth can be defined as , . So the dispersion relationship can be expressed as: Following Schäffer and Jonsson's approach, where α is the local angle of incidence. As is taken to be positive, Eq. (24) may be integrated to give: For uniform sloping beaches, Eq. (27) is proved to be equivalent with a much simpler form after some mathematical treatments, that is (28) Now this approach is applied to a case of a compound slope. The topography is depicted in Fig. 1.
Assuming that the caustic of the refracted wave is at h=h * >h 0 (If h * ≤h 0 , then the integration in Eq. (27) will stay the form and the calculated will be the same as that of a uniform sloping case) and following the derivation, we can obtain the following equation: where . Here, the definite integral in Eq.
(27) is replaced with a piecewise intergral in Eq. (29). According to the linear dispersion relationship, k varies from a smaller value in deep water to a larger value in shallow water. At the caustic h=h * , 30) or in the non-dimensional form, (31) In this way, the solution of Eq. (29) is confined to a small In this paper, we solve Eq. (29) with a simple but effective numerical approach. We consider as a constant parameter, then can be taken out of the unknown integral in Eq. (29) and will be cancelled by the left term of the equation. Set function F(b): where 0≤b≤1. F(b) is monotonic, so F(b)=0 can be solved numerically by zero-crossing method. With b and τ being known, the definite intergral can be calculated with an adaptive Gauss-Kronrod quadrature method. Results of this approximation method show its effectiveness.

Model verification and result analysis
For uniform sloping scenarios, three types of results are presented: "Analytical", "Numerical-Analytical" and "Numerical-Geometrical", or "A", "NA" and "NG" for short, which denote Ursell's analytical solution, integration solution with Ursell's dispersion relationship ( ), and integration solution with geometrical dispersion relationship ( ), respectively. The gridspacing k ∞ ∆y is set to be 0.005 and the calculation is terminated when one of the following two conditions is satisfied: Here denotes non-dimensional surface elevation beyond the computational domain, obtained from Ursell's theory.
The results are arranged in three groups, i.e. (1) same wave period, same mode, and different slope angles; (2) same wave period, same slope angle, and different modes; (3) same slope angle, same mode, and different wave periods, which are depicted in Figs. 2, 3 and 4, respectively. In general, Fig. 2 shows that NA results are in quite good coincidence with A results, especially in smaller angle scenarios. For larger angle cases (such as ), NA results deviate from A results at last, regardless of integration precision. This anomaly may partly result from Eq. (1), which is accurate only with bed slopes ranging from 0 to 1/3 (Booij, 1983), i.e. a slowlyvarying topography condition. Fig. 2 also shows that NG results have better performance in smaller angle cases. It is supposed to be related to the accuracy of dispersion relationship. For , 0.005, so NG results under this condition are proved to be highly reliable. Fig. 3 shows that NG results have slightly better performance with higher modes, given a certain slope angle. It should be noted that shallow water theory ( ) also yields accurate results on milder slope beaches. However, as fails with larger n, NG results perform much better, such as n=10 in Fig. 3, while shallow water or Boussinesq models almost fail at n≥8 (e.g. Fig. 2 in Sheremet and Guza, 1999). Fig. 4 shows that there is little difference among NG results with different wave periods. The reason is that throughout our derivation, there is no long-wave assumption, and the longshore wavenumber has been non-dimensionalized by the deep water wave number at the beginning. For simplicity, we choose the period to be 8 s for all the rest tests. Besides, Fig. 2. Non-dimensional elevation profile with same wave period (T=8 s), same mode (n=2), and different slope angles (β=π/10, π/22, π/32).  it should also be noted that NG approach leads to finite modes on certain topography, while shallow water theory leads to infinite modes, which is untrue.
For compound sloping scenarios, the grid-spacing is also 0.005, while the calculation is terminated when both and are satisfied.
The results are similarly arranged in three groups, i.e. (1) same wave period, same mode, same turning point position, and different slope angle ratios; (2) same wave period, same slope angle pair, same turning point position, and different modes; (3) same wave period, same slope angle pair, same mode, and different turning point positions, which are depicted in Figs. 5, 6 and 7, respectively. An obvious character is that all the edge-waves on compound sloping beaches decay to zero in the deep water zone. However, it has been derived that non-dimensional surface elevation of the highest mode (Ursell wave) in deep water is (-1) n for uniform sloping beach of critical angles (i.e. , n=1, 2,..., Zhang et al., 2009). One possible explanation is the property of the caustic. Further calculation shows that the caustic of all kinds of edge waves lies beyond the last extreme point of the elevation profile. In other words, edge waves directly decay to zero beyond the caustic, where the cross-shore wave number is imaginary (Wang et al., 2013). Ursell wave on a critical slope reaches its caustic at infinity that is why it will not decay. For compound sloping cases, a milder slope in deep water makes a finite caustic possible, so it decays at last. Fig. 5 shows that the profile is basically determined by β 0 and shifts a little due to different angle pairs. For both scenarios with β 0 /β 1 <1 and β 0 /β 1 >1, the profile extends to a marginally larger area with the decreasing ratio β 0 /β 1 . It may be explained directly by the trapping mechanism, where milder slope changes the direction of reflection wave more efficiently than steeper slope does, thus leading to smaller active area. N = π 4β 0 − 1 2 Fig. 6 shows the profiles of different modes with other parameters being fixed. An interesting finding is that higher mode (such as n=11 illustrated here), not allowed in fullwater analytical model, may also exist. Therefore, it may not be enough to determine the highest mode number for compound cases with . Fig. 7 shows the effect of different turning point positions. It seems that different positions would change the profile by "stretching" or "compressing", but would not change the overall shape.

Conclusions
(1) Based on the mild-slope equation and the periodic assumption for waves in the longshore direction, a second order ordinary differential equation is obtained for numerical simulation of coastal edge waves. A geometrical approach is applied in the determination of wavenumber on a compound slope. By using integration method, we are able to present a numerical model of edge waves on a slowlyvarying compound sloping beach. The new model yields good results consistent with analytical solutions on most uniform sloping beaches, especially milder slopes. Numerical precision is not sensitive to modes or wave periods as expected in shallow-water models. Higher modes may lead to better precision sometimes.
(2) For all compound slope cases, the results match well   with the evanescent assumption in the offshore direction. The surface elevation profile is primarily determined by the near-shore slope angle, secondly by the turning point position, and lastly by the angle ratio. The last two parameters shift the profile in a manner of "stretching" or "compressing".
(3) For compound slope cases with critical near-shore slope angle, it is shown that mode higher than Ursell wave may exist in cases with a milder connecting slope.