**A Numerical Implementation of a Nonlinear Mild Slope Model for Shoaling Directional Waves**

#### **Justin R. Davis, Alex Sheremet, Miao Tian and Saurabh Saxena**

**Abstract:** We describe the numerical implementation of a phase-resolving, nonlinear spectral model for shoaling directional waves over a mild sloping beach with straight parallel isobaths. The model accounts for non-linear, quadratic (triad) wave interactions as well as shoaling and refraction. The model integrates the coupled, nonlinear hyperbolic evolution equations that describe the transformation of the complex Fourier amplitudes of the deep-water directional wave field. Because typical directional wave spectra (observed or produced by deep-water forecasting models such as WAVEWATCH III™) do not contain phase information, individual realizations are generated by associating a random phase to each Fourier mode. The approach provides a natural extension to the deep-water spectral wave models, and has the advantage of fully describing the shoaling wave stochastic process, *i.e.*, the evolution of both the variance and higher order statistics (phase correlations), the latter related to the evolution of the wave shape. The numerical implementation (a Fortran 95/2003 code) includes unidirectional (shore-perpendicular) propagation as a special case. Interoperability, both with post-processing programs (e.g., MATLAB/Tecplot 360) and future model coupling (e.g., offshore wave conditions from WAVEWATCH III™), is promoted by using NetCDF-4/HD5 formatted output files. The capabilities of the model are demonstrated using a JONSWAP spectrum with a cos2s directional distribution, for shore-perpendicular and oblique propagation. The simulated wave transformation under combined shoaling, refraction and nonlinear interactions shows the expected generation of directional harmonics of the spectral peak and of infragravity (frequency <0.05 Hz) waves. Current development efforts focus on analytic testing, development of additional physics modules essential for applications and validation with laboratory and field observations.

Reprinted from *J. Mar. Sci. Eng.* Cite as: Davis, J.R.; Sheremet, A.; Tian, M.; Saxena, S. A Numerical Implementation of a Nonlinear Mild Slope Model for Shoaling Directional Waves. *J. Mar. Sci. Eng.* **2014**, *2*, 140-158.

#### **1. Introduction**

As waves propagate into shallow water, they change from almost sinusoidal in deep water to a sawtooth like shape in the surf zone. Troughs become wide and shallow; crests peak and lean forward, eventually overturning and breaking. In the spectral domain, this evolution is expressed in energy transfers from the spectral peak to peak harmonics and low frequency (between 0.001 Hz and 0.02 Hz) waves, as well as the development of phase correlations across the spectrum. Wave-shape evolution and the generation of zero-frequency motions (mean flow, wave setup) have significant effects on nearshore sediment transport and inundation.

Modeling nonlinear shoaling is challenging. Off-the-shelf finite-depth spectral models (e.g., SWAN [1]) are typically based on variance balance equations originally developed for deep-water waves [2], and therefore cannot account for phase correlation effects. Describing directional nonlinear wave interactions is problematic in intermediate depth. Shallow water spectra are typically wide (containing harmonics and low frequency waves) precluding the use of simpler weak dispersion approximations (cubic Schrodinger equation e.g., [3–6] or Boussinesq approximations, e.g., [7–10]).

The fundamental challenge of modeling nonlinear shoaling in the spectral domain resides in the character of wave interactions. The basis of the spectral representation is the decomposition of the wave field into statistically independent (in the leading order) Fourier modes. For a flat bottom (water depth *h* = constant), this representation is formally

$$\begin{aligned} \eta(\mathbf{x}, t) &= \sum\_{n,m} a\_n e^{i(\mathbf{K}\_{n,m}\mathbf{x} - \omega\_n t)} + c.c. \\ \omega\_n^2 &= g K\_n \tanh(\mathbf{K}\_n h) \end{aligned} \tag{1}$$

where ߟ is the free surface displacement, *a* is the complex modal amplitude, ߱ is the radian frequency, and "*c.c.*" stands for "complex conjugate". The sum (used here to denote symbolically any type of superposition, either discrete or continuous) is carried out over frequencies (indexed by *n*). Different directions of propagation are represented here by the wave number vector *Kn,m* which depends on both the frequency index and an additional index *m*, specifying, say, the propagation angle. For a given ߱ and a given depth *h*, Equation (1) constrains the modulus of the wave number vector of modulus *Kn*. The efficiency of wave nonlinearities [3,11] depends on the system of equations describing the resonance state of *N* interacting modes.

$$\sum\_{j=1}^{N} \omega\_j = 0\tag{2}$$

$$\sum\_{s=1}^{N} \mathbf{K}\_s = \mathbf{0} \tag{3}$$

Note that Equation (3) is a system of two equations for the components of the horizontal wave number vector. With the additional constraint (Equation (2)), only two of the three scalar equations (Equations (2) and (3)) can be independent. A set of *N* modes that satisfy Equations (2) and (3) is said to interact resonantly; those that do not are called non-resonant. In the wave evolution equation, the efficiency of resonant *N*-wave interactions scales like *O(İ<sup>N</sup>í<sup>1</sup> )*, where *İ* is the characteristic wave slope. Non-resonant effects are weaker and dynamically less relevant (produce higher order bound waves). Due to the form of the dispersion relation (Equation (1)), the smallest number of modes that can be resonant is *N* = 4 (quadruplet, or four-wave interaction); triad interactions (*N* = 3) are non-resonant in any water depth [12]; however, they approach resonance in shallow water.

The statistics of wave evolution can be described in terms of competing effects of dispersion and nonlinearity [13,14]: nonlinearity builds phase correlations and skews the statistical distribution of the wave-field; dispersion breaks them and restores the symmetry of the distribution.

In deep water, the dispersive terms of the evolution equation are of order *İ*, while the competing leading-order nonlinearity (resonant four-wave interactions) is of order *İ<sup>3</sup> .* Consequently, the wave field is Gaussian in the leading order, with its statistics completely determined by second order moments (variance, power spectrum; [2,15]), hence the suitability of models based on energy-balance equations.

As water becomes shallow, dispersion weakens to order *İ<sup>2</sup>* while nonlinearity strengthens. Near-resonant triad interaction (order *İ<sup>2</sup>* ) becomes the leading order nonlinear mechanism (e.g., [16–20] and many others). The evolution is characterized by the broadening of the spectrum, and the generation of significant phase correlations across the spectrum (wave crests peak, wave fronts steepen). The waves are no longer Gaussian: wave statistics are no longer completely determined by second order moments (power spectrum) alone, and higher order moments and spectra (e.g., bispectra) become important. Evolution depends on both local sea state and wave history (history of phase correlations).

The dynamics of triad interaction in shallow water are poorly (or not at all) implemented in existing numerical models. For example, SWAN [1] arguably the most advanced coastal spectral model, is essentially built on a WAM [15] energy balance structure [2]. It implements a crude and unrealistic triad interaction parameterization [21], limited to approximating collinear second harmonic generation exclusively, with depth dependent interaction coefficients alone (*i.e.*, accounting only for local effects, and not for wave history). Important processes such as infragravity (IG) wave generation, recurrence effects, and spectrum widening are also ignored.

A deterministic, unidirectional but complete triad interaction formulation was first introduced by [16] based on the Boussinesq approximation. Agnon *et al.* [20] proposed a generalization for arbitrary depth based on the Nonlinear Mild Slope Equation (NMSE, [22,23]). Limited directionality can be introduced using the parabolic approximation (e.g., [24–26]). Hyperbolic forms for nearly planar beaches were developed by [20,27,28].

This paper describes the modeling techniques used to implement the hyperbolic form of the NMSE developed by [20,29] for directional three-wave interaction. Development of the model is presented first followed by a demonstration of the model's capabilities for shore-perpendicular and oblique wave propagation. Finally, a summary of the work is presented along with a discussion of future enhancements to the model.

## **2. Model Development**

#### *2.1. Nearshore Directional Waves*

In Equation (1), the directionality of mode *n* is expressed by the direction of the wave number *K*. The wave number vector is an invariant of propagation in deep water, and can be used to label directional modes. The wave number is considered an independent variable, with *Ȧ* given by Equation (1), and modes are identified by the wave number components or, in polar coordinates, by the pair *(K,ș)*, where *ș* is the angle of propagation. Thus, directional modes are represented by a two-dimensional parameter space (e.g., indices *n* and *m* in Equation (1)).

In the nearshore, *K* is no longer invariant, but the wave frequency typically is. If the beach has straight and parallel isobaths, the alongshore wave number provides a second invariant that can be used to complete the two degrees of freedom necessary for describing directional waves. Therefore, in the nearshore, the Fourier representation of Equation (1) can be replaced by

**176**

$$\eta(\mathbf{x},t) = \sum\_{n,m} a\_{n,m} e^{i\left(\int^x k\_{n,m} dx + \kappa\_m y - \omega\_n t\right)} + c.c.\tag{4}$$

$$K\_n^2 = k\_{n,m}^2 + \kappa\_m^2 \tag{5}$$

where *x* and *y* are the cross- and alongshore coordinates. The independent parameters are, in the approach, the frequency *fn* (or *Ȧn* = 2Ɋ*fn*), and the alongshore wave number *țm*. The wave number modulus *K* depends on the frequency through Equation (1), and the cross-shore wave number *k* is a function of both *f* and *ț* through Equation (5). A mode is therefore defined as the pair *(fn,țm)*—indexing modes rather than the independent parameters—and mode *J* is defined as the pair

$$J = (f\_l, \kappa\_l) \tag{6}$$

From Equation (6), for a given frequency *f*, and at a given cross-shore location *x*, progressive modes satisfy the condition

$$
\kappa^2 \le K^2(f, \mathfrak{x}) \tag{7}
$$

where *K(f,ț)* is the local wavenumber modulus, given by the linear dispersion relation, Equation (1). Modes that do not satisfy this relation in some nearshore domain are called trapped modes. The location *x0* at which *k* = *K(f,x0)* is called the "turning point". For simple (e.g., monotonic) beach profiles, shoreward of the turning point, trapped modes can acquire oscillatory behavior since *K*ĺ as *h*ĺ0.

#### *2.2. A Hyperbolic Nonlinear Mild Slope Equation (NMSE)*

The numerical model described here implements the formulation proposed by [27] (see also [20]) for the nonlinear evolution of directional waves over a mildly sloping beach. The stationary nonlinear mild-slope equation can be written as

$$\frac{\partial B\_{I}}{\partial \mathbf{x}} = i \{k\_{I} + i d\_{I} \} B\_{I} + \sum\_{P,Q} \{-\delta\_{I,Q+P} + 2\delta\_{I,Q-P} \} i W\_{I,\pm P,Q} B\_{\pm P} B\_{Q} \tag{8}$$

$$\mathbf{a}\_{I} = \left[ \mathbf{c}^{1/2} a e^{i \int k dx} \right]\_{I} \mathbf{c}\_{I} = \left( \frac{k}{K} \mathbf{C} \right)\_{I} \tag{9}$$

where *J*, *P*, and *Q* are directional Fourier modes in the sense of Equation (6), and *cj* is the cross-shore component of the model group velocity *C*. The parameter *dj* represents dissipation and/or growth processes, such as breaking, wind input, bottom friction, and others. In Equation (8), ߜ is the Kronecker symbol, for example,

$$\mathcal{S}\_{f,Q\pm P} = \begin{cases} 1 & f = Q \pm P \\ 0 & f \neq Q \pm P \end{cases} \text{ with } \pm P = (\pm f\_p, \kappa\_p) \tag{10}$$

where the equality *J* = *P* has the regular meaning for ordered pairs, *i.e.*, *fJ* = *fP* and *ț<sup>J</sup>* = *țP*. Only modes that satisfy the selection criteria given by Equation (10) are allowed to contribute to the nonlinear terms. Triads satisfying *J* = *P* + *Q* ("sum" interaction) are responsible for transferring energy toward high frequencies; difference interactions *J* = í*P* + *Q* transfer energy toward low frequencies. An example of a sum-interaction triad is shown in Figure 1. With the notation ߪ ൌ ߱Ȁඥ݃, the interaction coefficient is

$$\mathcal{W}\_{I,\pm P,Q} = \frac{1}{8} \frac{\sqrt{g} \sigma\_{\parallel}}{\sigma\_P \sigma\_Q} \left( c\_I c\_P c\_Q \right)^{-1/2} \left( \pm 2 \mathcal{K}\_P \cdot \mathcal{K}\_Q + \mathcal{K}\_P^2 \frac{\sigma\_Q}{\sigma\_{\parallel}} \pm \mathcal{K}\_Q^2 \frac{\sigma\_P}{\sigma\_{\parallel}} + \sigma\_P^2 \sigma\_Q^2 \mp \sigma\_{\parallel}^2 \sigma\_P \sigma\_Q \right) \tag{11}$$

**Figure 1.** An example of a sum-interaction triad *J*(*j,s*) = *Q*(*q,v*) + *P*(*p,u*), *j* = *q* + *p* and *s* = *v* + *u* with *j* = 3, *q* = 1, *p* = 2 and *s* = 4, *v* = 3, *u* = 1.

Equation (8) represents the Nonlinear Mild Slope Equation (NMSE) model. The NMSE is hyperbolic and describes wave shoaling, refraction, and three-wave nonlinear interactions. The unknown function *BJ* is related to the energy flux in the cross-shore direction, หܤห <sup>ଶ</sup> ൌ หܽห ଶ ܿ. The linear part of the equation describes the conservation of the cross-shore component of the modal energy flux (the alongshore component is conserved trivially). The quadratic term represents the contribution of three-wave interaction to mode evolution and redistributes energy flux between modes.

The numerical implementation of Equation (8) is restricted only to triads that are close enough to resonance, as measured by the "detuning" parameter

$$\mu = \frac{\left|k\_I \pm k\_P - k\_Q\right|}{\left|k\_I\right|} = \frac{\left|\Delta\_{I,Q \pm P}k\right|}{\left|k\_I\right|}\tag{12}$$

The parameter *ȝ* compares the wavelength of the nonlinear term with the wavelength of mode *J* (on the left-hand side of the equation). If *ȝ* >> 1, the oscillations of the nonlinear term are fast and result in a small (second-order) "bound wave" correction to mode *J* that can be calculated approximately as

$$B\_I^{bound} = \sum\_{\substack{P,Q\\\Delta k = \mathcal{O}(k\_I)}} \left(\delta\_{I,Q+P} - 2\delta\_{I,Q-P}\right) \frac{W\_{I,\pm P,Q}}{\Delta\_{I,\pm P,Q}k} B\_{\pm P} B\_Q \tag{13}$$

This approximation becomes singular as *ȝ*ĺ0. This occurs as *h*ĺ0, *i.e.*, triads approach resonance as the water becomes shallow. In this case, the oscillation of the nonlinear term is slow and the equation has to be integrated numerically. In principle, the numerical solver should be able to handle triads with arbitrary values of *ȝ*. In practice, however, numerical calculations for *ȝ* = *O(1)* are slow because the model has to resolve fast oscillations that yield small contributions to the derivative. Controlling the errors becomes increasingly difficult for larger values of *ȝ* and the benefit of the effort becomes negligible. Because of that, an efficient numerical implementation of Equation (8) would limit the integration to triads characterized by *ȝ* < *ȝc*, for some critical value of *ȝc*, with bound waves computed using Equation (13). The numerical simulations shown here use *ȝ<sup>c</sup>* = 0.5, while the bound waves are ignored (will be included in future modifications of the code).

Equation (8) is valid strictly for progressive waves. Trapped modes are not allowed to interact in the spatial domain where their cross-shore structure is exponential, but are allowed in the domain where they have oscillatory behavior. The NMSE model is phase resolving, in that it requires initial values for both modal amplitudes and modal phases.

Equation (8) reduces to the unidirectional equation for a mild sloping beach [16,20,25] if all the modes propagate perpendicular to the shoreline, *i.e.*, for all *ț<sup>J</sup>* = 0. Numerical simulations using the unidirectional hyperbolic NMSE [27,29] have been extensively verified against both single-triad analytic solutions as well as laboratory and field observations.

In the current implementation of the model, the only dissipation mechanism used is depth-limited wave breaking, based on the frequency dependent parameterization developed by [29,30], with dissipation uniformly distributed over all directions.

#### *2.3. Model Discretization and Computational Grid*

In Fourier series representation, the frequency-alongshore wave number is discretized as

$$\begin{aligned} f\_j &= j \Delta f, \quad & 0 \le j \le N\\ \kappa\_s &= s \Delta \kappa, \quad & -M \le s \le M \end{aligned} \tag{14}$$

and a directional mode index *J* is a pair of indices *J* = *(j,s)*. For a triad of interacting modes *J*, *P*, and *Q*, the selection criterion given in Equation (10) can be written as indicial equations

$$\begin{array}{rcl}j & = & \pm p + q \quad \text{( $f$ )}\\s & = & u + v \quad \text{( $\kappa$ )}\end{array} \tag{15}$$

where

$$\begin{array}{rcl} \pm j &=& (\pm j, \mathfrak{s}) \\ \pm P &=& (\pm p, u) \\ Q &=& (q, v) \end{array} \tag{16}$$

for a given *f*, the effective *ț*-range of the allowable modes is limited by Equation (7), and can vary with depth. As the maximum extents of *f* and ț are known, a list of all possible triads can be created before shoreward marching of the solver begins. The matrix of triads involving a given mode *J* (in the left-hand side of Equation (8)) and all the allowable modes *P* and *Q* (right-hand side of Equation (8)) is

$$J = \begin{cases} (\mathcal{P}, Q)\_{11}^{J} & \text{(\(\mathcal{P}, Q\)}\_{12}^{J} & \cdots & \text{(\(\mathcal{P}, Q\)}\_{1n}^{J} & \text{(\(f = Q + P\)} \\ (\mathcal{P}, Q)\_{21}^{J} & \text{(\(\mathcal{P}, Q\)}\_{22}^{J} & \cdots & \text{(\(\mathcal{P}, Q\)}\_{2m}^{J} & \text{(\(f = Q - P\)} \end{cases} \tag{17}$$

Because the selection criteria (Equation (10)) are invariants of propagation, the interacting triads can be pre-computed for a given *(f,ț)* matrix.

#### *2.4. Solution Algorithm*

The NMSE represents a coupled system of ܰ ൈ ሺʹܯ ͳሻ complex ODEs, a hyperbolic initial value problem. These equations are solved using the Vode ODE Solver [31] using a non-stiff Adams method. Although the NMSE is written in complex form, for purposes of solving, the equation is split into real and imaginary components (doubling the number of equations to solve simultaneously) thus enabling the double precision (8 byte) real version of Vode to be used. An overview of the solution algorithm is shown in Figure 2.

**Figure 2.** Solution algorithm.

As waves propagate into shallower water, trapped modes (modes for which *k* = *K* at some depth) become active and participate in the interaction. A trapped mode is considered inactive (*i.e.*, not allowed to interact with other modes) in the domain where *k* > *K*, but becomes active if *k* < *K* (*i.e.*, shoreward of the turning point for monotonic beach profiles). Triads containing inactive modes are disabled; therefore, the maximum effective alongshore wave number (ߢ௫ ൌ ܯοߢ (depends on the local (cross-shore position) depth and frequency. The conditions that determine *M* are

$$M\_{local} = M\{f\_j, h\} = M\_j(h) = \left\lfloor \frac{|K\_j|}{\Delta \kappa} \right\rfloor \tag{18}$$

where *K* is determined by the dispersion relation (Equation (1)) and ہ ۂ refers to the integer value.

The maximum effective wave number increases with frequency and decreasing depth. The variation of *M* with depth and frequency can be handled using two different strategies: (1) Using the minimum depth and highest frequency, the maximum *M* for the entire domain can be determined; (2) Increase *M* as the solution marches toward the shore. The current implementation of the model uses Approach 1. This approach will result in sparse matrices (wasting some computer memory) but the triad interaction patterns can be defined once for all runs and there is no need to dynamically modify *M* as the solution marches toward the shore. For each evaluation of the derivatives, it is only necessary to determine whether all modes of the triad are active. Approach 2 is expected to result in dense matrices but also to significantly complicate coding as array sizes would vary as a function of cross-shore position. As with Approach 1, it is still necessary to determine whether all modes of the

triad are active. A pseudocode representation of how the model calculates the right hand side of the NMSE (Equation (8)) is shown in the Appendix.

#### *2.5. Model Input/Output*

An overview of how data is imported into and exported out of the model is shown in Figure 3. The code solves the NMSE as a Monte-Carlo simulation. Typically, available offshore wave information consists of directional spectral density of free-surface variance. Offshore modal amplitudes are provided in a simple text file which either contains the complex amplitudes (includes a phase for each mode) for each "realization" to be simulated, or a spectrum can be provided and the model will use a Random Phase Approximation (RPA) to generate phases for a user-defined number of realizations. Model output is provided in a NetCDF-4/HDF5 output format using NetCDF [32]. Metadata provided in the output file is compliant to CF-1.6 [33]. The variable defining the number of realizations being simulated has the NetCDF length "UNLIMITED". Thus, for a given set of simulation parameters and offshore wave conditions, realizations can be performed independently and their output files easily combined, e.g., with "ncrcat" [34]. It is also noted that while the input file is currently a simple text file, the model could easily be setup to read a NetCDF file using the same metadata convention as the output file.

The NMSE describes the shoaling transformation of a stationary directional wave field from deep into shallow water. The details of the discretization of the frequency and alongshore wave number spaces are user defined. In the current implementation, the model resolves the shortest wave with 10 points (spectral cutoff frequency is 1/5 the Nyquist frequency) and the alongshore wave field at a resolution of 5 m. The solver (Vode) used to integrate Equation (8) uses an adaptive algorithm that implicitly discretizes internally the cross-shore domain according to the accuracy requested for the solution. The user only controls the locations for the solution output. In the simulations presented, for the purpose of describing the details of each realizations (see Section 3), values of the solution are generated every 5 m in the cross-shore.

#### **Figure 3.** Model input and output.


**Table 1.** Synthetic scenario parameters.
