Simulations of nonlinear harmonic generation by an internal wave beam incident on a pycnocline

Internal wave beams generated by oceanic tidal flows propagate upward and interact with the increasing stratification found at the pycnocline. The nonlinear generation of harmonic modes by internal wave beams incident on a pycnocline has recently been demonstrated by laboratory experiments and numerical simulations. In these previous studies, the harmonic modes were trapped within the pycnocline because their frequencies exceeded that of the stratified fluid below. Here, two-dimensional numerical simulations are used to explore the effect of incidence angle on harmonic generation at a thin pycnocline. At incidence angles less than 30 degrees (typical of oceanic beams), the lowest harmonic mode freely radiates in the form of an internal wave beam rather than being trapped within the pycnocline. The results indicate that nonlinear refraction is the primary mechanism for harmonic generation at incidence angles exceeding 30 degrees, but that interaction of the incident and reflected beams is more important at smaller incidence angles. The simulations are compared to weakly nonlinear theory based on refraction at the pycnocline. The results yield good agreement for trapped harmonics, but weakly nonlinear theory substantially underpredicts the amplitude of the radiated harmonics.


Introduction
Internal waves in a continuously stratified fluid are ubiquitous in geophysical fluid dynamics.In the oceans, a broad spectrum of internal wave frequencies and wave numbers is observed (Garrett and Munk, 1979;Holloway, 1980;Garrett, 2003;Garrett and Kunze, 2007), bounded by the Coriolis and buoyancy frequencies.Energy is injected at the lowest frequencies, and nonlinear effects transfer energy to higher frequencies, establishing the observed spectrum.Oceanic internal wave beams are formed by tidal flow over topography (Bell, 1975;Lien and Gregg, 2001;Llewellyn-Smith and Young, 2003;Petrelis et al., 2006;Martin et al., 2006;Cole et al., 2009;Balmforth and Peacock, 2009;Echeverri and Peacock, 2010;Zhao et al., 2011;Gayen and Sarkar, 2011) and propagate upward toward the ocean surface, where they interact with the ocean pycnocline and mixed layer.Here, we investigate the nonlinear interaction of internal wave beams with a thin pycnocline, which results in the formation of harmonic modes.This effect may be relevant both to the evolution of internal wave beams and the transfer of energy between internal wave modes more generally.
There are many possible nonlinear effects which may be important for oceanic internal waves.Among known nonlinear effects is the excitation of harmonic modes which occurs during internal wave beam reflection from a solid surface (Peacock and Tabaei, 2005;Tabaei et al., 2005;Gerkema et al., 2006;Rodenborn et al., 2011;Gostiaux et al., 2013).The mechanism for harmonic generation is the nonlinear interaction between different wave numbers in the overlapping incident and reflected beams.A similar effect has been observed in numerical simulations during beam reflection from a free surface or a mixed layer (Zhou and Diamessis, 2013).Colliding beams with different frequencies can also exhibit nonlinear generation of internal wave beams with frequencies equal to the sum and difference of the original frequencies, as has been shown theoretically by Tabaei et al. (2005) and in the laboratory by Smith and Crockett (2014).The parametric subharmonic instability is a nonlinear instability leading to the creation of two subharmonic waves whose frequencies are about half the original wave frequency, and has Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
For a plane wave incident on a thin pycnocline (a thin, strongly stratified layer), Thorpe (1998) suggested that nonlinear effects would result in a harmonic mode with double the frequency and wave number of the incident wave.Recently, the nonlinear generation of harmonic modes by an internal wave beam incident on a pycnocline was demonstrated experimentally (Mercier et al., 2012;Wunsch and Brandt, 2012).Both experiments reported a series of harmonic modes, while Wunsch and Brandt also confirmed the approximate doubling of the incident wave number predicted by Thorpe (1998).Several recent numerical studies (Grisouard and Staquet, 2010;Grisouard et al., 2011;Gayen and Sarkar, 2013;Dossmann et al., 2013;Diamessis et al., 2014) have also demonstrated harmonic generation at a pycnocline, although not all of these studies were primarily concerned with harmonics.In most of these previous experimental and numerical results, the frequency of the observed harmonic modes exceeded the buoyancy frequency of the stratified fluid below the pycnocline.The harmonics were unable to propagate into the lower stratified layer and were therefore trapped within the pycnocline.The exception is a single experiment reported by Mercier et al. (2012) showing a radiated harmonic wave for an incidence angle of 15 degrees.Recent numerical simulations (Diamessis et al., 2014) with thicker pycnoclines also suggested that the harmonics formed due to nonlinear refraction through the rapidly increasing stratification at the pycnocline base, rather than due to overlap of incident and reflecting beams at the top of the pycnocline.The physical mechanism for harmonic generation might plausibly differ from that observed in previous studies of harmonics formed during reflection from a solid surface.Although harmonics have not been reported in ocean observations of internal wave beams, a recent study in the South China Sea demonstrated the presence of harmonics of the semi-diurnal tide in the vicinity of the pycnocline (Xie et al., 2013).
In the present work, fully nonlinear two-dimensional numerical simulations of internal wave beams incident on a thin pycnocline were performed for two incidence angles and various stratification profiles.One incidence angle is small enough to allow the lowest harmonic to re-radiate into the lower stratified layer.Hence both trapped and radiated harmonics are seen, depending on the angle of incidence.Rotation is not considered.
The outline of the present work is as follows: Sect. 2 reviews the weakly nonlinear theory for plane wave refraction at a thin pycnocline, and provides a straightforward extension to small incidence angles.Section 3 outlines the numerical method used.Section 4 presents results and comparisons with the weakly nonlinear theory of Sect. 2. Section 5 summarizes the conclusions and possible extensions for future work.

Weakly nonlinear internal wave refraction
No comprehensive theory of harmonic generation at a pycnocline presently exists.Treating the pycnocline as a thin interface in a two-layer fluid, Thorpe (1998) proposed that nonlinear effects due to the displacement of the interface by an incident plane wave results in the formation of a harmonic mode.A more general theoretical framework for nonlinear effects in two dimensions (Tabaei and Akylas, 2003;Tabaei et al., 2005) has recently been adapted to investigate the generation of trapped harmonics during refraction by the pycnocline (Diamessis et al., 2014).Here, this approach is generalized to include radiated harmonics.Although a number of simplifying assumptions are made, the analysis illustrates the basic physical processes involved during harmonic formation.An idealized stratification profile consisting of piece-wise linear segments is used.In this method, plane wave solutions for the primary frequency and lowest harmonic are constructed.Depending on the incident internal wave beam frequency ω and the lower layer buoyancy frequency N o , the lowest harmonic mode takes the form of either an interfacial wave trapped within the pycnocline (2ω > N o ) or a re-radiated plane wave (2ω < N o ).
In two-dimensional, inviscid, Boussinesq flow, internal wave motions are governed by Tabaei et al. (2005): where the stream function ψ(x, z, t) is the solenoidal (rotational) part of the velocity field, with the velocity components given by u = −ψ z and w = ψ x .The fluid density is ρ, the buoyancy frequency is N, and the gravitational acceleration is g.The density fluctuation δρ(x, z, t) is the difference between the local density and the mean background stratification profile.The Jacobian (J ) terms represent the nonlinear effects.In uniform stratification, nonlinear effects do not exist for single-frequency internal wave beam solutions of Eq. (1), as the Jacobian terms vanish (Tabaei et al., 2005).In strongly varying stratification, as is found at the ocean pycnocline, nonlinear refraction effects are possible.In the present weakly nonlinear analysis, a plane internal wave with frequency ω and horizontal wave number k (wavelength 2π/k) is incident on the pycnocline, and the nonlinear Jacobian terms calculated from this wave provide the forcing for a harmonic mode with frequency 2ω and wave number 2k.The wave stream function ψ and density perturbation δρ are expressed as where * denotes complex conjugation (so that the solutions are real).Here ψ 1 and δρ 1 represent the internal wave of frequency ω incident on the pycnocline, and ψ 2 and δρ 2 represent the harmonic correction to the solution, which is assumed to be small.Inserting Eq. (2) into Eq.( 1) and keeping only the lowest order terms proportional to exp(-iωt) yields the Taylor-Goldstein equation (Drazin and Reid, 1981) for the stream function of plane internal waves: The terms proportional to exp(iωt) give the complex conjugate of Eq. (3).Higher order terms (proportional to ψ 1 ψ 2 ) are assumed small and neglected in Eq. ( 3), so that the harmonic has no effect on the primary wave beam (the weakly nonlinear assumption).In the absence of any variation in the background stratification profile N, the nonlinear terms vanish and solutions which obey the Taylor-Goldstein equation (Eq. 3) are exact solutions of Eq. (1) (Tabaei and Akylas, 2005).In variable stratification, the nonlinear Jacobian terms are nonzero and oscillate in time and space with the harmonic frequency and wave number, as they are quadratic functions of the incident beam stream function.Combining these nonlinear terms with the harmonic correction terms in Eq. ( 2), which have the same dependence on kx − ωt, yields In constant stratification, Eq. ( 4) reduces to the Taylor-Goldstein equation (Eq. 3) for a plane internal wave with frequency 2ω and wave number 2k.The right-hand sides are the nonlinear source terms due to refraction of the primary internal wave ψ 1 in variable stratification.Equation (4) shows that internal wave refraction through any vertically varying stratification profile (non-zero dN 2 /dz) would be expected to generate harmonic modes with amplitude proportional to the gradient of N 2 .For the specific case of a thin pycnocline below a mixed layer with an incidence angle greater than 30 • , Diamessis et al. (2014) found solutions to Eq. ( 4) using piecewise constant N, with matching conditions at the boundaries between each piecewise constant region (Drazin and Reid, 1981).Here this solution is generalized to include smaller incidence angles for comparison with numerical simulations.The pycnocline consists of a sharp density change ρ located at depth z = 0, with constant buoyancy frequency N o below the pycnocline (z < 0) and zero stratification in the mixed layer above the pycnocline (z > 0).The stratified and mixed layers are both assumed to be semi-infinite.
A plane internal wave with frequency ω < N o and horizontal wave number k is assumed to be incident on the pycnocline from below.The stream function for the primary mode is given by Here A 1 is the amplitude of the incident wave and θ is its incidence angle.Solutions for the reflected wave amplitude B 1 and evanescent wave amplitude C 1 have been computed by Mathur and Peacock (2009) for ρ = 0 and generalized to include a density jump by Wunsch and Keller (2013) and Diamessis et al. (2014).
The solution for the harmonic mode ψ 2 with wave number 2k and frequency 2ω for the case 2ω > N o (incidence angles θ > 30 • ) was found by Diamessis et al. (2014): This solution corresponds to an interfacial wave with amplitude B 2 confined to the pycnocline.As is typical for weakly nonlinear theories, the harmonic mode amplitude B 2 is proportional to the square of the incident wave amplitude A 1 , which is assumed to be small.The constant of proportionality is a function of the parameter γ , which gives the dimensionless strength of the pycnocline in terms of the incident wave number and the underlying buoyancy frequency.A similar parameter was proposed to govern the coupling between incident internal waves and interfacial waves in Akylas et al. ( 2007) and Mercier et al. (2012).The amplitude B 2 is infinite (a singularity) at γ = 2 (1+α) sin 2 θ, where the denominator of Eq. ( 6b) vanishes.The singularity occurs when the harmonic mode lies on the dispersion curve for freely propagating interfacial waves on the pycnocline.Thorpe (1998) found a slightly different expression for the harmonic amplitude, but it also had a singularity at the same value of γ .Grisouard et al. (2011) offers a qualitatively similar explanation based on phase-speed matching for trapped pycnocline modes.Very close to the singularity, weakly nonlinear theory predicts arbitrarily large harmonic amplitudes, but its assumptions clearly fail in this case due to wave-breaking or other strongly nonlinear effects.In the simulations of Diamessis et al. ( 2014), the largest harmonic amplitudes were found for incident waves with values of γ in the vicinity of the singularity.In this work, results for numerical simulations with incidence angles θ = 40 • will be compared to Eq. ( 6).
For incidence angles θ less than 30 • (2ω < N o ), the solution is an internal wave which propagates into the lower stratified layer: This case has not been explored in previous numerical or laboratory studies.There is no singularity in the amplitude of the harmonic mode, which is always finite.In this work, results for numerical simulations with incidence angles θ = 25 • will be compared to the predictions of Eq. ( 7).

Numerical simulation method
Fully nonlinear two-dimensional numerical simulations were performed for two internal wave beam incidence angles and a variety of pycnocline characteristics.The two-dimensional equations of motion used in the simulations are Here x and z are the horizontal and vertical coordinates, respectively; u and w are the horizontal and vertical velocities; p is the pressure perturbation (deviation from hydrostatic); and ν and κ are the fluid kinematic viscosity and mass diffusivity, respectively.The total density ρ consists of three distinct components: a mean value ρ o , an imposed vertically varying stratification profile determined by N (z), and a density perturbation δρ which is advected according to Eq. (8c).These three density constituents obey the Boussinesq approximation, with the background variations due to N (z) being a few percent of the mean, and the fluctuations δρ/ρ o no greater than 10 −5 .Spatial discretization of Eq. ( 8) is carried out via a pseudospectral element method using Chebyshev polynomials, with time advancement executed via a third-order Runge-Kutta marching scheme (Ku et al., 1987(Ku et al., , 1989).Details of the numerical method can be found in the Appendix.Like all simulations, it frequently produces nonphysical oscillations when confronted with the steep density gradients that occur inside the pycnocline.To control such oscillations, Engquist et al. (1989) devised a nonlinear filtering algorithm for equally spaced grid points so that the desired portions (low frequency) are retained while the undesired (oscillatory high frequency) portions are eliminated.It locates local extrema and adjusts adjacent points, subject to the constraint that the summation of each field variable remains the same before and after the filtering.A modified version of the Engquist et al. (1989) filtering algorithm was applied here.The modification not only allows the scheme to be applied for unequally spaced allocation points, but also enforces a total variation diminished (TVD) scheme.The TVD-enforcing filter suppresses small-amplitude spurious oscillations near sharp boundaries and has been successfully implemented in Ku and Sibeck (1997).Here, the nonlinear filter was applied once every 60 time steps to the perturbed density δρ and velocity components u, w along a rectangular stripe enclosing the thin pycnocline to eliminate any nonphysical oscillations, if they exist.
The background buoyancy frequency profile (Eq.1d) used here is given by where h is the thickness of the pycnocline and λ is the internal wave beam width.The pycnocline has a buoyancy frequency that is larger than that of the lower stratified layer by a factor of r.Above it is an unstratified mixed layer, with thickness λ − h, sufficient to avoid any finite depth effects on pycnocline oscillations.At the pycnocline boundaries, the buoyancy frequency changes linearly over a short distance (0.1h) so that N(z) is continuous at the grid scale.The density change across the pycnocline is found by integrating Eq. ( 9): Figure 1  to resolve the large density gradient in the pycnocline, there are eight elements within the pycnocline and two more on each pycnocline boundary.The horizontal resolution is uniform except in the vicinity of the beam source, where the resolution is more refined.Simulations were run for 20 internal wave beam forcing periods to reach a steady state.
The internal wave beam is generated by imposing the oscillating boundary condition at the left (x = 0) boundary, within the range −4λ < z < −3λ.This boundary condition is intended to mimic the oscillating plate wave source with frequency ω f and amplitude A f which is commonly used in laboratory experiments (Gostiaux et al., 2007;Mercier et al., 2012).The source velocity is A f ω f at the center of the source aperture (z = −3.5λ)and is reduced quadratically toward the edge.Elsewhere on the side and bottom boundaries, no slip boundary conditions are applied.Free slip boundary conditions are applied to the upper boundary.
At the bottom and right edges of the simulation domain, Rayleigh absorbing layers (Klemp and Durran, 1983) prevent internal wave reflections from the boundaries.Within these absorbing layers, an artificial damping term is added to the right-hand side of equations of motion and mass: Equation (12a) applies in the right damping layer (x d < x < x L ), and Eq.(12b) applies in the lower damping layer (z L < z < z d ).Here β = 28 A f ω f /λ is a damping coefficient.The right damping layer thickness is 1.2λ, and the lower damping layer thickness is 3.3λ.Simulation parameters were selected to represent laboratory-scale flows, such as those reported in Wunsch and Brandt (2012) and Mercier et al. (2012).Parameters common to all simulations are λ = 10 cm, h = 0.5 cm, N o = 0.7 rad s −1 , A f = 0.25 cm, ν = 10 −2 cm 2 s −1 , κ = 10 −4 cm 2 s −1 , and ρ = 1 g cm −3 at the top of the domain.The only parameters which vary are the beam incidence angle θ and the pycnocline buoyancy frequency ratio r.These are listed in Table 1.The incident beam contains a continuum of horizontal wave numbers, but the wave number where the kinetic energy spectrum peaks is listed as k o in Table 1.The corresponding pycnocline parameter (Eq.6d) for this wave number is also given, although the incident beam consists of a spectrum of wave numbers, each associated with a different value of γ (k).r = 5.To extract the incident beam amplitude and wave number spectrum, a cross section of the beam was extracted at a depth z = −27 cm.An example power spectrum of the kinetic energy in the incident beam as a function of the crossbeam wave number k η and the frequency ω is presented in Fig. 3 (left).The spectrum is nearly monochromatic in frequency but has significant power for wave numbers up to ∼ 1.5 cm −1 .The absence of harmonic content in the incident beam is important, since harmonic generation is the focus of the present study.Wave number spectra filtered to isolate the frequency ω f and its lowest harmonics are presented in Fig. 3 (right).These incident spectra are nearly identical for all simulations.The crossbeam wave number k η where the forcing frequency spectrum peaks (0.65 cm −1 ) is converted to a horizontal wave number using the relation k = k η sinθ to compute the peak wave number k o listed in Table 1.Using the observed horizontal velocity and peak wave number, the wave steepness (Diamessis et al., 2014), defined here as max(uk η /2πN o ), was 0.01 for all simulations.The numerical simulations have been designed to explore the thin pycnocline configuration described in the weakly nonlinear refraction section.The piece-wise linear buoyancy frequency profile is less realistic than the continually varying profiles of Grisouard et al. (2011) and Diamessis et al. (2014), but it allows for precise definitions of the pycnocline thickness h and density change ρ.To satisfy the assumptions of the previous weakly nonlinear calculation, the pycnocline characteristics must satisfy the conditions k o h 1 and r 2 1 (Diamessis et al., 2014).In the simulations, the dimensionless thickness of the pycnocline k o h is 0.13 and 0.21 for the θ = 25 • and θ = 40 • cases, respectively.With the possible exception of the smallest r values in Table 1, the simulations presented here should satisfy the thin pycnocline assumptions.

Results
Simulations were performed to elucidate how harmonic generation varies with beam incidence angle and pycnocline characteristics for a thin pycnocline.Two harmonic generation mechanisms were expected to contribute: refraction through the pycnocline, as described in Sect. 2 and in the previous simulations of Diamessis et al. (2014), and overlap of the incident and reflected beams below the pycnocline, as in the harmonic generation at solid surfaces studied by Peacock and Tabaei (2005), Tabaei et al. (2005), Gerkema et al. (2006), Rodenborn et al. (2011), andGostiaux et al. (2013).Both effects were observed here.

Overall flow behavior
Figure 4 illustrates the results typical for simulations with an incidence angle θ of 25 degrees and r = 5.It shows the amplitude of the kinetic energy per unit mass in the incident wave frequency (ω f , left) and lowest harmonic frequency (2ω f , right) modes.These were computed by Fourier transforming the time series of the velocities at each point in space and keeping only the frequencies within 10 % of the frequency of interest (0.9ω f < ω < 1.1ω f for the primary mode, and 1.8ω f < ω < 2.2ω f for the harmonic).To avoid transient phenomena, the first eight wave periods were excluded.The magnitude of the kinetic energy (|u| 2 + |w| 2 )/2 for each band is presented on a logarithmic scale in Fig. 4. The path of the incident and reflected beams are clearly exhibited in red/orange in the left panel of Fig. 4. The right panel shows the generation of a harmonic beam (orange/yellow) which propagates downward from the pycnocline with angle ϕ, as given by Eq. (7c).The harmonic mode appears strongest (orange) at two distinct generation sites.One is located at the pycnocline (z ∼ −10 cm) and is presumably associated with refraction.The other is located in the stratified layer (z ∼ −15 cm) where the buoyancy frequency is constant but the incident and reflected beams overlap.Here the harmonic generation is presumably analogous to that observed in reflections from solid surfaces.Since both sites lie along the harmonic beam path, they may both contribute to its amplitude.Other simulations with different values of the pycnocline buoyancy frequency ratio r exhibit results very similar to Fig. 4. The amplitude of the harmonic beam varies only slightly when the value of r is changed.
Figure 5 shows the filtered kinetic energy per unit mass from a simulation with θ = 40 degrees and r = 4.The primary frequency band (left) is very similar to the θ = 25 • case (Fig. 4), apart from the steeper beam angle.The harmonic band (right) exhibits generation (red/orange) at approximately the same locations as the θ = 25 • case (Fig. 5), namely at z ∼ −10 cm due to refraction, and at z ∼ −15 cm due to beam overlap.However, because the harmonic frequency 2ω f exceeds the buoyancy frequency N o in the lower layer, no propagating harmonic beam is able to form.Instead,  an interfacial wave is seen propagating along the pycnocline (z ∼ −10 cm) away from the refraction generation region.The beam overlap generation site is not spatially connected to the pycnocline interfacial wave and may not contribute significantly to its formation.Results for other values of r exhibit interfacial waves with different amplitudes; Fig. 5 shows one of the strongest cases.However, all of the simulations exhibit two distinct harmonic generation sites, with only the pycnocline refraction site appearing to be connected to the interfacial wave.
To further illustrate the characteristics of the harmonic mode, Fig. 6 presents the vertical velocity at the end of simulation with θ = 25 • and r = 5.In addition, time-filtered velocity fields for the primary and harmonic modes are shown in the middle and lower panels, respectively.These were computed using the same Fourier band-pass method as Figs. 4 and 5.A harmonic beam originating from the pycnocline incidence region is clearly seen in the bottom panel.The harmonic velocity amplitude is about 10 % of the incident beam amplitude, while its horizontal wavelength is smaller.Results for other values of r with θ = 25 • are very similar to Fig. 6.
Figure 7 presents the total and filtered vertical velocities (as in Fig. 6) for the case with θ = 40 • and r = 2.5.In this case, the harmonic mode is confined to the vicinity of the pycnocline, originating in the primary beam incidence region and propagating to the right.The maximum amplitude of the harmonic velocity field is about 10 % of the primary beam.
Figure 8 compares the trapped pycnocline harmonic mode for θ = 40 • and r = 2.5, 4.0, and 5.0.The filtered vertical velocity (as in the bottom panel of Figs. 6 and 7) is shown on a domain which is restricted to focus on the pycnocline.The vertical velocity is largest within the pycnocline (z = −9.5 to −10 cm) and propagates horizontally.The observed modes lack significant vertical structure within the pycnocline, unlike the simulations of Diamessis et al. (2014) with thicker pycnoclines.The horizontal wavelength of the harmonic increases with r, consistent with the weakly nonlinear result that the preferred harmonic interfacial wave number is inversely proportional to the pycnocline density jump ρ (Eq.6) for fixed θ.This preferred wave number lies on the  dispersion curve for an interfacial wave with a frequency that is twice the incident wave frequency.The amplitude is also smaller for larger values of r (note the different color scales).The trapped harmonics decay with distance from the incidence region, presumably due to viscous dissipation.This effect appears strongest for the shortest wavelength modes (r = 2.5).

Comparison with weakly nonlinear theory
For incidence angles exceeding 30 degrees, the weakly nonlinear refraction calculation suggests that the amplitude of the trapped harmonic will be largest for the wave number which lies on the interfacial wave dispersion curve for the given harmonic frequency.To evaluate this conjecture, twodimensional spectra of kinetic energy were computed along the pycnocline downstream of the internal wave beam incidence location.The sampled depth interval was from z = −10.5 cm to z = −10.0cm (just below the pycnocline).The selected horizontal interval extended from x = 20 cm (near beam incidence) downstream to x = 60 cm for the smallest value of r and to x = 200 cm for the largest.Longer horizontal intervals were used for larger values of r to capture the longer wavelengths (as shown in Fig. 8), while shorter values were used for smaller r to minimize the effect of viscous dissipation on the shorter wavelengths.The first 10 wave periods were excluded, to allow time for the interfacial wave to form and reach a steady state, and Welch windowing was used.
Figure 9 presents the results for three different values of r.The dispersion relation of an interfacial wave is shown by dotted white lines in Fig. 9, which crosses ω = 2ω f at wave number 0.9 cm −1 for r = 2.5, 0.4 cm −1 for r = 4, and 0.2 cm −1 for r = 6.These wave numbers correspond to the value of γ resulting in a singularity in the harmonic amplitude in Eq. ( 6).(For θ = 40 • , γ ∼ 1.35.)The harmonic amplitude is also proportional to the square of the incident beam amplitude, so formation of a strong harmonic at the natural interfacial wave number requires significant incident energy at half this wave number.The incident beam kinetic energy spectrum peaks at 2k o ∼ 0.85 cm −1 , which lies closest to the natural interfacial wave for r = 2.5 (0.9 cm −1 ).As expected, the strongest harmonic amplitude was observed for r = 2.5, with weaker harmonics for larger values of r.The harmonic spectrum peaks near the white curves for r = 2.5 and r = 4, which lie within the band of significant incident energy.For r = 6, the white curve has moved far outside this band, and no peak is found at the natural interfacial wave number due to the lack of incident energy.All of these observations are consistent with the selection of the trapped harmonic wave number predicted by the weakly nonlinear theory.
Figure 10 presents a more quantitative comparison of the harmonic amplitudes observed in the numerical simulations and the weakly nonlinear refraction calculation.The black (θ = 40 • ) and blue (θ = 25 • ) lines show the weakly nonlinear prediction for the harmonic kinetic energy of a plane wave as a function of γ (k) = r 2 kh, normalized by the incident beam kinetic energy.These are computed from Eqs. (6b) and (7b) using the definition of the stream function, yielding Comparable simulation results for the normalized kinetic energy are shown by symbols in Fig. 10.The value of γ for each simulation is based on the wave number k o of the peak in the incident energy spectrum (Table 1), and is equal to r 2 k o h.The incident wave beam kinetic energy was found by applying the same algorithm to the crossbeam spectra, such as those shown in Fig. 3.The incident beam amplitude was the same for all simulations with the same incidence angle (7.0 × 10 −5 cm 2 s −2 for θ = 25 • and 1.6 × 10 −4 cm 2 s −2 for θ = 40 • ), so the dependence on γ (k o ) in Fig. 10 is due entirely to changes in the harmonic amplitude.For the θ = 40 • cases (black circles), the harmonic amplitude was computed from the spectra shown in Fig. 9 by filtering to keep frequencies within 10 % of the harmonic 2ω f and then integrating over all wave numbers to obtain a total kinetic energy amplitude for the harmonic frequency.For simulations with θ = 40 • , the normalized harmonic amplitude is in reasonable agreement with the weakly nonlinear result (Eq.13a).This is in spite of the fact that Eq. ( 13a) assumes a plane wave with a single wave number, while the numerical simulations have an incident beam with a wave number spectrum peaked at k o .The strongest harmonic interfacial wave is found for γ (k o ) = 1.3, which corresponds to the approximate location of the singularity in the weakly nonlinear refraction result (Eq.13a).The actual amplitude in the numerical simulation result is finite, as the weakly nonlinear calculation neglects fully nonlinear effects that are important at this wave number.
The harmonic amplitude decreases for other values of γ (k o ) as the harmonic of the incident wave number peak moves away from the interfacial wave dispersion curve.For θ = 25 • , the kinetic energy in the harmonic beam was computed from a cross-beam profile in the stratified lower layer, analogous to the calculation of the incident beam kinetic energy.The normalized results are shown by blue squares in Fig. 10.Except for the smallest value of γ (k o ), the harmonic amplitude is independent of the pycnocline buoyancy frequency and much larger than the weakly nonlinear prediction of Eq. (13b).This suggests that, in these cases, nonlinear harmonic generation in the beam overlap region (which is not accounted for in the plane wave calculation leading to Eq. 13b) dominates over refraction in forming the radiated harmonic beam.Hence the natural extension of the weakly nonlinear theory in Diamessis et al. (2014), which yields good agreement with simulations of trapped harmonics, fails to explain the large harmonic amplitudes found here for radiated harmonics.

Discussion and future work
The fully nonlinear simulations of an internal wave beam incident on a thin pycnocline presented here demonstrate two distinct mechanisms of harmonic generation.For beam incidence angles exceeding 30 degrees, refraction at the pycnocline generates a harmonic interfacial wave which is trapped within the pycnocline.The harmonic mode is strongest when the incident beam kinetic energy spectrum peaks at a wave number k o such that its harmonic (2ω, 2k o ) lies on the dispersion curve for a freely propagating interfacial wave.The numerical results approximately agree with a weakly nonlinear calculation for plane wave refraction at the pycnocline.For smaller incidence angles, a harmonic beam is radiated into the lower stratified layer.The overlap of the incident and reflected beams below the pycnocline appears to be the dominant mechanism for the generation of the harmonic beam.Although pycnocline refraction is also capable of generating a harmonic beam, this mechanism appeared to be less important than beam overlap for an incidence angle of 25 degrees.Hence the weakly nonlinear theory of trapped harmonics (Diamessis et al., 2014) fails to extend to smaller incidence angles with radiated harmonics.
Oceanic internal wave beams generally have incidence angles much less than 30 degrees.Nonlinear energy transfer to higher frequency radiated harmonics may play a role in the degradation of such beams as well as to the observed spectrum of internal waves in the oceans.For very small incidence angles, the normalized kinetic energy of the harmonic in Eq. ( 13b) reduces to a constant (|KE H |/|KE inc | 2 = 32k 2 /ω 2 ), independent of incidence angle or pycnocline characteristics.This is comparable to the harmonic amplitude attributed to beam overlap in the simulations.An oceanic internal wave beam might have a velocity scale of 0.1 m s −1 , as observed in the Hawaii Ocean Mixing Experiment (Cole et al., 2009), a wave number k ∼ 10 −4 m −1 , and frequency ω ∼ 2 × 10 −4 s −1 .In this example, both the weakly nonlinear analysis and the extrapolation of the simulation results yield a harmonic kinetic energy that is ∼ 10 % of the incident kinetic energy.This estimate suggests that harmonic generation might be significant in the oceans.Whether refraction or beam overlap is more important for small incidence angles is a topic of future work.
Excitation of a pycnocline interfacial wave by the lowest harmonic is only possible for incidence angles greater than 30 degrees, which is generally not the case for oceanic internal wave beams.Although the lowest harmonic will be radiated for smaller incidence angles, higher harmonics might be trapped within the pycnocline if their frequency exceeds the buoyancy frequency below the pycnocline.If a higher harmonic matched the dispersion relation for an interfacial wave, analogous to the weakly nonlinear analysis presented here for the lowest harmonic, it might attain significant amplitude.As in the simulations presented here, pycnocline refraction may be the dominant mechanism for generation of higher harmonic interfacial waves.Trapped harmonic modes have been conjectured to play a role in the local generation of internal solitary waves on the pycnocline (New andPingree, 1990, 1992;da Silva et al., 2007da Silva et al., , 2009;;Grisouard and Staquet, 2010;Grisouard et al., 2011).Estimation of the amplitudes of harmonics beyond the lowest for small incidence angles is left to future work.
Harmonics of the semi-diurnal tides have recently been observed near the pycnocline in the South China Sea (Xie et al., 2013).Tidal oscillations differ from internal wave beams, but refraction at the pycnocline would still be expected to generate harmonics according to Eq. ( 4).Future studies of semi-diurnal tides in the presence of a pycnocline are needed to evaluate the magnitude of harmonic generation by refraction in this case.Rotation may also be important here, so three-dimensional simulations would likely be needed for quantitative understanding of this effect.

Figure 1 .
Figure 1.Layout of the numerical simulations.
Figure1is a schematic of the simulations.An internal wave source with oscillation frequency ω f and thickness λ is located at the left edge of the simulation domain, in the range −4λ < z < −3λ.It generates a beam which travels upward through a uniformly stratified layer of fluid with buoyancy frequency N o until it passes through the pycnocline and reflects from the base of the mixed layer.The simulations used two different frequencies to give incidence angles θ of 25 and 40 • .The domain depth is 7.5λ, while its length is 21λ for the θ = 40 • simulations and 25λ for the θ = 25 • simulations.For the 40 • incidence angle simulations, 79 × 54 elements (6 points per element) are used, while the 25 • simulations use 93 × 54 elements.The vertical resolution varies with depth; to resolve the large density gradient in the pycnocline, there are eight elements within the pycnocline and two more on each pycnocline boundary.The horizontal resolution is uniform except in the vicinity of the beam source, where the resolution is more refined.Simulations were run for 20 internal wave beam forcing periods to reach a steady state.The internal wave beam is generated by imposing the oscillating boundary condition

Figure 2 Figure 2 .
Figure 2. Example image of the horizontal velocity u(x, z) from the simulation with θ = 25 • and r = 5.

Figure 3 .
Figure 3.An example incident internal wave beam 2-D kinetic energy (log scale) spectrum (left), which is strongly peaked at the forcing frequency, and the spectrum filtered for primary and harmonic frequencies (right).Results from the simulation with θ = 40 • and r = 4.

Figure 4 .
Figure 4. Kinetic energy amplitude in the incident frequency (ω f , left) and lowest harmonic (2ω f , right) bands for the simulation with θ = 25 degrees and r = 5, illustrating the radiated harmonic beam.

Figure 5 .
Figure 5. Kinetic energy amplitude in the incident frequency (ω f , left) and lowest harmonic (2ω f , right) bands for the simulation with θ = 40 degrees and r = 4, illustrating the trapped harmonic mode.

Figure 6 .
Figure 6.Snapshot of the vertical velocity field (top) for the simulation with θ = 25 degrees and r = 5.Time-filtered velocity fields for the ω f (middle) and 2ω f (bottom) bands show the radiated harmonic wave beam.

Figure 7 .
Figure 7. Snapshot of the vertical velocity field (top) for the simulation with θ = 40 degrees and r = 2.5.Time-filtered velocity fields for the ω f (middle) and 2ω f (bottom) bands show the trapped harmonic mode.

Figure 8 .
Figure 8.The trapped harmonic mode for r = 2.5 (top), r = 4 (middle), and r = 5 (bottom).As r increases, the wavelength of the preferred pycnocline mode increases and the amplitude decreases.

Figure 9 .
Figure 9. Kinetic energy spectra in the pycnocline for simulations with θ = 40 degrees and three values of r.The white line shows the dispersion curve for interfacial waves in each case.

Figure 10 .
Figure 10.Comparison of the kinetic energy amplitudes in the harmonic mode between the simulations (symbols) and Eq.(13) (dotted lines) from the weakly nonlinear refraction calculation.