## Abstract

Large-amplitude mode-2 nonlinear internal waves were observed in 250-m-deep water on the Australian North West shelf. Wave amplitudes were derived from temperature measurements using three through-the-water-column moorings spaced 600 m apart in a triangular configuration. The moorings were deployed for 2 months during the transition period between the tropical monsoon and the dry season. The site had a 25–30-m-amplitude mode-1 internal tide that essentially followed the spring–neap tidal cycle. Regular mode-2 nonlinear wave trains with amplitudes exceeding 25 m, with the largest event exceeding 50 m, were also observed at the site. Overturning was observed during several mode-2 events, and the relatively high wave Froude number and steepness (0.15) suggested kinematic (convective) instability was likely to be the driving mechanism. The presence of the mode-2 waves was not correlated with the tidal forcing but rather occurred when the nonlinear steepening length scale was smaller than the distance from the generation region to the observation site. This steepening length scale is inversely proportional to the nonlinear parameter in the Korteweg–de Vries equation, and it varied by at least one order of magnitude under the evolving background thermal stratification over the observation period. Despite the complexity of the internal waves in the region, the nonlinear steepening length was shown to be a reliable indicator for the formation of large-amplitude mode-2 waves and the rarer occurrence of mode-1 large-amplitude waves. A local mode-2 generation mechanism caused by a beam interacting with a pycnocline is demonstrated using a fully nonlinear numerical solution.

## 1. Introduction

Mode-1 nonlinear internal waves (NLIWs) are often seen in the ocean, and a growing number of studies have also observed mode-2 waves in shallower ocean regions. These include observations in the northern South China Sea (Yang et al. 2009), the New Jersey shelf (Shroyer et al. 2010), and the Mascarene Ridge (da Silva et al. 2015). Mode-2 NLIWs have also been observed in Lakes Biwa and Kinneret (Boegman et al. 2003). Mode-2 NLIWs are characterized by either locally diverging isotherms forming a bulge or locally converging isotherms forming a contraction. The isotherms thus have either a locally convex or concave appearance to an observer in the upper water column (Yang et al. 2010) and hence the terms convex or concave waveforms. Mode-2 NLIWs are of particular interest as they can potentially create ocean conditions susceptible to higher rates of local dissipation and vertical mixing in the pycnocline. Furthermore, they can more easily form recirculating cores that efficiently transport tracers over large horizontal distances in the pycnocline (e.g., Deepwell and Stastna 2016).

Mode-2 waves lead to enhanced mixing through a combination of both shear and kinematic instability mechanisms [see, e.g., Helfrich and Melville (2006) and Thorpe (2018) for reviews on the topic of internal wave stability]. By virtue of their smaller horizontal and vertical wavenumber, mode-2 NLIWs create higher vertical shear *S* and wave steepness than mode-1 waves of equal amplitude. This, in turn, creates lower gradient Richardson numbers (Ri), where and *N* is the buoyancy frequency, increasing the likelihood of local shear instability (e.g., Thorpe 2018). Overturning via kinematic (also termed as convective) instability (Lamb 2002) has also been shown to occur within mode-2 NLIWs in numerical experiments (Deepwell and Stastna 2016). This mechanism has been observed in mode-1 NLIWs in the ocean when the wave induced velocity exceeds the phase speed resulting in isopycnal overturning and a trapped recirculating core (e.g., Zhang and Alford 2015). Shroyer et al. (2010) found higher dissipation rates for mode-2 NLIWs in the coastal ocean. They inferred dissipation time scales of h for mode-2 waves versus h for mode-1 solitary waves based on ship-tracked observations on the New Jersey shelf.

As well as breaking, mode-2 waves can also interact with other waves and lose energy. Stastna et al. (2015) studied the interaction between mode-1 and mode-2 NLIWs using a fully nonlinear numerical solution of the Euler equations. They show that large-amplitude mode-2 waves [with , where is the amplitude and *H* is the water depth] lose coherence (wave form) during overtaking and colliding interactions with the mode-1 NLIWs. Stastna et al. (2015) suggests this degeneration is due to the locally modified stratification and vertical shear seen by the mode-2 wave during the interaction phase.

Several generation mechanisms for mode-2 NLIWs have been proposed in the literature (see the review by Yang et al. 2009). Many of these mechanisms involve mode-1 NLIWs interacting with topography and either scattering or breaking. Another mode-2 generation mechanism is the release of a density front into a pycnocline, and this mechanism is often used to generate mode-2 waves in the laboratory and direct numerical studies (e.g., Brandt and Shipley 2014; Deepwell and Stastna 2016). Grisouard et al. (2011) proposed another mechanism, which they name local generation, where internal beams hit a sharp pycnocline and directly trigger a high-mode (mode-2 and greater) solitary wave train. This mechanism is dependent on the geometry of the beam and the pycnocline thickness.

The two main goals of this study are, first, to characterize the occurrence and variability of mode-2 (and mode-1) NLIWs on the continental shelf of Northern Australia and, second, to identify the relevant environmental conditions that led to large-amplitude waves forming. Additionally, we will characterize the wave properties of some of the larger-amplitude events. In section 2 we outline the weakly nonlinear wave theory and the concept of the nonlinear steepening length that will be used to interpret the observations. An overview of the study site, numerical modeling results, and the field experiment is given in section 3, while section 4 outlines the analysis techniques applied to the observations. In section 5, new observations of mode-1 and mode-2 NLIWs are presented along with details of the environmental parameters. In section 6 we discuss the relationship between the steepening length and large-amplitude internal wave occurrence in the region and how these results can be applied to other regions. Finally, a nonlinear model is used to demonstrate the mode-2 generation mechanism.

## 2. Weakly nonlinear theory

The evolution of a small-amplitude, unimodal disturbance in a continuously stratified fluid of depth *H* can be described by weakly nonlinear (WNL) theory with the Korteweg–de Vries (KdV) equation (cf. Lamb and Yan 1996):

where

and

are the nonlinear steepening and dispersive coefficients, respectively. The modal structure functions are dependent on the density stratification via the hydrostatic, shear-free Taylor–Goldstein equation (Gill 1982):

with boundary conditions and the normalization condition . The modal structure functions and the linear phase speed *c* are eigenvector and eigenvalue solutions to Eq. (4). Various extensions to WNL theory have been proposed that account for vertical shear, higher-order nonlinearity, variable topography, rotation, and dissipation [see the review papers by Apel et al. (2006) and Ostrovsky and Stepanyants (1989)]. We present the equations for *α*, *ϕ*, and *c* that include the effect of vertical shear in the appendix.

Under WNL theory, the wave amplitude function is related to the two-dimensional (2D) velocity and buoyancy fields by multiplying it with the vertical structure functions and adding nonlinear terms in an asymptotic fashion. Velocity fields are derived from continuity and using the 2D streamfunction *ψ*, where , , and , and the dots denote the additional nonlinear terms (i.e., terms containing higher powers of *A*; see, e.g., Lamb and Yan 1996; Holloway et al. 1999). The buoyancy perturbation and other quantities such as isopycnal displacement and horizontal shear are derived in a similar manner. In section 4, we will use a method to extract *A* from the buoyancy perturbation derived from moored temperature observations.

The propagation speed of a nonlinear wave or wave train will vary depending on its form. For example, an analytical expression for the nonlinear phase speed , which is amplitude dependent, exists only for the special case of a finite-amplitude solitary wave solution of permanent form, that is, when nonlinear steepening balances dispersion (Helfrich and Melville 2006). For small-amplitude waves that are evolving, a dispersion relationship is obtained from WNL theory by inserting plane wave solutions of the form into a linear version of Eq. (1), resulting in (cf. Ostrovsky and Stepanyants 1989)

where is the nonlinear phase speed. It is apparent from Eq. (5) that one must know both the wave frequency *ω* and the horizontal wavenumber *k* to estimate the nonlinear phase speed. Furthermore, the wave amplitude plus the balance between dispersion and nonlinearity also needs to be considered before applying Eq. (5) to observations.

There are various limitations of the WNL theory that should be noted. First, it is only applicable for small-amplitude perturbations, defined to be when . Second, interactions between modes are assumed small and hence ignored; that is, different wave modes are assumed to evolve independently. Last, the velocity predictions are often in poor agreement with fully nonlinear solutions (Lamb and Yan 1996). Fully nonlinear models do exist [e.g., the steady-state Dubreil–Jacotin–Long (DJL) equation] that do not make the small-amplitude assumption and have shown good agreement with mode-1 wave observations when describing individual steady-state solitary waves (e.g., Lien et al. 2014). The distinct advantage of the WNL theory, however, is that it includes time-dependent effects, thus making it useful for understanding nonlinear wave evolution in observations [see Helfrich and Melville (2006) for a thorough overview of the different nonlinear internal wave models and their limitations].

An internal tide—a sinusoidal wave with tidal frequency *ω* and horizontal wavenumber *k*—requires sufficient time, or space, to evolve from its initial sinusoidal form at generation into a nonlinear wave train. Following Horn et al. (2001), this steepening time can be estimated using WNL theory by comparing the unsteady and steepening terms in Eq. (1) and noting the maximum steepness is , where is the initial amplitude and is the horizontal wavenumber, to obtain

The steepening length scale for an internal tide can thus be written as

and indicates whether nonlinear wave steepening will occur within a given horizontal distance from the generation site.

## 3. Location and observations overview

### a. Field experiment overview

The field site was situated on a broad section of the northern portion of the Australian North West shelf known as the Browse Basin (Fig. 1). The moorings were deployed in 250-m-deep water, and the surrounding region was relatively flat with an average topographic slope of roughly 0.2% for at least a 40-km radius around the site. The barotropic tidal ellipse, calculated from the measurements, was oriented in a northwest (NW)–southeast (SE) direction and had a peak amplitude of approximately 0.40 m s^{−1}. Previously, Rayson et al. (2012) presented observations of significant mode-1 and mode-2 internal tides at a 550-m-deep mooring roughly 150 km southwest (SW) of the current site.

Three through-water-column moorings (NP250, WP250, and SP250) were deployed in a triangular configuration with a horizontal spacing of roughly 600 m (Fig. 1). The choice of horizontal spacing was based on initial estimates of a wavelength of 1.2 km for the mode-1 internal solitary wave using numerical solutions to the KdV equation [Eq. (1)] with representative environmental conditions. The moorings were deployed from 2 April to 22 May 2017, although reliable data were obtained for the first 5 weeks only. All three moorings were configured to sample water temperature [Seabird Electronics (SBE) 56 and 39 models] at roughly 10-m vertical spacing from the seabed to 20 m below the free surface. Velocity was measured using either 75-, 150-, or 300-kHz Teledyne RDI acoustic Doppler current profilers (ADCPs) configured to sample at 4–4.5-m vertical resolution and to cover most of the water column. Additional pressure sensors (SBE39-TP) were attached to each mooring to measure knockdown. Four conductivity sensors (SBE37) were spaced through the water column on NP250 to measure salinity. Temporal sampling resolutions were 0.5, 10, and 25 s for the SBE 56, 39, and 37 instruments, respectively. ADCPs were all configured to sample at 60 s and collected data in beam coordinates. Additionally, turbulence packages were mounted onto SP250, and a bottom-mounted frame instrumented to measure boundary layer turbulence was deployed during the experiment (SP250-Lander), but these measurements are discussed elsewhere. See Fig. 2 for details of the instrument layouts and sampling duration for each mooring.

### b. Regional internal tide

A three-dimensional hydrodynamic model, the Stanford Unstructured Nonhydrostatic Terrain-Following Adaptive Navier–Stokes Solver (SUNTANS; Fringer et al. 2006), was run in hydrostatic mode to provide insight into the regional internal tide dynamics. We used an unstructured mesh and forced with realistic tides and representative density stratification for April [see Rayson et al. (2018) for details of the model configuration and validation]. The model was run for the period 23 March–10 April, and model variables from the spring tide period, 1–3 April, were harmonically decomposed into M_{2}, M_{4}, and M_{6} tidal frequency components. We then computed internal tide energy quantities in harmonic form to extract the mean energy density and energy flux magnitude and direction (see Rayson et al. 2018). The minimum horizontal resolution was 1 km, and the model solved the hydrostatic momentum equations only. In section 6, we use a 2D, 50-m resolution, nonhydrostatic scenario to investigate the nonlinear internal wave generation and evolution processes over a representative shelf region.

The numerical solution identified several steep topographic features on the shelf that generated large-amplitude internal tides (Fig. 3a). To the south of the measurement site, a northward-facing bank rapidly drops from 100- to 200-m water depth, emitting northward-propagating internal waves. The continental slope edge, indicated by the 400-m isobath, is approximately 75 km NW of the site and was also an internal tide generation region owing to the incident angle of the barotropic tide on the local slope. There is thus a complicated three-dimensional internal tide climatology (standing wave-like pattern) resulting from the various generation sites within the region (Rayson et al. 2011) making it unclear a priori from which direction nonlinear wave trains would arrive at the mooring site. Last, note that a vertical slice of the vertical velocity divergence field from the numerical model indicated that high-mode internal tides (internal beams) were generated at the offshore shelf break (Fig. 3b). In section 6, we apply the same code, but with the nonhydrostatic pressure solver, to a two-dimensional idealized cross section to demonstrate a mechanism for high-mode nonlinear wave generation at a shelfbreak region.

## 4. Analysis techniques

### a. Stratification

Calculating the environmental parameters *α* and *β* in Eq. (1) requires full-water column background density information. Moored temperature data were first converted to density using a nonlinear equation of state assuming a constant salinity (34.6 psu). Four CTDs deployed on NP250 showed that salinity was 34.6 ± 0.1 psu at all depths over the duration of the mooring deployment. Density data were then mapped from the moving pressure coordinates to fixed vertical levels via spline interpolation, where the average depth of each instrument was chosen as a representative fixed depth. We then fitted an analytical profile to the background density profiles to give smooth, monotonic, full-water column background density profiles.

A parametric double-pycnocline model was used to represent the density profile

where is the density difference, is the thickness, and is the midpoint of each pycnocline. Here the vertical coordinate *z* sign convention is zero at the surface and positive going up; that is, all depths are negative. A robust least squares routine was used to determine the seven free parameters in Eq. (8) that best fit a given measured density profile.

We used a backward-in-time low-pass filter derived from a finite impulse response differential equation to compute the background density

where is the filter time scale, set to 34 h, and *ρ* and are the instantaneous and background densities, respectively. The background density was computed at each observational time step *i* via

where is the sampling interval (60 s). Missing values were handled by substituting the previous (in time) good observation at . This filter was used because it utilizes the density for all times leading up to the arrival of a large-amplitude wave at the site.

### b. Modal amplitude fitting

The wave amplitude was extracted from the mooring data by least squares fitting both the velocity and density data to modal structure functions, similar to the procedure described by Rayson et al. (2012). Note that the weakly nonlinear theory only allows expansion into a single mode and does not account for interactions between modes. Furthermore, nonlinear forms of the velocity and buoyancy functions are not linearly separable, meaning only a single mode can be extracted from field measurements via fitting under the weakly nonlinear assumption. We fit the different vertical mode amplitudes to the buoyancy (temperature) observations via

where *n* is the mode number. Fitting nonlinear waves using only the superposition of linear modes is fairly common practice in both modeling (e.g., Yuan et al. 2018) and observational studies (e.g., Colosi et al. 2018).

In a unidirectional wave field with waves propagating in the *x* direction, horizontal and vertical velocity components are then given by (ignoring rotation)

meaning that, in theory, one could calculate horizontal and vertical velocity from the buoyancy perturbation amplitude alone. In reality, and in particular in the Browse Basin, the internal waves are not unidirectional so Eq. (13) cannot be used. We primarily use the buoyancy amplitudes calculated from Eq. (11) for internal wave characterization and analysis throughout this paper.

### c. Internal tide identification

While an internal tide is broadly defined as an internal wave with tidal frequency, in practice this definition does not capture the fact the tide is composed of multiple frequencies and also that internal tide phases are modulated (Doppler shifted) by the background flow and temporally evolving stratification. Furthermore, nonlinear steepening causes the internal tide to evolve into higher harmonics prior to solitary wave train formation.

To capture the nonstationary effects of steepening and Doppler shifting in the interpretation of the internal tide from our observations, we adopt the approach of least squares fitting a semidiurnal frequency plus its first two harmonics (i.e., M_{2}, M_{4}, and M_{6}) to 2-day segments of each amplitude record with a 50% overlap. We define the amplitude as the M_{2} component of the short-time harmonic fit (STHF), which is dominated by linear waves (see below).

### d. Wave speed and direction detection

The nonlinear internal wave speed and direction were derived from the triangular mooring array by using the lagged correlation of the modal buoyancy amplitude [Eq. (11)] signal between each pair of sites. Lag time , defined as the time when correlation peaked, was then used to form the system of equations

where and are the wave velocity components, is the distance between site pairs *i* and *j*, and *ε* is an error term. The triangular mooring configuration yielded a system of three nonlinear equations with two unknowns and . The fitted phase speed and direction were then computed from these vector components. We used an optimization routine to minimize the error term in the nonlinear system of equations given by Eq. (14). This wave direction detection method was applied to 6-h segment amplitudes and only large-amplitude mode-1 waves where and the amplitude m were retained. For mode-2 waves, only events where and m were retained for nonlinear wave direction analysis. The nonlinear phase speed was then directly estimated by removing the background barotropic velocity; that is, .

## 5. Results

### a. Internal tide amplitude

The mode-1 internal tide amplitude , defined here as the semidiurnal component of the short-time harmonic fit to the buoyancy amplitude, had a maximum of 30 m and tended to follow the spring–neap cycle of the cross-shelf barotropic tidal current (Figs. 4a,b). The largest mode-1 internal tides occurred in early April when the maximum barotropic forcing was 0.3 m s^{−1}. A lagged correlation between the RMS tidal amplitude and mode-1 harmonic amplitude revealed a 54-h lag between peak spring tide and peak mode-1 internal tide amplitude with *r*^{2} = 0.64 (not shown). Despite the barotropic velocity being 25% larger in May, the internal tide amplitude was roughly 30% smaller. This suggests another process, likely wave–wave interactions, was contributing to the greater mode-1 internal tide amplitude during early April and thus accounts for some of the unexplained variance in the relationship with the surface tide.

The mode-2 internal tide amplitude showed little correlation with the barotropic tidal forcing. The short-time harmonic fit to the buoyancy amplitude derived from NP250 captured 84% of the mode-1 variance but only 36% of the mode-2 variance. Stationary tidal harmonic fitting using eight tidal constituents (M_{2}, S_{2}, N_{2}, K_{2}, K_{1}, O_{1}, P_{1}, Q_{1}) applied to the full time record (34 days) captured 71% of the mode-1 but only 18% of the mode-two variance (see Fig. 5 for a visual depiction of the fit quality for both mode-1 and mode-2 amplitude signals). This suggests that the mode-1 response was predominantly due to linear internal tides, whereas the mode-2 response was predominantly due to NLIWs.

### b. Mode-2 NLIWs

Large-amplitude mode-2 NLIWs were seen throughout the observation period with amplitudes of 30–50 m—similar to the mode-1 internal tide amplitude (Figs. 4c,e). Convex nonlinear mode-2 waves occurred throughout the 5-week observation record (three examples are shown in Fig. 6). Note the mode-2 wave amplitude at WP250 should be treated with caution owing to a gap in thermistors between 30 and 190 m above sea bed (ASB).

The first example of a large mode-2 NLIW event (18 April 2017) had a 20–30-m amplitude leading wave with a weak trailing wave signal. Amplitudes of the leading wave were roughly equal at all moorings for this event. The second example (30 April) was characterized by a smaller leading wave and larger-amplitude trailing waves at both NP250 and SP250. Visually, this event had the greatest spatial variability with a 45-m amplitude disturbance observed at NP250 and only a 25-m maximum at SP250. Isotherm overturning was observed in the core of the large-amplitude trailing wave at both NP250 and SP250. The largest-amplitude mode-2 NLIW was observed on 6 May 2017 with an amplitude of 50 m (Fig. 6l). This event had a large leading wave followed by three to four smaller-amplitude trailing waves. This event passed NP250 and WP250 first and then reached SP250 roughly 15 min later.

The 6 May event had the largest isotherm displacement amplitude and also some of the strongest horizontal velocities observed at the site with *u* = 0.80 m s^{−1} measured at SP250 (Fig. 7a). Vertical velocities during this event were up to 0.25 m s^{−1} (Fig. 7b). High-frequency (>*N*) temperature oscillations, detected with the SBE56 thermistors sampling at 2 Hz, were also evident in the core of this mode-2 wave, as well as isotherm overturning, for example, at 0747 UTC around a depth of 140 m. Other cases of isotherm overturning in the core of mode-2 waves were observed during smaller-amplitude waves as well, for example, on 20 April 2017 (Fig. 8). In both cases, overturning occurred in the rear core of the wave and super-*N* frequency isotherm oscillations persisted for about 10 min after the passage of the core at 120–150-m depth within the pycnocline. These high-frequency oscillations indicate a midwater column region of enhanced turbulent dissipation and vertical mixing likely driven by internal wave breaking.

### c. Mode-1 NLIWs

A small number of large-amplitude mode-1 waves of depression were observed between 2 and 5 April 2017. The three largest events were characterized by an initial depression with amplitude of 50–70 m followed by smaller-amplitude trailing waves with a period of 10–15 min (Fig. 9). The period (trough to trough) of the leading wave was 20–30 min, approximately twice that of the trailing waves for all three cases. The cleanest example of this type of behavior was observed at SP250 on 2 April (Fig. 9g). Each event resulted in a similar buoyancy amplitude at all three moorings. The NLIWs always reached SP250 last (Figs. 9j–l), suggesting the waves were propagating toward the SE. Nonlinear mode-1 waves of depression or elevation were not observed outside of the short time period in early April. Usually, the mode-1 buoyancy amplitude was well characterized by the short-time harmonic fit; that is, it was usually linear. Nonlinear mode-1 waves of elevation were not observed during the observation period.

### d. Environmental (KdV) parameters

The moorings were deployed during the summer to winter (monsoon to dry season) transition period, and the vertical density profile transitioned from having a double thermocline structure at 50 and 175 m around the beginning of April, to a deeper and sharper single thermocline by 30 April (Fig. 10a and Fig. 11a). The corresponding buoyancy frequency profile had a double peak around 3 April (*N* = 0.015 s^{−1}) and a single peak around 30 April (*N* = 0.022 s^{−1}) or buoyancy periods of 7 and 5 min, respectively (Fig. 10b and Fig. 11a). Linear mode-1 phase speeds increased from 1.16 to 1.28 m s^{−1} owing to the sharpening thermocline, while the linear mode-2 phase speed decreased from 0.58 to 0.48 m s^{−1} (Fig. 11b).

Thermocline deepening toward late April 2017 resulted in deeper peaks for both the mode-1 and mode-2 vertical structure functions, and importantly this changed both the sign and magnitude of the nonlinear steepening parameter [according to Eq. (2)]. The depth of the mode-1 structure function peak deepened from slightly shallower than 125 m (middepth) to 140 m by the end of April, while the depth of the mode-2 structure function peak deepened from 60 to 80 m over the corresponding period (Figs. 10c and 10d, respectively). Over the entire period, the mode-1 nonlinearity parameter *α* thus changed sign from (conducive to waves of depression) to m^{−1} (conducive to waves of elevation). Note also there was roughly a two-week period beginning 5 April when the mode-1 *α* was close to 0. The mode-2 nonlinearity parameter *α* was initially negative at but was positive over most of the period and increased in magnitude to m^{−1} at the end. Importantly, the mode-2 parameter was typically 2–5 times larger in magnitude than the mode-1 parameter throughout the period, meaning mode-2 waves could steepen faster than mode-1 waves of equivalent initial amplitude.

### e. Background shear effect

The influence of background shear on the environmental parameters was assessed using moored ADCP velocity measurements and computing the first-order nonlinear coefficient with shear (see appendix). ADCP measurements from the SP250 75-kHz RDI Long Ranger were low-pass filtered, similar to density, and then rotated into a plane along the approximate direction of wave propagation (−50° from east as shown below). Velocity data were smoothed and suitably extrapolated near the boundaries by fitting fourth-order Chebychev polynomials and mapping onto the same vertical grid as . The background velocity peaked at 0.30 m s^{−1}, although it was generally from −0.15 to 0.15 m s^{−1} (Fig. 12a). The velocity in the thermocline was positive during April, aligned with the direction of wave propagation, and only reversed direction in early May.

The mode-1 and mode-2 nonlinear coefficients ( and , respectively, as defined in the appendix) were calculated with the WNL theory using the background velocity rotated onto the approximate plane of wave propagation (Fig. 12). For mode-1 the range of with shear was similar to the range of no-shear values over most of the observation period but was sometimes up to 0.03 s^{−1} less than the no-shear values. For mode-2 the range of with shear was again similar to the range of no-shear values, but the differences between the two were greater. Toward the end of the record, for example, was 0.05 s^{−1} less than the no-shear case.

### f. Nonlinear wave speed and direction

Using the lagged correlation to determine the arrival time at each of the three moorings (section 4d), we were able to estimate the phase speed and direction of the largest amplitude mode-1 and mode-2 NLIWs (Table 1; see also Figs. 9 and 6j–l as examples of the wave amplitude signals used for the fitting). The nonlinear internal waves observed at the mooring were primarily directed toward the SE quadrant, indicating their origin was at the shelf break to the NW of the site. Large-amplitude mode-1 waves of depression mainly occurred on ebb tide and all occurred when s^{−1}. Large-amplitude mode-2 waves occurred more intermittently but also frequently throughout the record.

The deviation between the nonlinear and linear phase speed , where , varied from −0.13 to 0.14 m s^{−1} and from 0.07 to 0.36 m s^{−1} for the mode-1 and mode-2 NLIWs, respectively (Table 1). The WNL dispersion relationship [Eq. (5)] predicts that the nonlinear contribution results in larger phase speeds, that is, , given the positive values of the dispersion term for both modes ( and m^{3} s^{−1}). The negative was estimated for two mode-1 waves (1947 UTC 3 April 2017 and 2049 UTC 4 April 2017), indicating a deviation from WNL theory or a large error in the optimization to get . Including shear in the linear phase speed calculation (see the appendix) does not account for this large deficit. Shear actually increases the linear phase speed from 1.2 to 1.4 m s^{−1} during this time period. Estimation of the representative barotropic tide and background velocity add additional uncertainty to the calculation. It is beyond the present scope to explain this deviation from WNL theory, but we note this highlights some of the challenges with directly estimating NLIW propagation speed in the ocean.

Wave steepness and Froude number are used as diagnostics to determine the stability of large-amplitude internal waves (Zhang and Alford 2015; Thorpe 2018). The nonlinear wavelength was determined from our estimate of via Eq. (5). Mode-1 NLIWs were 1030–1870 m long, while mode-2 waves were 220–550 m. The wavelength of the large-amplitude mode-2 NLIWs on 6 May 2017 (Fig. 7) was 330 m resulting in a wave steepness of 0.15, while the event on 20 April 2017 (Fig. 8) was 0.10. The wave Froude number was greater than unity for both of these waves, suggesting they may have been undergoing overturning via a kinematic instability. The maximum steepness of mode-1 waves was 0.07 and no overturning events were observed in these waves (Table 1). Zhang and Alford (2015) observed increased steepness when mode-1 waves were undergoing kinematic instability in agreement with our observations of large steepness (0.15) during unstable mode-2 waves.

An important implication of the finding that waves were directed in a SE direction is that the internal tides generated at the shelf break south-southwest of the site did not evolve into detectable nonlinear wave trains by the time they reached the observation site (see Fig. 3). Both mode-1 and mode-2 nonlinear internal waves thus originated from a direction that was different from the primary time-averaged energy propagation direction as predicted by the 3D hydrostatic numerical model. We argue that the differences in steepening length scale explain this observation.

## 6. Discussion

### a. Steepening length scale calculation

A goal of this study was to identify the background conditions that promote an increased likelihood of occurrence of large-amplitude NLIWs. Key to this is determining the steepening length scale, influenced by the initial wave amplitude and wavenumber, phase speed, and steepening parameter, and these latter quantities are strongly dependent on the background density stratification in the domain. We assumed is represented by the amplitude of the semidiurnal (M_{2}) component of the short-time harmonic fit to either the mode-1 or -2 amplitude signal (see Fig. 4). The maximum amplitude was calculated by finding either the maximum or minimum modal amplitude ( or ) in 6-h blocks. We retained the sign to determine the polarity of the waves. For particular wave events we used the representative density stratification at that time to compute the steepening parameter and wave speed in the domain to compute from Eq. (7).

Mode-2 waves typically had a shorter steepening length scale than the mode-1 waves. Mode-2 waves followed an inverse relationship with the largest-amplitude NLIWs ( m) having steepening length scales km, a distance shorter than the length from the offshore shelfbreak generation site to the observation site (Fig. 13). While depends on both the initial amplitude and the nonlinear steepening parameter, the latter typically had the greatest influence. For example, the largest mode-2 wave occurred around 6 May when was largest even though was relatively small (5 m; is shown by color in Fig. 13).

Mode-1 waves of depression typically had a longer steepening length scale than mode-2 waves with only the very largest waves having km. This long mode-1 steepening length scale was likely the reason why nonlinear mode-1 waves were rarely observed. Assuming mode-2 waves originated at the outer shelf break as a sinusoidal plane wave (as indicated by Fig. 3b and discussed below), then the relatively short steepening time (and length scale) would allow them to evolve and lose coherence more rapidly than mode-1 waves. This is also likely why only a small fraction of the mode-2 variance was captured via harmonic fitting: 36% for mode-2 versus 84% for mode-1. Continued steepening could cause these waves to become unstable and eventually break. The lack of mode-2 observations in the global ocean has been attributed to this highly dissipative character, and the observations shown here are in a region where mode-2 waves underwent the first stage of this evolution (nonlinear steepening) before significant breaking and dissipation occurred.

While shear effects could potentially be included in the determination of , by using as defined in the appendix (see Fig. 12), we have chosen not to do this for several practical reasons. First, choosing a representative steady background flow that acts uniformly over the entire domain through which the wave evolves requires velocity measurements across the width of the shelf (100–150 km in our case). A shelf-scale eddy, for example, would cause velocity variations along the wavelength. Second, one must know the exact wave propagation direction a priori in order to rotate in that plane. Third, extracting a suitable continuous profile that spans the entire water column is also challenging given that long-range ADCPs suffer from acoustic interference close to boundaries (approximately the upper 20 m in our case). We used Chebychev polynomials to smooth the profile and give a reasonable extrapolation. Numerical solutions to and *c* using Eq. (A4), however, were found to be sensitive to this extrapolation. For these reasons we restrict our attention to a first-order predictor for large-amplitude NLIW occurrence.

### b. Spatial variability of the steepening length

We used the 3D SUNTANS numerical model solution to examine the spatial variability of the steepening length scale. The internal tide amplitude was calculated by first fitting harmonic amplitudes of density (buoyancy) perturbations, as described in detail in Rayson et al. (2018) and discussed in section 3b. The modal eigenvalue problem [Eq. (4)] was then solved at each horizontal grid cell using the mean density from the April simulation. Last, modal amplitudes were fit via least squares to the complex harmonic amplitudes using Eq. (11).

The largest predicted mode-1 internal tide had an initial amplitude of m near the 200-m shelf break, roughly 50 km south of the field site (Fig. 14a). There was a standing wave pattern over much of the shelf in depths 200–500 m due to the interaction of waves from multiple generation sites (Rayson et al. 2011, 2012). Mode-2 internal tides also exhibited a standing wave pattern, although the wave lengths were shorter and the maximum amplitudes were roughly 50% smaller with m (Fig. 14b).

In the 200–500-m shelf region, the mode-1 steepening length scale was 50–300 km, whereas for mode-2, was in the range 20–50 km (Figs. 14c,d). This indicates that the shelf region, between 200- and 500-m water depth, corresponding to a region of about 100 km in width, was most conducive to nonlinear mode-2 internal wave formation. This spatial trend is supported by the temporal analysis of the mooring data where nonlinear waves accounted for more of the mode-2 buoyancy amplitude signal variance, as discussed in section 5a. Seasonal variations in stratification should also be considered, although our results are focused on the April observation period.

### c. Fully nonlinear numerical simulations

We solved the nonhydrostatic Reynolds-averaged Navier–Stokes (RANS) equations over an idealized, 2D shelfbreak topography using the SUNTANS code to examine the evolution of the internal tide beam into mode-2 NLIWs. Early May 2017 stratification conditions (where the thermocline was around middepth) were used. The shelf-slope region was idealized as a 2D (*x*–*z*) section where the bottom depth is given by

where *H* is the offshore water depth, is the shelf height, is the location of the shelf break, and is the shelf width. Representative values of the topography for the Browse Basin region were selected to ensure a supercritical shelf break and a 250-m-deep shelf (see Table 2).We ran the numerical solution with stratification conditions representative of 30 April 2017 and forced the boundary with an oscillatory barotropic flux of 110 m^{2} s^{−1} and M_{2} frequency (1.42 × 10^{−4} s^{−1}) that was representative of spring tide conditions (see Fig. 10 and Table 2). The nonlinear coefficients indicated a greater tendency for mode-2 steepening according to WNL theory, with values of 0.0017 and 0.0089 s^{−1} for mode-1 and mode-2, respectively.

The model time step was 5 s, while the horizontal and vertical grid resolutions were 50 and 10 m, respectively. The horizontal resolution satisfied the criteria for the model nonhydrostatic dispersion of nonlinear internal waves to be greater than numerical grid-induced dispersion as outlined in Vitousek and Fringer (2011). A free-slip bottom boundary condition was applied, while sponge conditions were applied to the horizontal momentum equation at the lateral boundaries to absorb baroclinic motions. The horizontal and vertical eddy viscosity were set to 1.0 and 1 × 10^{−4} m^{2} s^{−1}, respectively. The simulation was run for 10 tidal cycles.

Mode-2 NLIWs evolved in the isopycnal field over a tidal cycle across the shelf region (Fig. 15). A virtual mooring at a distance of 110 km shows a time series of isopycnal displacement with a mode-2 NLIW evident at the start and toward the end of the tidal cycle (Fig. 15a; time = 5.5 and 17.5 h, respectively). This wave train evolved from a linear mode-2 wave that was generated near the shelf break and was most pronounced at a distance of 83 km (Fig. 15b; mode-2 waves are characterized by a middepth peak in the vertical velocity divergence ). A convex mode-2 shock-like feature then formed around a distance of 92 km (Fig. 15e). This shock-like feature evolved into a wave train, reminiscent of the observations (Fig. 6), toward the end of the tidal cycle (Figs. 15f,g). An internal beam emanated from the shelf break and was most pronounced in Figs. 15d–f and, as discussed below, is the source of NLIWs on the shelf.

This fully nonlinear numerical solution, with its relatively simple forcing and topography, demonstrates not only the generation of a mode-2 NLIW but also the complex interactions between different modes in the ocean. The shape of the mode-2 shock stretched and contracted as a mode-1 internal wave passed through it (Figs. 15e–g) (mode-1 waves are characterized by reversing from top to bottom). This mode-1 internal wave subsequently steepened into a nonlinear wave train of elevation at a distance of 120 km (Fig. 15f). This wave of elevation had not evolved into a wave train at a distance of 110 km owing to the longer mode-1 steepening length and was therefore not observed at the virtual mooring unlike the mode-2 NLIW. This numerical example demonstrates how the combined effects of steepening length and observation location both affect the observed wave form at a fixed point.

### d. Necessary conditions for mode-2 NLIWs

Grisouard et al. (2011) make a geometric argument stating that for a particular vertical mode to be excited, the phase lag between a beam interacting with the bottom and then the top of a pycnocline has to be equal to the travel time for that given mode. For example, for mode-2 excitation, the top and bottom isopycnal displacement lag must be one half-wavelength out of phase. The equation for determining the excited mode in a single pycnocline is (cf. Grisouard et al. 2011)

where is the wave characteristic slope in the pycnocline, is a representative pycnocline buoyancy frequency, is the pycnocline thickness, and is the planar wavelength. For late April conditions (see Fig. 10), s^{−1} and = 22 000 m, thus requiring a necessary pycnocline thickness for mode-2 excitation, in reasonable agreement with the observed thickness in Fig. 10a and Table 2. This mode-2 generation mechanism is also evident in the numerical solution between 70 and 90 km where the slope generated beam displaces the base of the pycnocline while the change in beam propagation through the pycnocline results in a convex, mode-2 bulge (Figs. 15b–g).

Note that Eq. (16) is a linear excitation mechanism and the steepening length determined from WNL theory is still necessary to determine whether a wave has enough time to steepen by the time it passes a given observation point. The possibilities for high-mode excitation increase with a double pycnocline like Eq. (8), suggesting that the stratification-dependent evolution, not just the initial generation, can be an important factor for determining mode-2 NLIW occurrence.

The only comparable ocean site where mode-2 waves have been reported in the literature is the continental shelf region of the northern South China Sea (SCS), which also has a relatively wide shelf section with intermediate water depths between 200 and 500 m (Yang et al. 2009, 2010). Yang et al. (2009) observed both convex and concave mode-2 waves at a mooring in 300 m deep water, with the former occurring more regularly. On the Australian North West shelf (NWS), we mainly observed convex waves during our 5-week experiment.

Subsequent modeling studies in the SCS using a three-layer stratification revealed that the middle layer had to be thicker than half the water depth to support concave waves (Yang et al. 2010). For continuous stratification, the nonlinear parameter must be negative to support concave mode-2 NLIWs (Guo and Chen 2012). Guo and Chen (2012) used fully nonlinear numerical simulations to show that the concave waves on the SCS slope arrived from offshore and only formed into convex wave trains, like those observed on the Australian NWS, once they reached a shallower shelf region. The sign of in the observations presented here was usually positive except at the very beginning of April (Fig. 11c). During our period of study, the Australian NWS was therefore likely supportive of convex mode-2 waves in water shallower than 250 m, and concave waves were more likely to occur in deeper waters although there was no evidence of mode-2 waves offshore in the 2D numerical solutions. In the SCS, mode-2 NLIWs form at the Luzon Ridge and propagate across the ocean basin. Similar local generation at tall, sharp topography also occurs at the Mascarene Ridge (da Silva et al. 2015). In contrast, we propose that mode-2 internal waves were generated, at the northern Australian NWS location, via a local beam–pycnocline interaction mechanism near the shelf break.

The last condition necessary for mode-2 waves to be identified is that they cannot have already decayed in amplitude due to turbulent dissipation (e.g., Shroyer et al. 2010) or due to wave–wave interactions (e.g., Stastna et al. 2015). Shroyer et al. (2010) reported a turbulent dissipation time scale of h on the New Jersey shelf in relatively shallow water (<100 m). The observations presented here indicate likely regions of isotherm overturning and enhanced turbulent activity (see Figs. 7 and 8) although we have not calculated any turbulence quantities. A key to future successful prediction of mode-2 wave occurrence and amplitude will be to better estimate the dissipation rates and associated time scales, such as in Boegman et al. (2003), in order to properly parameterize this process in both weakly or fully nonlinear models of NLIW evolution.

## Acknowledgments

This work was supported by the Australian Research Council Industrial Transformation Research Hub for Offshore Floating Facilities (IH140100012). Ship time was partly funded through the Indian Ocean Marine Research Centre partners: the University of Western Australia, Commonwealth Scientific and Industrial Research Organisation, and the Australian Institute of Marine Science (AIMS). We are grateful to the captain and crew of the R/V *Solander* plus John Luetchford and Simon Spagnol from AIMS for their assistance in the deployment and recovery of the moorings. Cynthia Bluteau and Andrew Zulberti designed the moorings, and we thank them plus other UWA staff and students for their efforts at sea. Observational data are archived on the UWA library research repository (https://doi.org/10.4225/23/5afbf8fc55ed1).

### APPENDIX

#### Background Shear Equations

In the presence of a mean background flow , the eigenvalue problem for internal wave modes reduces to the Taylor–Goldstein equation (see, e.g., Stastna and Lamb 2002):

with boundary conditions . The subscript *z* denotes a derivative in the vertical direction. The nonlinear parameter in the presence of shear is

where the superscript *s* denotes the inclusion of shear and

We use a shooting method to solve Eq. (A1) for a given profile of and (cf. Hazel 1972). The first step is to break up Eq. (A1) into a system of (two) first-order ordinary differential equations (ODEs) by first setting that leads to

The goal is to find values of and *c* so that the boundary condition is satisfied. Numerically, we use a root-finding routine to find and *c* such that by repeatedly integrating Eq. (A4). In practice, we use the root function with a Krylov solver in the scientific python optimization library (www.scipy.org) and the odeint function to solve each ODE iteration. The solution to the shear-free eigenvalue problem [Eq. (4)] is used as an initial guess for *c*. Note that this numerical method for solving Eq. (A4) does not converge when for all depths or when is large, that is, for nonsmooth velocity profiles.

## REFERENCES

*Atmosphere-Ocean Dynamics*. Academic Press, 662 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).