Spectral methods for internal waves : indistinguishable density profiles and double-humped solitary waves

Internal solitary waves are widely observed in both the oceans and large lakes. They can be described by a variety of mathematical theories, covering the full spectrum from first order asymptotic theory (i.e. Korteweg-de Vries, or KdV, theory), through higher order extensions of weakly nonlinear-weakly nonhydrostatic theory, to fully nonlinear-weakly nonhydrostatic theories and finally exact theory based on the Dubreil-Jacotin-Long (DJL) equation that is formally equivalent to the full set of Euler equations. We discuss how spectral and pseudospectral methods allow for the computation of novel phenomena in both approximate and exact theories. In particular we construct markedly different density profiles for which the coefficients in the KdV theory are very nearly identical. These two density profiles yield qualitatively different behaviour for both exact, or fully nonlinear, waves computed using the DJL equation and in dynamic simulations of the time dependent Euler equations. For exact, DJL, theory we compute exact solitary waves with two-scales, or so-called double-humped waves.


Introduction
Internal solitary-like waves (henceforth ISWs) are a commonly observed feature of many natural bodies of water including both the deep and coastal oceans as well as large lakes (Bogucki et al., 1997;Hosegood and van Haren, 2004;Bogucki et al., 2005;Carter et al., 2005;Moum and Smyth, 2006;Moum et al., 2007).Being both nonhydrostatic and finite amplitude, naturally occurring ISWs cannot be accurately represented by hydrostatic numerical models (including all models presently used for global ocean dynamics and climate simulations) or linear wave theories.ISWs can Correspondence to: M. Dunphy (mdunphy@math.uwaterloo.ca)be described mathematically by both weakly nonlinear (see the reviews by Grimshaw, 1997 andHelfrich andMelville, 2006) and exact theories (Turkington et al., 1991).While observed waves often lie outside its formal domain of applicability, weakly nonlinear theory is commonly used to interpret field measurements (Moum and Smyth, 2006;Trevorrow, 1998).ISWs induce currents throughout the water column and hence have implications for physical processes such as small scale mixing and sediment resuspension, as well as chemical and biological transport.Hence, the understanding of their roles in physical processes, and the knowledge of appropriate mathematical models for a given process, are important first steps in constructing suitable parametrizations of the effects of these waves in larger scale models.
The various eigenvalue and boundary value problems in the theoretical description of ISWs, as well as the integral expressions used to obtain the coefficients in model wave equations describing their horizontal propagation, can be computed by a variety of techniques.In this article we describe the spectral and pseudo-spectral approaches, which yield highly accurate results at moderate grid resolutions.Moreover this approach can be used to rapidly search parameter space and thus test various hypotheses.We illustrate how the accuracy of the spectral approach can yield results that are essentially impossible to compute using low-order (e.g.finite difference methods) for both approximate and exact theory.Spectral methods have been used to solve the nonlinear wave equations of weakly nonlinear theory for ISWs in the past (e.g.Grimshaw and Smyth, 1986), building on the seminal work of Fornberg and Whitham (1978).
For leading order (i.e.KdV) weakly nonlinear theory we demonstrate that two markedly different density profiles yield coefficients that are essentially identical, as far as the KdV theory is concerned.We demonstrate that in these cases, exact solitary waves have a qualitatively different upper bound for the two stratifications, with only one of the two yielding wave overturning with the possible formation of trapped cores.We subsequently confirm that the fully nonlinear theory result is consistent with time dependent simulations of wave formation over isolated topography.In the Appendix we discuss the relation of the fully nonlinear theory to the weakly nonlinear theory with quadratic nonlinearity that leads to the Gardner equation (reviewed in Helfrich andMelville, 2006 andGrimshaw et al., 2007).
For fully nonlinear theory, we extend the results of Lamb and Wan (1998) to explicitly construct exact solitary waves with multiple length scales, or so-called double-humped solitary waves.Double-humped solitary waves have been considered from an asymptotic point of view in the recent article by Makarenko et al. (2009) as well as in the Russian language literature (Borisov and Derzho, 1990 and several papers discussed in Makarenko et al., 2009).In the context of multi-layer (as opposed to continuously stratified) fluids double-humped waves are similar to limiting solitary waves with overhangs (see Rusas and Grue, 2002 and the references therein), though we are not aware of an analogue for continuously stratified fluids (possibly because most DJL theory makes the Boussinesq approximation and hence restricts the magnitude of the density change).For both case studies we discuss why spectral methods are vital to carry out the presented computations.

Descriptions of ISWs
We consider a nonrotating, incompressible, inviscid fluid under the Boussinesq and rigid lid approximations in a fixed frame of reference.The origin is found at the ocean surface, the x-axis runs parallel to the flat ocean surface and the zaxis points upward ( k is the upward pointing unit vector).The governing equations read, where we have divided the momentum Eq. ( 1) by the constant reference density ρ 0 and absorbed the constant into the pressure, P , as is conventional.Throughout, we assume the ocean bottom is flat, found at z = −H , and that the fluid motion is two-dimensional with no background current (this is not necessary, but it does simplify some of the algebra in what follows).The incompressibility of the fluid implies that there exists a streamfunction, ψ so that (u,w) = (ψ z ,−ψ x ) where subscripts denote partial derivatives.
If one looks for traveling waves of permanent form and switches into a frame moving with the unknown wave speed, the Euler equations can be reduced to a single equation (a nonlinear, elliptic, eigenvalue problem) for the isopycnal displacement, η(x,z).Denoting the far upstream, or background density as ρ(z) we write the density equation as ρ = ρ(z − η).Some algebra (Turkington et al., 1991) then yields the Dubreil-Jacotin-Long equation where is the definition of the buoyancy frequency squared (the reference density, ρ 0 , has been scaled out) and the propagation speed, c, is to be determined as part of the solution.Once η and c are known, the wave-induced velocities can be computed from the relation ψ = cη.Unless the background density profile is linear (i.e.N 2 (z) is constant), to the best of our knowledge, there are no explicit solutions of the DJL equation.
Weakly nonlinear theory (henceforth WNL), following the original derivation of Benney (1966), albeit in the notation of Lamb and Yan (1996) (also see this latter article for details of the derivation), considers an expansion in two small parameters, for amplitude and µ for the square of the aspect ratio.The first order, mode-1 streamfunction is given by where and The latter is the well known KdV equation with solitary wave solutions given by It can be noted that the solitary wave properties (propagation speed, c, and half-width, λ) depend not only on the wave amplitude b 0 and the linear longwave speed c lw , but also the so-called nonlinearity (r 10 ) and dispersion (r 01 ) parameters.These are given implicitly from the background density profile through the relations Note that KdV theory does not provide an error bound.Furthermore, it yields an amplitude at which waves break for any stratification, including those for which solutions of the DJL equation broaden and reach the limit of a flat centered solitary wave well below breaking.

Numerical methods
The linear eigenvalue problem used for WNL theory is solved using a Chebyshev discretization and the standard, direct eigenvalue solver in MATLAB (Weideman et al., 2000).
The DJL equation is solved iteratively on the domain 0 ≤ x ≤ L, −H ≤ z ≤ 0, following the procedure described in (Turkington et al., 1991), with the initial guess obtained from WNL theory and A, the scaled available potential energy held fixed.The next iteration (given the current iteration η k and c k ) is computed as follows: where λ k = gH c k 2 , and 2. Compute λ k+1 by where and 3. Define the new estimate of isopycnal displacement by η k+1 according to and the new estimate of the wave propagation speed according to This procedure is repeated until convergence criteria are met.
The efficiency and accuracy of the solution procedure is thus determined by how accurately the multiple Poisson solves in step 1 can be carried out, and the integrals in step 2 can be computed.Note that neither the wave amplitude, nor the wave propagation speed are specified.Instead, the kinetic energy of the disturbance is minimized under the constraint that the available potential energy is held fixed.
We have discretized the Laplacian in Eq. ( 18) using second order finite differences, as well as Chebyshev pseudospectral and Fourier spectral methods.For the Chebyshev and finite difference cases, we have solved the resulting matrix problem problem employing GMRES with various preconditioning strategies.For moderate sized problems an incomplete LU preconditioner worked well, though this required considerable storage overhead for larger grids.Finite differences are thus an impractical strategy for problems requiring repeated solution with high accuracy.A sine transform based solver, on the other hand, allows for spectral accuracy even using moderate grids, and scales well for larger grids as it does not require any explicit matrix construction.It further allows for a simple and effective refinement strategy in which the solution is iterated to convergence and then the grid resolution is doubled.The zero-padded Fourier transform of the coarse solution is employed to initialize the fine solution, which is subsequently iterated to convergence.Finally, considerable optimization of the FFT algorithm is available, for example as implemented in MATLAB, allowing the solver to take advantage of built-in parallelization strategies.
It is possible to test the scalability of the solution of the algorithm in various ways.After some experimentation, we have settled on using the change in predicted propagation speed, c, as the number of degrees of freedom is increased.This is effectively a Cauchy criterion for convergence.

Weakly nonlinear theory
It has been noted in the literature (Llewellyn Smith and Young, 2002) that the longwave eigenvalue problem Eq. ( 9) can yield nearly identical sets of eigenvalues for different background density profiles.For example, choosing a single pycnocline stratification as the target we find a number of candidate stratifications with a sum root mean square error less than 2 % for the first ten mode speeds, with the surface trapped density ρ(z) = 1 + 0.0601tanh z 0.07415H being one extreme example.However, numerical experimentation shows that a near match between longwave speeds is no guarantee that the weakly nonlinear coefficients r 10 and r 01 match, or in fact, are even close for mode-1 waves.We have conducted numerical experiments to determine if it is possible to find qualitatively different density profiles that yield nearly identical coefficients for weakly nonlinear theory.The experiments involved repeated (on the order of 10 000 simulations) solution of (Eq.9) using a Chebyshev pseudospectral method (Trefethen, 2000) while varying a set of parameters in the functional form for the density profile.The pseudospectral method is ideal for this type of application as it yields highly accurate eigenvalues and eigenspectra with a relatively low (< 200) number of grid points.Since we were interested in the first eigenmode and largest eigenvalue, only, the well-known issue of spurious roots in spectral discretizations (Trefethen, 2000) is not pertinent to the present calculations.Denoting the target set with the superscript T , we have used the error measure and have been able to find matches between several different sets of qualitatively different density profiles (see Fig. 2a for the example discussed in the following) to e < 0.01, well below what would be discernible from field measurements.Both of the density profiles shown have pycnoclines well away from the mid-depth and hence higher order nonlinearities leading to the Gardner equation are not expected to be important (Grimshaw, 1997).The spectral accuracy of the Chebyshev based method is vital to being able to perform the repeated calculations in a reasonable time (approximately 20 hours on desktop PowerMac).Comparable calculations with second order finite difference methods would take months.
With the two density profiles shown in Fig. 2, we have solved the DJL equation to compute two sets of fully nonlinear waves and to ascertain the nature of the upper bound on wave amplitude.The profile indicated by a dashed line in Fig. 2a has the wave amplitude bounded above by the onset of streamline overturning, while the profile indicated by the solid line in Fig. 2a has a wave amplitude bounded above by wave broadening to a limiting flat-centered wave (see Fig. 2b).The single pycnocline stratification (solid line) is labeled as "target" while the best fit profile, which exhibits significant stratification near the surface is labeled "fit".In Fig. 2b we show the ratio between the maximum waveinduced horizontal velocity and the wave propagation speed.The critical value of one, at which streamline overturning, or Nonlin.Processes Geophys., 18, 351-358, 2011 wave breaking, sets in, is indicated by a dot-dashed horizontal line.Weakly nonlinear theory with a quadratic coefficient of nonlinearity is another manner in which the behaviour of larger amplitude waves may be discussed and this is carried out in the Appendix.
We subsequently performed time-dependent simulations of the field Eqs.(1-3) using a variable time-step, secondorder projection technique described in, for example, Stastna and Lamb (2008), where the strengths and weaknesses of the model, along with typical resolution tests are described.The simulations considered forced flow over topography with a constant velocity u = (U,0) far upstream of the obstacle and were thus similar to past work on the resonant generation of ISWs (Stastna and Peltier, 2005).The numerical experiments were designed to demonstrate that for supercritical inflows (U = 1.32 c lw is used for Fig. 3) significant overturning occurs only for one of the two density profiles, even though they are nearly identical as far as weakly nonlinear theory is concerned.For subcritical inflows (U ≤ 1.2 c lw ) no breaking occurs, but the resonantly generated, upstream propagating waves are markedly different for the two cases (not shown).Figure 3a and 3b show six shaded density contours for the two cases, with the single pycnocline case in panel a.Both panels are for t = 250 where we have scaled time by the advective time scale H /U .In Fig. 3c and d we show Hovmöller plots leading up to the two states whose snapshots are shown in panels a and b, respectively.The Hovmöller plot is based on the wave-induced velocity at the surface, with darker colors indicating large positive values.It can be seen that for the single pycnocline case the topographically trapped disturbance is essentially steady after t = 150, while for the case with a significant near surface stratification breaking sets in around t = 70 and takes the form of billows ejected from the main wave body.

Two-scale exact solitary waves
Lamb and Wan (1998) have pointed out that for density stratifications with multiple pycnoclines it is possible to compute flat-crested waves of elevation and depression, as well as socalled conjugate flows that correspond to the limiting vertical profiles of isopycnal displacement, propagation speed and wave-induced horizontal velocity at the wave crest for either the positive or negative wave polarity.We consider a two pycnocline stratification For this stratification, two conjugate flows are possible, with conjugate flow speeds c + = 1.104c lw (c − = 1.164c lw ) for the waves of elevation (depression).Both polarities of ISW are bounded above by the so-called conjugate flow limit, of broad, flat-crested waves and hence regions of closed streamlines do not occur.To compute approximations of a multi-scale solitary wave, we first compute a flat crested ISW whose propagation speed matches the smaller of the two conjugate flow speeds (c + ) to eight decimal places.We subsequently compute a wave of depression, making sure that its propagation speed matches that of the flat-crested wave to eight decimal places.This wave will not be flat-crested since c − > c + .The resulting isopycnal displacement and wave-induced profiles are then used to construct the twoscale solitary wave.Note, that by construction, the density profile Eq. ( 25) does not specify the density profile far upstream.Indeed the upstream profile, as well as a background shear current profile, is specified implicitly from the wave solution, via the conjugate flow with speed c + .
In Fig. 4a we show ten isopycnals for a two scale solitary wave (only the right half of the wave is shown).In Fig. 4b we show vertical profiles of isopycnal displacement at the wave crest, or the left boundary of the figure (solid), and at the vertical dashed line in panel (a) (dashed) with respect to the far upstream density field (far right of the figure).Note that at the vertical dashed line in the figure, the horizontal, waveinduced velocity is constant.The double-humped wave was used as an initial condition in the time-dependent solver and was found to propagate without changing form in a stable manner (not shown).This is quite unlikely to be the case for double-humped waves computed from asymptotic theory, such as those in Makarenko et al. (2009).Moreover, Makarenko et al. (2009) suggest that their Fig.7 is similar to the measurements of Duda et al. (2004).However, the measurements involve such long time scales that the Earth's rotation is almost certain to play a role in the wave shape and evolution.
Again the spectral accuracy of the solver used for the DJL equation is absolutely vital.The repeated solution of the DJL equation necessary to find a wave with a matching propagation speed would involve prohibitively large grids were low order finite difference methods used.While results matching to eight decimal places are used to create Fig. 4, matching to as high as ten decimal places proved possible in less than two hours of computation time on a desktop iMac.

Conclusions
We have implemented solvers with spectral accuracy for both the weakly nonlinear and exact theories of internal solitary waves.While the primary role of these solvers is to provide spectrally accurate initial conditions for time dependent problems of ISW propagation and degeneration, we have used these models to derive two new mathematical results.First, we have demonstrated that weakly nonlinear, KdV theory for internal solitary waves exhibits a peculiar degeneracy in the sense that two very different background density profiles can yield nearly identical coefficients in the KdV equation.This means that weakly nonlinear theory would predict waves with the same structure and propagation speed.For the two density profiles constructed, the fully nonlinear solitary waves, which are exact solutions of the full stratified Euler equations, have a qualitatively different upper bound on amplitude.In particular, breaking and the possible formation of trapped cores occurs in only one case.This is consistent with time dependent simulations of wave generation by flow over topography in which only one of the two stratifications yields significant wave breaking.
The search algorithm involved repeated solution of the eigenvalue problem Eq. ( 9).Using standard low order finite difference algorithms, the search with even a moderate matching accuracy would take prohibitively long.Using Chebyshev pseudospectral methods it is possible to attain six digit accuracy in the eigenvalue and eigenfunction with tens of grid points.Hence the matching procedure can be carried out in a reasonable amount of time even for searches that require on the order of 10 4 iterations.
Second, for fully nonlinear ISWs, which are solutions of the DJL equation, we were able to compute exact double humped solitary waves, generalizing the asymptotic results of Makarenko et al. (2009).This procedure involved the matching of the ISW speed to a preset accuracy (eight digits in the above, though tests with ten were successfully carried out).Again, the feasibility of this procedure depended on getting an accurate solution for grids with moderate resolution.In the above, this was accomplished using a numerical method that employs the fast sine transform.Doublehumped ISWs involve a delicate balance of wave speeds and density stratifications with multiple pycnoclines.Moreover, the solution procedure for such waves (unlike the standard algorithm for internal solitary waves) determines the upstream density and background current profiles implicitly.Future work should thus address whether an efficient means of generation for these exotic mathematical objects occurs in the ocean.
The utility of spectral methods for internal wave computations is, of course, not restricted to exact solitary waves.Indeed, as mentioned above, the primary utility of the above described methods is in providing initial conditions for time dependent simulations.We are presently involved in the construction of a pseudospectral model with bottom topography, which will employ the above described methods for initialization, and both the model and novel results made possible by it, will be reported on in the near future.

Appendix A
Following the notation of Helfrich and Melville (2006), the Gardner equation is written as This equation is not the general, full second order extension of the first order KdV since it neglects both second order dispersive, and nonlinear-dispersive terms.It is the correct first order governing equation when α 1 is much smaller than the dispersive term and, for a single pycnocline stratification, this occurs when the pycnocline is centered near the middepth.For the degenerate KdV theory example presented in the main text α 1 and β are the same for both stratifications.Indeed the "target" stratification is chosen so that a priori α 2 should be neglected in a formal first-order scaling argument.
Nevertheless, an a posteriori check reveals that α 2 is different for the two stratifications (by 52.4 %).Note that as pointed out in Lamb and Yan (1996) (Eq.3.17 and the related discussion) the computation of α 2 requires an arbitrary multiple of the linear eigenfucntion to be specified (α 1,0 in their notation).We chose α 1,0 = 0 though it is far from clear that this is the best choice.
It is possible that a larger optimization problem in which both the background density, ρ(z) and current, U (z) are varied could yield a situation in which α 1 , α 2 and β all match.An efficient numerical formulation and solution of this problem remains as an outstanding problem for future research.
The solitary wave solutions of Eq. (A1) yield a flat-crest, or table top, limit (assuming that α 2 has the appropriate sign).However, the amplitude of this limiting solution generally does not match the limiting exact solitary wave (unless the pycnocline lies near the mid-depth).This is evident from Fig. 5 in Helfrich and Melville (2006) which compares the propagation speed and wave half-width for the KdV, Gardner and fully nonlinear-weakly nonhydrostatic Miyata-Camassa-Choi (MCC) theory.The target stratification is similar to the two leftmost panels.
Whether the Gardner theory would yield the same result as that found via DJL theory in the main text is thus far from clear.It is true that the α 2 coefficient in the Gardner theory is different for the "target" and "match" stratifications, but whether the maximum wave-induced velocity exceeds the propagation speed, and the wave breaks (before broadening to the limiting wave occurs) would require the vertical velocity profile to be worked out.This is a non-trivial task in the second order theory, as discussed in detail by Lamb and Yan (1996).
On the other hand the spectral implementation of the DJL equation solver presented above can answer the same question with no ambiguity, and given the efficiency of our implementation, in a manner of minutes.

Fig. 1 .
Fig. 1.The change in the predicted propagation speed as a function of degrees of freedom.Chebyshev pseudospectral method -solid, Fourier sine transform -dashed, Second order finite difference -dot-dashed.

Fig. 2 .
Fig. 2. Left Panel: The scaled N 2 (z) profiles for the two profiles that yield nearly identical coefficients in the KdV theory.Right Panel: The ratio of maximum wave-induced horizontal velocity to the wave propagation speed as a function of dimensionless isopycnal displacement for the two cases shown in the Left Panel.Waves break when u/c exceeds one (the dot-dashed, horizontal line).

Fig. 3 .Fig. 4 .
Fig. 3. (a) Six shaded isopycnals for the supercritical response over the topography for the case labeled "target" in Fig. 2. The flow is nearly steady for this case.(b) Six shaded isopycnals for the supercritical response over the topography for the case labeled "fit" in Fig. 1.The flow is dominated by a breaking region near the surface, and is highly unsteady for this case.(c) Hovmöller plot of the wave-induced horizontal velocities for the case pictured in (a).Stronger upstream directed velocities are indicated by darker regions.(d) Hovmöller plot of the wave-induced horizontal velocities for the case pictured in (b).Downstream propagating billows are visible as alternating bands emanating from around x = 0.