**3. Nonlinear Shoaling of Two JONSWAP Spectra: Shore-Perpendicular and Oblique Propagation**

We demonstrate the capabilities of the model with two shoaling tests over a plane beach of 0.03 slope. The offshore spectra are standard directional JONSWAP spectra propagating shore-perpendicular in the first test, and obliquely in the second. The parameters used here (Table 1) are typical for long Eastern Pacific swells; however, the directional spread is probably exaggerated in the simulations. The JONSWAP spectral shape (maybe not entirely realistic for representing the incoming waves at the deep end of the simulation domain) is used here solely to illustrate the capabilities of the current model and the use of the RPA for simulating the shoaling transformation of a deep-water variance density distribution. Because the JONSWAP spectrum does not contain any variance in the infragravity band, the second-order bound spectrum associated with the deep water swell was computed using Equation (13). A summary of the simulation parameters for both scenarios is shown in Table 1.

For input-output purposes, the numerical model requires mapping the directional wave information between the model *(f,ț)* grid and the standard frequency-angle representation *(f,ș)* (e.g., as used in WAVEWATCH III™ [35]). The existence of turning points makes the implementation of the mapping procedure sensitive to the bathymetric profile. For a given frequency, turning points (Section 2.1) are cross-shore locations where additional modes are introduced into the system (become active). Their effect on the geometry of the computational grid is illustrated in Figure 4. The active computational grid is limited to the band defined by *|ț|* < *K(f,h)*, where *K(f,h)* is given by Equation (1). In deep water, this is a narrow band (widening toward higher frequencies, Figure 4a). As the water depth decreases, the band widens (additional modes become active), and becomes triangular in shape (Figure 4c) and extends into higher alongshore wave numbers as the shallow water boundary approaches *|ț|* < *Kshallow* = *Ȧ(gh)*<sup>í</sup>1/2. As the limiting alongshore wave number increases, the frequency-angle representation degrades slightly (in Figure 4d, the mapped grid does not cover the entire available *(f,ș)* domain).

**Figure 4.** Directional characteristics of the frequency-alongshore wave number *(f,ț)* representation in comparison with the standard (deep-water) frequency-angle *(f,ș)* representation, in 15 m water depth (upper panels) and in 3 m water depth (lower panels). (**a**,**c**) The *(f,ț)* representation; contours show the corresponding angles of propagation with respect to shore-parallel. The shaded area in (**a**) marks trapped-wave modes (*i.e.*, modes that have the turning point between 15 m and 3 m water depth). The nodes of the *(f,ț)* grid are marked by blue points.

Designing the computational grid for applications poses thus the additional challenge of balancing the conflicting needs for resolving wide propagation angles (large limiting *ț*) at high resolution (small *ț* increments), and for keeping the number of triads described reasonably small for numerical integrations. The need for wide angles is non-trivial: for example, it is straightforward to check that directional difference triads containing two nearly collinear, shore-normal swell modes can excite a low-frequency wave that propagates nearly parallel to the shoreline. Note also that a significant fraction of the computational grid is never used.

Mapping the directional spectral distribution between *(f,ț)* and *(f,ș)* space*s,* shown in Figure 5, consists of two steps in each direction: a direct mapping of the modal amplitudes from the uniform grid (or angle, Figure 5a or alongshore wave number, Figure 5c) onto a non-uniform one in the complementary space (Figure 5b or 5d), and a re-sampling (interpolation) of the non-uniformly spaced values into the uniform grid. All transformations are designed to preserve the frequency spectrum (*i.e.*, the directional spectrum integrated over either angles or wave number).

**Figure 5.** Illustration of the mapping of the directional JONSWAP spectrum from the *(f,ș)* space onto the *(f,ț)* space, and back for the shore-perpendicular spectrum. Upper panels: directional spectra in different representations; lower panels: corresponding frequency spectrum. The transformation preserves the frequency spectrum. (**a**) Standard *(f,ș)* representation of the JONSWAP directional spectrum at 15 m isobaths; (**b**) Direct map from *(f,ș)* to *(f,ț)*. The resulting grid in *ț* is not uniform; (**c**) Spectrum re-sampled in the uniform *ț* grid used for computation; (**d**) Spectrum directly mapped back to *(f,ș)* space. The angle grid is not uniform. Units of the variance density contour plots are arbitrary.

The evolution of a total of ܰ ൈ ሺʹܯ ͳሻ ൌ ͳͲͲ ൈ ͳ ൌ ͳͲͲ possible (however, some high-frequency trapped modes never become active) directional modes are simulated (Figure 4). For each scenario, simulations are performed with both the full, and the linearized version of Equation (8). Note that the present implementation of the model only includes the linear and triad nonlinear evolution "engine" and wave breaking, with no additional physics (e.g., wave setup), that would be essential for realistic modeling of wave propagation in the nearshore.

The simulations shown here test the representation capability of the *(f,ț)* grid as well as illustrate the directional effects of nonlinear shoaling. The initial spectra at the deep-water end of the domain (15 m water depth) are shown in Figures 6a and 7a. Linear runs (Figures 6b and 7b) show the expected refraction effect of decreasing directional spread, with modes slowly turning around toward shore-perpendicular propagation. The main nonlinear effects (clearly visible in Figures 6c and 7c) are energy transfers from the peak to (a) peak harmonics, and (b) low-frequency infragravity modes. For oblique propagation, artifacts of the resolution of the *k*-grid are visible in the deep-water spectrum (Figure 7a), but become less severe as the waves refract and the grid coverage in the frequency-angle space increases with decreasing water depth. Note that infragravity waves (frequency < 0.05 Hz) are significantly more directionally spread (approximately 60 degrees) than the rest of the spectrum (approximately 30 degrees for swell and 15 degrees for the shortest waves represented).

**Figure 6.** Evolution of the directional shore-perpendicular JONSWAP spectrum (see Table 1). (**a**) Initial spectrum in 15 m water depth; (**b**) Linear evolution (3 m water depth); (**c**) Nonlinear evolution (3 m water depth). Simulations are averages over *N* = 100 random phase realizations.

**)**

**Figure 6.** *Cont.* 

**Figure 7.** Evolution of the directional oblique JONSWAP spectrum (see Table 1). (**a**) Initial spectrum in 15 m water depth; (**b**) Linear evolution (3 m water depth); (**c**) Nonlinear evolution (3 m water depth). Simulations are averages over *N* = 100 random phase realizations.

**)**

Figures 8 and 9 show the free surface elevation corresponding to one of the realizations used to estimate the spectra in Figures 6 and 7. The figures illustrate the change in the wave shape caused by the excitation of the phase-correlated harmonics of the spectral peak. A comparison of the linear (Figure 8a) and nonlinear (Figure 8b) oblique wave field clearly shows the steepening of the wave front. Both the shore-perpendicular and the oblique propagation realizations generate a significant infragravity field, with heights between 0.2 and 0.4 m. This effect is mainly a nonlinear shoaling effect (the linear shoaling of the initial bound infragravity band accounts for about 5 cm of the heights).

**Figure 8.** Contours of the simulated free surface elevation field for shore-perpendicular propagation, corresponding to the spectrum shown in Figure 6 at a fixed (arbitrary) time. The parameters of the offshore spectrum are given in Table 1. (**a**) Linear model; (**b**) Nonlinear model; (**c**) Infragravity waves (*f* < 0.05 Hz) generated during nonlinear shoaling.

**Figure 9.** Contours of the simulated free surface elevation field for oblique propagation, corresponding to the spectrum shown in Figure 7 at a fixed (arbitrary) time. The parameters of the offshore spectrum are given in Table 1. (**a**) Linear model; (**b**) Nonlinear model; (**c**) Infragravity waves (*f* < 0.05 Hz) generated during nonlinear shoaling.

#### **4. Summary and Conclusions**

Based on the equations derived in [27], a Fortran 95/2003 code was written to solve the NMSE. The model (as of Version 1-38) contains approximately 32,000 lines of code, 60% of which is the Vode ODE solver and 8% are testing routines. The model has been compiled using several different Fortran compilers (GNU 4.8 and Intel 2013) and executed under several different LINUX platforms (Ubuntu and RHEL). This implementation of the NMSE is a phase-resolving, spectral *(f,k)* model that describes wave evolution over a beach with straight and parallel isobaths. The model solves set of hyperbolic equations for a two-dimensional surface gravity wave field approaching a beach accounting for non-linear, quadratic triad interactions, shoaling and refraction. The capabilities of the model are illustrated for shore-perpendicular and oblique propagation of a directional JONSWAP spectrum. The results show wide-angle infragravity (IG) generation as well as significant energy

transfers toward high-frequencies modes. The runs demonstrate the essential role directionality plays in nonlinear shoaling, especially in the generation of directionally spread infragravity waves.

With the "engine" of a wave model developed, future research goals include: (1) Reformulating the governing equation from its current mathematically relevant form to one which is more numerically adept. By solving an equation better suited to floating point arithmetic (e.g., a form which minimizes rounding error), we should be able to improve the model's stability, accuracy and energy flux conservation; (2) Adding modules for relevant physics such as setup (feedback to the total depth) and wind effects; (3) Performing additional numerical tests. Comparing the model to known analytic solutions (e.g., individual triad interactions) as well as observational directional wave data will provide for better model verification and validation. Additional IG scenarios will let us better understand how they are generated. Finally, future applications of the model will include making it publicly available to the community as well as incorporating it into phase-averaged deep-water models such as WAVEWATCH III™ using the RPA mechanism described here.

## **Acknowledgments**

This research is supported by U. S. Office of Naval Research Grants N00014-10-1-0805, N00014-10-1-0389, N00014-12-1-0220, and N00014-13-1-0620, the U. S. Strategic Environmental Research and Development Program (SERDP) Grant N00173-13-1-G012, and the University of Florida. The authors also wish to acknowledge the University of Florida Research Computing for providing computational resources and support that have contributed to the research results reported in this publication. URL: http://researchcomputing.ufl.edu.

#### **Conflicts of Interest**

The authors declare no conflict of interest.

#### **References**


#### **Appendix: Pseudocode for Nonlinear Derivatives**

```
######################################################################### 