Internal wave energy flux from density perturbations in nonlinear stratifications

Internal gravity wave energy contributes significantly to the energy budget of the oceans, affecting mixing and the thermohaline circulation. Hence it is important to determine the internal wave energy flux $\mathbf{J} = p \, \mathbf{v}$, where $p$ is the pressure perturbation field and $\mathbf{v}$ is the velocity perturbation field. However, the pressure perturbation field is not directly accessible in laboratory or field observations. Previously, a Green's function based method was developed to calculate the instantaneous energy flux field from a measured density perturbation field $\rho(x,z,t)$, given a $constant$ buoyancy frequency $N$. Here we present methods for computing the instantaneous energy flux $\mathbf{J}(x,z,t)$ for a spatially varying $N(z)$, as in the oceans where $N(z)$ typically decreases by two orders of magnitude from shallow water to the deep ocean. Analytic methods are presented for computing $\mathbf{J}(x,z,t)$ from a density perturbation field for $N(z)$ varying linearly with $z$ and for $N^2(z)$ varying as $\tanh (z)$. To generalize this approach to $arbitrary$ $N(z)$, we present a computational method for obtaining $\mathbf{J}(x,z,t)$. The results for $\mathbf{J}(x,z,t)$ for the different cases agree well with results from direct numerical simulations of the Navier-Stokes equations. Our computational method can be applied to any density perturbation data using the MATLAB graphical user interface EnergyFlux.


Introduction
Ubiquitous internal gravity waves are generated in the oceans by tidal flow over bottom topography and by surface storms (Munk & Wunsch 1998;Alford 2003;Wunsch & Ferrari 2004). The internal waves transmit energy from generation sites to large distances, and ultimately the energy is converted into small-scale mixing. Internal waves are a major contributor of the energy budget of the oceans, and the mechanism for this contribution can be better understood through the energy flux field. In this paper, we examine the baroclinic energy flux J in the presence of a stable background density stratification, where p is the pressure perturbation from the background hydrostatic pressure field, and v is the velocity perturbation from the background velocity flow field. Determining the energy flux requires the simultaneous measurement of both the pressure and velocity perturbation fields. In numerical investigations of internal waves, these fields are computed (Lamb 2004;Niwa & Hibiya 2004;Zilberman et al. 2009;King et al. 2009King et al. , 2010Gayen & Sarkar 2010, 2011Rapaka et al. 2013), enabling a direct calculation of the energy flux. However, In laboratory and field studies simultaneous measurements of the velocity and pressure perturbation fields are usually not possible.
Laboratory observations of internal waves have been made by particle image velocimetry (PIV) (Echeverri et al. 2009;King et al. 2009King et al. , 2010Paoletti et al. 2014), synthetic schlieren (Sutherland et al. 1999;Dalziel et al. 2000;Clark & Sutherland 2010;Allshouse et al. 2016), and light attenuation (Dossmann et al. 2016). In PIV, neutrally buoyant particles scatter incident laser light, and movies of the scattered light field yield the time-varying velocity field. In the schlieren method, a patterned mask is placed behind a tank that contains the internal waves, and the time-varying distorted image formed by light transmitted through the tank in the transverse direction gives the timedependent density perturbation field (Sutherland et al. 1999;Dalziel et al. 2000). The light attenuation technique uses dye intensity that varies with depth and moves with the fluid to measure the density field as a function of time (Dossmann et al. 2016). These experimental techniques yield the velocity field in the case of PIV and the density perturbation field in the case of synthetic schlieren and light attenuation measurements. No experimental technique directly yields the pressure perturbation field needed for the energy flux calculation.
Given the importance of the energy flux and the far-field radiated power obtained by integrating the flux, multiple efforts have been made to calculate the energy flux from experimental measurements. One method assumed a constant buoyancy frequency and calculated the energy flux (averaged over a tidal period) given only the stream function, thus eliminating the need to measure the pressure perturbation field (Balmforth et al. 2002). The stream function method was subsequently extended to a buoyancy frequency varying exponentially with depth Lee et al. (2014). Another method applied the polarization relations to density perturbation data to obtain estimates for the velocity and pressure perturbation amplitudes (Clark & Sutherland 2010). This method provided the energy flux time averaged over a tidal period for a monochromatic plane wave, which is not representative of the complex ocean internal wave field, which has many natural frequencies and spatially varying wave packets. Finally, an ocean observation analysis technique calculated the pressure perturbation field for measured density profiles, assuming a hydrostatic pressure field; this together with simultaneous velocity measurements allowed the calculation of the energy flux field (Nash et al. 2005). This method assumed that the contribution of dynamic pressure is negligible, which is often not the case in experiments and in high velocity events in the ocean.
The method presented here for obtaining the instantaneous energy flux from laboratory and field data differs from previous methods that determined the time-averaged energy flux or a global energy conversion rate. Recently a Green's function method was used to calculate the instantaneous energy flux field from a density perturbation field (Allshouse et al. 2016) for a fluid with constant buoyancy frequency, . Two buoyancy frequency profiles from a World Ocean Circulation Experiment data set: (a) A buoyancy frequency squared profile (black) fit to a tanh profile (red). (b) A buoyancy frequency profile (black) fit to a linear profile (red). The insets show regions 1000 km × 1000 km that contain the locations (red dots) where the measurements were made. The mean buoyancy frequency (squared for (a)) of bins of stratification measurements is plotted as a function of depth (black dots) with error bars representing two standard deviations from the mean. but in the oceans N varies significantly with depth, as figure 1 illustrates with data from two locations.
Recognizing the strength of the Green's function method, we expand that method to accommodate linear N and tanh N 2 profiles. These two profiles were selected due to their mathematical properties and their presence in ocean stratifications. The tanh N 2 profile is often a good approximation in the ocean when two near-constant buoyancy-frequency zones are separated by a pycnocline, as figure 1(a) illustrates. The linear N profile can occur in the ocean when there is no strong pycnocline, and the buoyancy frequency near the surface decays to a near constant value in the depth, as illustrated in figure 1 However, as many buoyancy frequency profiles cannot be adequately approximated by either a linear N or a tanh N 2 profile as was done in figure 1, the general N (z) case must be treated separately. To account for this, we present a numerical method for computing the instantaneous energy flux field and the integrated far-field internal wave power solely from a density perturbation field, which can have an arbitrary N (z) profile.
The present work provides a tool for laboratory experiments and field measurements: the calculation of the instantaneous energy flux field from density perturbation data. The method can be applied to ocean density perturbation space-time data when such data becomes available. The methods presented here and in Allshouse et al. (2016) provide the instantaneous rather than time-averaged energy flux field. Thus the resultant energy flux and integrated far-field power include all spectral components, while previous methods provided only the global conversion rates or monochromatic results.
The paper is organized as follows. An outline of our method for obtaining the energy flux from the density perturbation field in a tanh N 2 and linear N stratification is presented in §2. This method is then verified with numerical simulations in §3. A finite difference method for calculating the energy flux for an arbitrary buoyancy frequency profile is presented in §4, and is applied to an ocean-inspired stratification. Lastly, conclusions and potential applications of this work are presented in §5.

Theoretical development
As this work builds off the theoretical foundation presented in Allshouse et al. (2016), we present the general equations in §2.1. These equations provide the velocity perturbation field as a function of the density perturbation, and a functional relationship between the density and pressure perturbation fields is established. For an analytic tanh N 2 and linear N , we calculate the Green's function in §2.2 and §2.3, respectively.

Generalities
Our goal of obtaining the energy flux (1.1) from the density perturbation field alone requires calculating the pressure perturbation and components of the velocity perturbation from the density perturbation field. The details of these calculations were given in Allshouse et al. (2016) for a uniform N profile, but here we present a condensed version of the pressure perturbation calculation, as needed for the calculations for the tanh N 2 and linear N profiles.
We begin with the linearized Euler equations for perturbation about a hydrostatic background, where u and v are the horizontal and vertical components of the velocity perturbation, respectively, p is the pressure perturbation, ρ is the density perturbation, ρ 0 is the hydrostatic background density profile, and N is the buoyancy frequency. By manipulating (2.1) and (2.2) we obtain an equation for the pressure perturbation in terms of the density perturbation, (2.3) First, we solve this equation for p assuming we have the measured ρ. Equation (2.3) is brought into a convenient form by first applying the following transformation: and then Fourier expanding in x, yielding Here F (z; k) and Q(z; k) denote the Fourier components of and q(x, z), respectively. We solve (2.6) for Q given F by obtaining the Green's function for the Fourier components, which satisfies with a no-flux condition in the z direction at the top and bottom of the domain, (2.9) and the Green's function matching conditions, Thus, given profiles for the buoyancy frequency N and source term f , the pressure perturbation is given by the following expression: where k = 2πn/l, l is the width of the system, and n is a positive integer. Next, we obtain the components of the velocity perturbation. The vertical component follows by inverting the first equation of (2.2) yielding, and the horizontal component is obtained by using the vertical velocity perturbation (2.13) and the incompressibility condition, the second equation of (2.2), which gives the differential equation (2.14) None of these calculations depend on the particular form of the buoyancy frequency profile, so it is possible to perform all the necessary expressions for calculating the energy flux from ρ alone in a general stratification. To calculate analytically the Green's function for (2.8), it is necessary that the functional form of the buoyancy frequency profile be specified. Allshouse et al. (2016) investigates the particular case where the buoyancy frequency is constant resulting in a Green's function that is exponential. In §2.2 we present the calculations for obtaining the pressure perturbation for the tanh N 2 profile, and in §2.3 we present the analogous calculation for the linear N profile.

The tanh profile
The buoyancy frequency profile we assume in this section is given by because this gives a convenient form for N dN/dz. Here α controls the transition width between the two buoyancy frequency values N 1 and N 2 , and z t is the midpoint of the transition. Note, for large α (2.16) approximates a two-layer N 2 profile, which we will investigate in §3.2. We assume that N 4 /4g 2 in (2.16) is negligible. For low mode numbers, outside of the transition region, k 2 (∼ 10 0 − 10 1 in MKS) is much larger than N 4 /4g 2 (∼ 10 −2 ), and near the transition region k 2 is roughly the same order as (N/g)(dN/dz). For larger mode numbers k 2 is the dominating term. Thus for simplicity we keep k 2 and (N/g)(dN/dz) and drop N 4 /4g 2 for all modes. Upon substituting the analytic stratification into the Green's function equation, (2.8) becomes Equation (2.17) is of the form of a well-studied (e.g. Epstein (1930); Pöshl & Teller (1933); Lekner (2007)) time-independent Schrödinger equation.
With the dimensionless coordinate transformation where the dimensionless Green's functionḠ is given by and the parameters ν and µ are given by Thus the transformation takes (2.17) into the associated Legendre equation (2.19), which has the two linearly independent solutions P µ ν (y) and Q µ ν (y), the associated Legendre functions of the first and second kind, respectively. Then, solving (2.19) with the boundary conditions (2.9) and the matching conditions (2.10) and (2.11) gives the Green's function, ( 2.22) Here Note, the transformation factor T (z) for this case is given by Internal wave energy flux from density perturbations in nonlinear stratifications 7 2.3. The linear profile The calculations for the linear N profile are similar to those of §2.2, so we only highlight the important differences. The linear profile for the buoyancy frequency is given by where z t is now the location where the buoyancy frequency becomes zero. We again neglect N 4 /4g 2 (∼ 10 −2 ) in comparison to k 2 and (N/g)(dN/dz) (∼ n 2 and ∼ 1, respectively) as we insert (2.27) into (2.8). For the linear N profile, instead of (2.17) we obtain (2.28) With the coordinate transformation where once again y is a dimensionless coordinate variable, equation (2.28) becomes which is the Airy equation with the two independent solutions Ai(y) and Bi(y), the Airy functions of the first and second kind, respectively. Then, the dimensionless Green's function is given bȳ which when given dimensions becomes G(z(y)) = g 1/3 (N ) −2/3Ḡ (y). (2.32) Here where z 0 and z h are the coordinates of the bottom and top of the domain, respectively, and y 0 and y h are the corresponding transformed coordinates. The transformation factor T (z) in this case is given by (2.36)

Analysis verification
To verify the Green's function analysis in §2, we compare those predictions with results for the energy flux obtained from direct numerical simulations of the Navier-Stokes equations. The simulations are described in §3.1. The simulated velocity perturbation, pressure perturbation, and energy flux fields of internal waves in a stratified fluid are compared with the predictions from the analyses for a tanh N 2 profile in §3.2 and for a linear N profile in §3.3.

Simulation of the density perturbation field
To verify the Green's function method, we perform direct numerical simulations of the Navier-Stokes equations in the Boussinesq approximation. These simulations provide the density perturbation field needed to calculate the velocity perturbation, pressure perturbation, and energy flux fields. The simulations use the CDP-2.4 algorithm, which is a finite volume solver that implements a fractional-step time-marching scheme (Ham & Iaccarino 2004;Mahesh et al. 2004). This code has previously been used to simulate internal waves and has been validated with experiments (King et al. 2009;Lee et al. 2014;Dettner et al. 2013;Zhang & Swinney 2014;Paoletti et al. 2014;Allshouse et al. 2016).
Our two-dimensional simulations span the domain x ∈ [−1.5, 3] m and z ∈ [0, 1.5] m. The simulation solves for the total density ρ T , pressure p T , and velocity u T : where ρ 00 = 1000 kg/m 3 (density of water), ν w = 10 −6 m 2 /s (kinematic viscosity of water at 20 o C), and κ s = 2 × 10 −9 m 2 /s (the diffusivity of NaCl in water). The system is initially at rest and the prescribed density field is unperturbed. The initial density field is analytically derived from the buoyancy frequency profiles presented in figure 2(a). The boundary conditions at the bottom and top are no slip and free slip, respectively. The left and right boundaries are set to be periodic; however, Rayleigh damping is used along the perimeter of the domain (gray region in figure 2(b)), thus forcing the velocity to be negligible at the left and right boundary. The internal wave beam is produced by using a momentum source that forms a rectangle with height 0.15 m and width 0.04 m, centered at (−0.02, 0.8) m and rotated to match the internal wave beam angle corresponding to the buoyancy frequency at z = 0.8 m. The wave beam velocity imposed is with an amplitude profile given by where the lengths are in meters, the rotated coordinates x and z correspond to the beam tangent and normal coordinates centered at x = (−0.02, 0.8) m, respectively, ω = 2π/13 rad/sec and k z = 8245 m −1 . A time step size δt = 0.0025 s (5200 steps per period) is sufficient for temporal convergence. Spatial convergence is obtained using an unstructured mesh with resolution δx ≈ 0.0014 m inside the region x ∈ [−0.8, 1.80] m, y ∈ [0.5, 1.1] m. This high resolution region contains the beam generation, the density gradient transition for the tanh N 2 profiles, and generation of any additional beams. The resolution outside of this region grows to δx ≈ 0.0025 m near the boundaries. Changes in the velocity field are less than 1% when spatial and temporal resolutions are doubled. The density perturbation field for the case where we have a rapid change in buoyancy frequency (blue line in figure 2(a)) is presented in figure 2(b). The internal wave beam is generated at (−0.02, 0.8) m and produces a beam propagating to the right that is the focus of our studies and a beam propagating to the left which is damped out by the Rayleigh damping. The beam propagating down to the right reaches the interface at z = 0.6 m at which point three beams are produced: a reflected beam to the top right at the same angle to the horizontal as the incoming beam, a transmitted beam to the bottom right that has a different angle, and a reflected second harmonic beam at approximately twice the incoming angle. This particular snapshot is shown after 23.06 periods of forcing, which is sufficient for the beam in the region of interest to reach steady state.

Tanh N 2 profile analysis verification
The vertical (2.2)a and horizontal (2.14) components of the velocity and the pressure perturbation calculated from the density perturbation using the Green's function (2.22) for the tanh N 2 profile are verified by comparison with the direct numerical simulations described in §3.1. For large α the tanh N 2 profile can be approximated as a two-layer N system, as illustrated in figure 2(a), where α = 4, corresponding to a transition thickness of 0.01 m for a 95% change in N 2 ; this is at least an order of magnitude smaller than the the beam width and domain height. Henceforth the large α case is called the "narrowtransition" tanh N 2 profile.
The difference between the density-perturbation-based velocity perturbation and the simulated velocity perturbation for the narrow-transition tanh N 2 profile are presented in figure 3; in most of the domain, the difference is less than 3% (with respect to the beam amplitude), except in the transition region at z = 0.6 m where there is a significant amount of nonlinearity. Because the horizontal velocity perturbation is found by solving an ODE on constant z levels (2.14), the patch of large vertical velocity perturbation error is propagated horizontally from the reflection site; thus the region of error is larger for u. Despite this nonlinearity, the error is small in most of the domain.
Next we investigate how well the Green's function method calculates the pressure perturbation field from the density perturbation field. Figure 4(a) shows the simulated pressure perturbation field, and figure 4(b) shows the pressure perturbation calculated using the Green's function method. Despite the nonlinearities in the narrow transition layer, the Green's function method, which is based on the linear equations, yields accurate estimates of the pressure perturbation for the reflected, transmitted, and second harmonic beams, as figure 4(b) illustrates. The beam-normalized percent difference between the calculated and simulated pressure perturbation fields is presented in figure 4(c). The calculation is accurate to within 5% over most of the domain, and to better than 10% everywhere except within 0.02 m of where the beam enters the domain. Near the top of the domain, the Green's function method overestimates the pressure perturbation by 4 − 6%, which causes some distortion in the second harmonic, as can be seen around (1.25,0.8) m in figures 4(a) and 4(b). The Green's function method underestimates the pressure perturbation in the center of the domain, but the error is less than 5%.
Finally, we use the calculated velocity and pressure perturbation fields to compare the energy flux J directly from the numerical simulations with the flux computed from the Green's function analysis. The magnitude of the energy flux from the simulations is presented for the narrow tanh N 2 transition regions in figures 5(a). For the narrow transition region case the energy flux for the reflected beam is higher than for the transmitted beam and an order of magnitude greater than in the second harmonic. The beam-normalized percent difference of the horizontal and vertical energy flux are presented in figures 5(b) and (c), respectively. Outside of the immediate vicinity of the interface region at z = 0.65 m the percent difference is less than 3%. The accumulated error from multiplying the calculated velocity and pressure perturbation to obtain the flux components is large at the narrow transition interface as a consequence of error in the horizontal velocity perturbation, which is compensated to some extent by a more accurate pressure perturbation calculation at the interface; the error is smaller for J z (x, z) than for J x (x, z). In the lower half of the domain the magnitude of the energy flux is underestimated due to underestimation of the pressure perturbation. The insets show that the error along three beam transects is mostly smaller than 3% for the narrow transition simulation.
We also simulate a tanh N 2 profile with a broader transition thickness layer of 0.31 m. We omit the comparison of the velocity and pressure perturbations for brevity and instead examine the energy fluxes, as shown in figure 5. The energy flux field for the broader tanh N 2 profile is presented in figure 5(d). The internal wave beam passes through the broad transition without reflection because there are no rapid changes in buoyancy frequency. This smooth transition reduces the nonlinearities so there are significantly smaller errors in the velocity perturbation field and thus the energy flux field as compared to the narrow transition tanh N 2 . The magnitude of the energy flux decreases as the beam widens in the bottom of the domain and then increases again as the beam narrows after reflection. Beam-normalized percent differences are presented for the horizontal and vertical energy flux in figures 5(e) and (f), respectively. There is a change of overestimating the energy flux in the top of the domain to underestimating the energy flux in the bottom of the domain. This is most clearly seen at (0.7,0.7) m where the bands of constant phase change from red to blue and vice versa. This change is due to the pressure perturbation again being over estimated in the top of the domain and underestimated in the bottom of the domain. The two insets show that within the beam the percent difference is consistently less than 5%. Figures 5(a)-(c) demonstrate that our method can handle rapid changes N 2 , while the broad N 2 transition thickness in figures 5(d)-(f) is more representative of ocean stratifications. Figure 5(d) shows the energy flux amplitudes and reveals that the broadening of the transition layer eliminates the reflected and second harmonic beams. Further, the error in the broad transition region in figures 5(e)-(f) is much smaller than in the Figure 6. (a) The energy flux magnitude from the numerical simulation for the linear N profile. The beam-normalized percent difference between simulation and the Green's function method for the (b) x and (c) z components of the energy flux; the difference is less than 3% for most of the beam, as illustrated by insets showing the difference for two beam transects. narrow transition region figures 5(b)-(c). The errors are less than 5% except in the narrow transition region.

Linear profile analysis verification
To verify that the theory for the linear N profile of §2.3 is valid, we perform simulations analogous to those in §3.2. The energy flux field in figure 6(a) demonstrates that the internal wave beam bends more gradually for the linear N profile (figure 6(a)) as compared to the tanh N 2 profiles discussed in §3.2 ( figure 5(a) and (d)). This slower change is due to the smaller gradient of the buoyancy frequency for the linear N profile. Again, because there are no rapid changes in N there are no reflection depths other than the bottom of the system, so nonlinearities are limited to the reflection point at (1.6, 0) m.
We present only the errors in the energy flux calculation; the errors in the velocity and pressure perturbation calculations are qualitatively the same as the results in figures 3 and 4. Figures 6(b) and (c) show the percent difference of the horizontal and vertical components of the energy flux, respectively. As with the tanh N 2 profile comparisons, the errors in the energy flux field are confined to the internal wave beam. Throughout the beam the difference between the simulation and Green's function method is less than 5%, as illustrated by the beam transects in the insets of figures 6(b) and (c); the largest errors occur where the beam enters and leaves the domain and where it reflects off the bottom boundary. The transition from pressure perturbation overestimation to underestimation is highlighted by the change from red to blue and vice versa near (0.5, 0.75) m.

Arbitrary stratification analysis
Implementation of the Green's function method works for systems with particular stratifications only when an analytic representation of the Green's function exists. While some stratifications in the ocean and laboratory may approximately fit to these particular stratification profiles as we show in figure 1, making this density-perturbation-based calculation more general is necessary for most applications. To accomplish this generalization, we use a finite difference method to determine the pressure perturbation field. We present the method in detail along with a comparison between the Green's function method and the finite difference method in §4.1. Then, we apply the finite difference method to an ocean-inspired stratification in §4.2.

Finite difference method
Since the velocity perturbation calculation does not depend on having an analytic stratification, only the calculation of the pressure perturbation field requires modification for application to general stratifications. This is accomplished by implementing a numerical solver of the second order differential equation (2.6). The boundary conditions for this differential equation are analogous to (2.9): We solve equation (2.6) using a second-order finite difference method. The Robin boundary conditions are calculated to second order by adding ghost points to the top and bottom of the domain. This numerical method is applied to both the real and imaginary components for every Fourier mode. After the calculation of Q(z; k) using the finite difference method, the dependence in the x-direction is accounted for by multiplying by the particular Fourier mode just as it is done for the Green's function method. Finally, the transformation (2.5) is performed to determine the contribution to the pressure perturbation field by that particular mode. Applying this strategy to the previous analytic stratifications provides a baseline for comparison to the Green's function method. The percent difference of the pressure perturbation fields relative to the Green's function results are presented in figure 7. This figure shows that the pressure perturbation fields calculated using the finite difference method are everywhere less than 5% different for the tanh profile and less than 1% different for the linear buoyancy frequency profile when compared to the Green's function pressure perturbation. The only major discrepancy between the two methods is near the narrow transition in the tanh profile shown in figure 7(a). In this region, the Green's function method is consistently more accurate than the finite difference method. This is highlighted in figure 7(c) by comparing pressure perturbation profiles just above the transition layer. The discrepancy is likely due to the Green's function's accurately accounting for the rapid change in the buoyancy frequency when it modifies the coefficients in the calculation of the Legendre functions. The length scale of the transformed coordinate variable y (2.18), is set by the steepness coefficient α. This increases the spacial resolution at rapid transitions. The finite difference method can only account for variations on the scale of the original data set step size, which, in the case of the narrow transition, is too coarse.

Verification of the finite difference method
To further validate the finite difference method, we apply the method to a stratification that does not have a simple analytic function as has been the case in §2 and §3. The stratification we simulate is based on a density profile measured in the ocean during the World Ocean Circulation Experiment (WOCE). The particular profile presented in figure 8(a) was measured at 165 • W, 51.5 • N on September 20 th , 1994. This profile features two layers of large density gradient similar to the transitions of the tanh N 2 profiles. The first, more abrupt layer is centered at 30 m below the surface and the second layer is centered at 100 m. The full profile extends to a depth of 1000 m, but there is little variation in the buoyancy frequency below 200 m.
In order to simulate the beams in a similar domain and time scale as the analytic stratifications, we rescale the vertical coordinate and the density. We note that this is done to mimic the actual ocean profile and use it as an inspiration, rather than to model it accurately. Because the length scale of the transition layers in the simulation are small compared to the length scale of the beams the rescaled simulation done here stress-tests the method.
The first adjustment we make is to provide additional vertical space above the stratification features so that the internal wave beam is fully developed and the resulting reflection off the top of the first transition is visible. The vertical coordinate is scaled from 200 m to 0.8 m in the simulation. The density is also modified to increase the buoyancy frequency, so that the values of the buoyancy frequency are comparable to those used in the Green's function verification. The minimum buoyancy frequency of the scaled density profile is N = 0.55 rad/s and the maximum value is N = 2.40 rad/s. Finally, we shift the location of forcing to be at (0.2,1.2) m to have the internal wave beam enter from the top to demonstrate the flexibility of the domain of measurement. The time scale and forcing periodicity match the previous simulations.
The magnitude of the energy flux field is presented in figure 8(b). There are a number of reflections and transmissions due to the more complicated density profile. For the first transition layer, the internal wave beam produces reflections off the top and bottom of the pycnocline layer, which can be seen at (0.5, 0.8) m and (1.0,0.8) m, respectively. In addition to the reflected energy, some of the internal wave energy is trapped in the pycnocline layer and is transported to the right (e.g. (1.25, 0.7) m). A large fraction of the energy however is transmitted through the layer. Very little energy is reflected off the second layer, allowing the rest of the energy to reflect off the bottom of the domain.
The finite difference method is applied to the modified ocean density profile, and the beam-normalized percent difference of the energy flux magnitude is presented in figure 8(c). The largest errors occur near the more abrupt transition layer. The maximum percent difference in this region 28.1%. There is no consistent trend with regards to under or over estimating the energy flux. Outside of the immediate region of the sharper transition, the percent difference is generally within 5%. It is also important to note that the method is able to capture and accurately determine the energy flux in the reflected, transmitted, and trapped internal waves outside of the highly nonlinear first reflection region.

Conclusions
We have presented two methods for calculating for internal waves the instantaneous energy flux field using only density perturbation field data. Both methods are applicable to nonlinear stratifications: the first method, a Green's function method, uses convenient analytic density stratification profiles, while the second, a finite difference method, applies to arbitrary stratification profiles.
Using our Green's function method we obtained the instantaneous energy flux from density perturbations for two buoyancy frequency profiles: one linear in z and the other where N (z) 2 ∝ tanh(z). The difference between the Green's function method and our direct numerical simulations is less than 3% outside of regions containing significant nonlinearity. Despite the Green's function method being based on linear theory, it accurately predicts the energy flux in the transmitted, reflected, and second harmonic beams, which involve significant nonlinearities.
With our finite difference method we showed how to capture the energy flux in an internal wave field containing nonlinear interactions, wave beam reflections, and second harmonic beams for any buoyancy frequency profile N (z). This method was compared with the Green's function method and direct numerical simulations, and again the errors are less than 3% for most of the domain.
The two methods presented here and in Allshouse et al. (2016) allow detailed studies of the entire instantaneous energy flux field for internal wave field data, as contrasted with methods that yield a single global conversion rate or a time-averaged result. Our methods can be used to determine the instantaneous velocity perturbation, pressure perturbation, and energy flux fields from density perturbation data obtained in experiments using synthetic schlieren or light attenuation measurements. We emphasize that the methods require only the density perturbation field over time and the background buoyancy frequency profile. Application to ocean observations will be possible provided a timevarying density perturbation field can be measured. The methods assume the flow is two dimensional, but future work could extend the method to weakly three-dimensional flows as in ocean applications.
The Matlab GUI "EnergyF lux" developed in Allshouse et al. (2016) is extended to include the methods discussed in this paper. The GUI requires density perturbation data, domain coordinates, time step size, and the N (z) profile. A manual and tutorial that reproduces the results in this work is provided to make possible straightforward applications of the methods presented here. and PJM were supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, under Award Number DE-FG02-04ER-54742.