A model for large-amplitude internal solitary waves with trapped cores

Large-amplitude internal solitary waves in continuously stratified systems can be found by solution of the Dubreil-Jacotin-Long (DJL) equation. For finite ambient density gradients at the surface (bottom) for waves of depression (elevation) these solutions may develop recirculating cores for wave speeds above a critical value. As typically modeled, these recirculating cores contain densities outside the ambient range, may be statically unstable, and thus are physically questionable. To address these issues the problem for trapped-core solitary waves is reformulated. A finite core of homogeneous density and velocity, but unknown shape, is assumed. The core density is arbitrary, but generally set equal to the ambient density on the streamline bounding the core. The flow outside the core satisfies the DJL equation. The flow in the core is given by a vorticity-streamfunction relation that may be arbitrarily specified. For simplicity, the simplest choice of a stagnant, zero vorticity core in the frame of the wave is assumed. A pressure matching condition is imposed along the core boundary. Simultaneous numerical solution of the DJL equation and the core condition gives the exterior flow and the core shape. Numerical solutions of time-dependent non-hydrostatic equations initiated with the new stagnant-core DJL solutions show that for the ambient stratification considered, the waves are stable up to a critical amplitude above which shear instability destroys the initial wave. Steadily propagating trapped-core waves formed by lock-release initial conditions also agree well with the theoretical wave properties despite the presence of a “leaky” core region that contains vorticity of opposite sign from the ambient flow. Correspondence to: K. R. Helfrich (khelfrich@whoi.edu)


Introduction
Oceanic observations of internal solitary-like waves show clearly that waves of extremely large amplitude are possible, even common.For example, waves off the Oregon coast have amplitudes of up to 25 m where the upper layer depth is only 5 m (Stanton and Ostrovsky, 1998).The measure of nonlinearity α=a/H ≈5, where a is the wave amplitude and H is a depth scale (say H =h 1 , the mean upper layer depth), is well beyond the KdV limit α 1.While the Oregon observations are at the extreme end, values of α=O(1) have been observed in numerous other locations such as the South China Sea (Orr and Mignerey, 2003;Duda et al., 2004) and the New Jersey shelf (Shroyer et al., 2009).
An interesting and important aspect of such large waves is the possibility of trapped, or vortex, cores.Klymak and Moum (2003) and Scotti and Pineda (2004) found largeamplitude solitary-like waves of elevation with trapped cores propagating along the ocean bottom.The Morning Glory in northeast Australia is an atmospheric example of a largeamplitude internal undular bore that has been observed to have trapped cores beneath individual crests (Clarke et al., 1981).Doviak and Christie (1989) and Cheung and Little (1990) report similar trapped-core waves propagating in the atmospheric boundary layer.Both of these observations are significant since they also have coincident measurement of potential temperature that show that the cores contained fluid with densities near those of the ground-level ambient environment.
The first experimental study of internal solitary waves with trapped cores were reported by Davis and Acrivos (1967).They produced mode-two solitary waves with recirculating cores riding on a thin interface between two uniform density layers.Maxworthy (1980), Stamp and Jacka (1995) and, more recently, Sutherland and Nault (2007) have studied Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
experimentally the properties of mode-two solitary waves with trapped cores.However, there is a subtle, but important issue with the mode-two waves.Regardless of whether a trapped core is present, localized, stationary mode-two solitary wave solutions can not be found because of a resonance with finite-wavenumber mode-one waves (Akylas and Grimshaw, 1992).From a physical point of view, these mode-two waves do exist, but only for a finite time before the mode-one radiation destroys them.Grue et al. (2000) and Carr et al. (2008) carried out experimental studies of large-amplitude, mode-one internal solitary waves of depression.The waves were generated by a partial-depth lock-release into a two-layer ambient stratification.The lower layer had uniform density, the upper layer was linearly stratified and the fluid released by the lock has an upper layer with density equal to the density of the free surface.The focus of these experiments was on the presence of instability, both shear and convective, in the waves.However, they also found some waves with regions of fluid velocity u>c, i.e. trapped cores.
Despite the atmospheric, oceanic and laboratory observations of large-amplitude internal solitary waves with trapped cores there remain many unaddressed issues.Foremost among these is the need for a physically consistent theory for waves with trapped cores.As discussed below, models for waves with trapped cores have, with one exception, been based on a questionable approach.
In (1) the buoyancy frequency N 2 = −(g/ρ 0 )d ρ/d z is found from the ambient density distribution ρ(z).The density of the fluid at any point in the flow is ρ(x,z) = ρ (z − η(x,z)).
The wave is held stationary by a uniform oncoming flow u=−c.Here c is the phase speed of the wave and the horizontal and vertical velocities in the frame of the wave are, respectively, u=c(η z − 1) and w=−cη x .The total streamfunction ψ=c (η − z).
In what follows, the variables are non-dimensionalized using the depth H for x and z and g H for u, w, and c.Here g = g(ρ b − ρ t )/ρ 0 is the reduced gravity based on the densities at the bottom, ρ b , and top, ρ t , of the ambient stratification.The densities are scaled by a reference density ρ 0 and N 2 by g /H .
Mode-one solitary wave solutions of (1) are found for phase speeds c>c 0 , where c 0 is the linear long-wave phase speed (found by solving (1) for ∂/∂x = 0 and N 2 = N 2 (z)).The wave amplitude increases with c until the solutions end in one of three outcomes (Lamb, 2002).One is a low Richardson number shear instability limit.If the shear limit does not occur the solutions may reach the limiting flat-top wave, or conjugate state, solution.In some cases the solutions reach a breaking limit defined by the presence of incipient overturning η z =1, or u=0 in the frame of the wave, at z=0 (z=1) just beneath the wave crest (trough) for waves of elevation (depression).The phase speed at the point of incipient breaking is denoted c * .A necessary condition for overturning is that the ambient stratification have finite N 2 at the bottom (top) for elevation (depression) waves (cf.Lamb, 2002).In the following only waves of elevation will be considered since the fluid is assumed to be Boussinesq.Waves of depression can be obtained by symmetry.
Solutions to the DJL equation can be found beyond the critical point (c>c * ) if one assumes that the function for ρ(z) continues to increase smoothly for negative arguments (i.e., outside the physical domain).That is, N 2 (z − η) is defined for z = z − η < 0 (e.g.Davis and Acrivos, 1967;Tung et al., 1982;Brown and Christie, 1998;Fructus and Grue, 2004).These solutions develop closed recirculation zones (trapped cores) within a finite region bounded by the bottom and the core boundary h(x), |x| ≤ L c , shown schematically in Fig. 1.Because the streamlines within this core have z − η < 0, they are nominally linked to streamlines that "originate" beneath the bottom and have densities outside the ambient range.This is a consequence of the violation of the assumption in (1) that all streamlines in the fluid extend to x=±∞.Furthermore, trapped-core solutions found by extending N 2 for z−η<0 have the same relation between vorticity, stratification and streamfunction inside the core as in the exterior flow.While extension of ρ for z<0 is mathematically acceptable and solutions to (1) can be found, the presence of densities outside the ambient range, statically unstable core density structure, and the core vorticity-streamfunction relationship from the ambient upstream flow make these solutions physically questionable.
One solution to this problem is to specify a uniform core density ρ c ≥ρ b .This follows from Derzho and Grimshaw (1997) who derived a nonlinear, long-wave theory for a background stratification with N≈constant and with ρ c =ρ b .This ambient stratification is a special case since the maximum height of the trapped core, h c =h(0), remains small, permitting analytical progress.For more general stratifications there is no guarantee that the core height will remain small.
In this paper a model of internal solitary waves with arbitrary stratifications with trapped cores is developed.Once c>c * and a trapped core forms, the model consists of solving the DJL equation outside the core and matching this solution to another model of the core structure.We will follow Derzho and Grimshaw (1997) and assume a finite core with uniform density ρ c .For generality, we allow ρ c ≥ρ b , but the physically consistent choice for steady motion at long times in flows with weak diffusion is for the core density to homogenize to ρ c =ρ b by the Prandtl-Batchelor theorem (Batchelor, 1956;Grimshaw, 1969).
In general, the flow in an inviscid recirculating core with uniform density obeys (Batchelor, 1956) where f (ψ) is some unknown function.Frequently in problems of this type, the Prandtl-Batchelor theorem is used to argue that vorticity within the closed recirculation region homogenizes to the value of the vorticity of the ambient fluid on the bounding streamline (e.g., Rhines and Young, 1982).However, from (1), the vorticity of the ambient fluid flowing along the core boundary (η=h(x)), is not constant and the Prandtl-Batchelor argument can not be invoked.Just inside the homogeneous core, the vorticity on the bounding streamline is constant regardless of f (ψ) since ψ is constant on a streamline.Thus the vorticity in any inviscid, uniform density core model will be discontinuous across the core boundary.
In the Derzho and Grimshaw (1997) model the core is shown to have zero vorticity (i.e. was stagnant in the frame of the wave) to leading order and the vorticity discontinuity could also be ignored.This was confirmed in subsequent numerical studies (Aigner et al., 1999;Aigner and Grimshaw, 2001).However, with arbitrary N 2 (z) and inviscid dynamics there is no similar theoretical limit on the core circulation.In principal any f (ψ) is possible.Consideration of the formation process, viscous effects and stability may constrain the possible solutions.
To make progress we will make the assumption that the core has zero vorticity.In the reference frame moving with the wave phase speed c, the core is stagnant.The full model for internal solitary waves with stagnant, uniform density trapped cores is given in the next section where the matching conditions on the core boundary are given.The numerical procedure for finding the steady solutions is presented in Sect.3. The solitary wave solutions are explored for a particular class of ambient stratification that lead to trapped cores in Sect. 5.Because the stagnant-core assumption is just one of many possible choices, these new model solutions are compared to two-dimensional numerical solutions of the time-dependent, non-hydrostatic Euler equations (described in Sect.4) initiated with the new theoretical solutions and to solitary waves produced by a lock-release initial condition.

The model
As already stated, the flow in the ambient fluid outside the core is governed by the DJL equation ( 1) with the boundary conditions When there is no core (c<c * ), h(x) = 0 (and L c =0) so that the last condition is redundant.When a core is present, the last condition is a statement of the kinematic constraint that fluid parcels flowing along z=0 upstream of the core must remain adjacent to the core boundary.If h(x) were a topographic feature, rather than a internal fluid boundary requiring dynamical considerations, then the problem statement would be complete.
In general, inside the core (z≤h(x) for |x|≤L c ) the flow is governed by (2) with ψ=0 on the unknown boundary z=h(x).Along the boundary the core and ambient pressures must be equal.Consider the ambient flow streamline that originates upstream at x=∞ along z=0.The dimensional Bernoulli constant on this streamline is where is the (hydrostatic) pressure at z=0 as |x|→∞.The pressure on the ambient bounding streamline is then given by The subscript b is used to denote this ambient streamline that originates upstream at z = 0 where ρ=ρ b .
The Bernoulli function on the bounding streamline in the core is equal to the pressure at the stagnation point x=L c Thus the pressure on this bounding streamline is The subscript c indicates core quantities and u c and w c are the core fluid velocities in the general case where the core fluid is recirculating.The pressure must be continuous at the core boundary z = h(x) (for |x|<L c ).Thus from ( 5) and ( 6) the matching condition is This can be simplified using the kinematic boundary condition w=uh x on z = h(x) for both the core and ambient flow to give When ρ c =ρ b , the horizontal velocity is continuous at the core boundary.
With the assumption of a stagnant core, u c =0, (8) becomes, after using the non-dimensionalization introduced with (1), or Here S = (ρ b − ρ t )/(ρ c − ρ t ) ∈ [0, 1].The negative root was chosen so that u b ≤0, avoiding overturning in the free stream.Note that when ρ c =ρ b , S = 1 and (10) give η z =1.In this case the velocity in the exterior flow adjacent to the core is zero (in the wave frame) and just at the overturning limit.For ρ c >ρ b , u b <0 along the core.
When c 0 ≤c≤c * there is no core and the DJL equation must be solved for η(x,z) subject to (3) with L c =0.For c>c * the unknown core boundary h(x) is found by solution of (1) subject to (3) along with the core condition (10).These stagnant-core solitary wave solutions should exist from c * up to the limiting conjugate state solutions with stagnant cores found by Lamb and Wilkie (2004) for ρ c =ρ b and White and Helfrich (2008) for ρ c ≥ρ b .The conjugate state solutions are found in a one dimensional version (∂/∂x=0) of the DJL problem, along with additional conditions that impose energy and momentum flux conservation between the flows upstream and over the uniform section of the wave.The conjugate states have core heights h cs and speeds c cs .Therefore the solitary wave solutions have core heights 0≤h c ≤h cs for c * ≤c≤c cs .
3 Stagnant-core DJL solution method Turkington et al. (1991) developed an efficient variational technique that is frequently employed for solving the DJL equation in the absence of a core (c<c * ), or when the ambient stratification is extended for z − η<0 (c>c * ).However, this method is not easily adapted to the current problem with uniform density cores.Numerical solutions will instead be found using the Newton-Raphson technique (Press et al., 1986).In order to solve the DJL equation outside the trapped core, which is for the moment assumed to be known, (1) is re-written in a coordinate system where the physical coordinates (x,z) outside the core are mapped to a new coordinate system σ =σ (x,z) and ξ = ξ(x,z) (i.e., x=x(ξ,σ ) and z=z(ξ,σ )).With this transformation (1) becomes where and is the Jacobian of the transformation.The subscripts indicate partial differentiation.Here a standard boundary-fitted system is employed.The function (x) is introduced to allow a stretched grid in x.
The mapped DJL equation ( 11) is solved in a domain 0≤σ ≤1 and 0≤ξ ≤1.Care is taken that L, the domain length in x, is sufficiently large to minimize effects of a finite length domain.The wave is symmetric about the crest so that η x =h x =0 at ξ =0.From (3) and using η where it is understood that h(ξ )=0 for ξ ≥ ξ c = L c /L. Finally, in the mapped coordinate system the core pressure matching condition (10) becomes, using η Nonlin.Processes Geophys., 17, 303-318, 2010 www.nonlin-processes-geophys.net/17/303/2010/ Treating the core boundary z = h(x) as known (in each iteration of the Newton method), ( 11)-( 16) are approximated using second-order finite-differences.The grid has σ discretized into M+1 uniformly spaced points in σ and N+1 uniformly spaced points in ξ , with a ξ -grid point located at ξ c .The core boundary h i is then defined at N c grid points where h>0.The discretized system (11)-( 16) results in (M − 1)(N − 1) + N c equations for the interior values η ij (i = 2 − N, j = 2 − M) and the N c non-zero values of h i .This set is solved by the Newton-Raphson method.After a Newton update, the new h i (i≤N c ) will have changed and a new estimate of the end point ξ c is found by a quadratic extrapolation of the last three (i = N c − 2 to N c ) core points.The solution is then interpolated onto a new grid with a grid point located at the new L c .The ξ grid is arranged such that there is approximately the same grid spacing inside and outside the core (i.e.ξ ≈ 1/N everywhere).Thus the number of points, N c , that define the core varies as the core length changes.For very small cores care is taken that N c >3.
This solution and remapping procedure is repeated until the sum of absolute value of the variable increments and the sum of the absolute values of function residuals each drop below specified tolerances.Particular care is taken that the core length L c and height h c (=h(0)) are not changing.For solutions discussed here M=50 and N ≈50 L. Thus in the absence of a core the grid is approximately isotropic.In the presence of a core the grid provides sufficient resolution in x and z to resolve the core stagnation point at L c .This solution procedure is not sophisticated; however, it does produce solutions that, as will be shown below, remain essentially unchanged when used as initial conditions in time-dependent, non-hydrostatic numerical calculations.The exceptions are waves that appear to be unstable due to a physical shear instability.
A trapped core solution c>c * found using the Newton-Raphson procedure with N 2 (z − η) extended for z − η < 0 provides an initial guesses for h(x), L c , and the exterior η(x,z) for the same stratification and c.Once a converged solution is obtained, it is used as the initial guess for a new c, leading to a family of solutions for c>c * with a given ambient stratification.

Time-dependent non-hydrostatic numerical model
The theoretical solutions obtained with the trapped-core model will be tested using numerical solutions of the inviscid, two-dimensional non-hydrostatic equations of motion Here u and w are the horizontal and vertical velocities in the laboratory frame, s = (ρ − ρ 0 )/(ρ b − ρ t ) and p is the pres-sure (less the hydrostatic pressure from ρ 0 ).The u, w, x, and z have been non-dimensionalized as (1).Pressure is scaled with ρ 0 g H and time t with (H /g ) 1/2 .These equations are solved using the finite-volume, second-order projection method of Bell and Marcus (1992).Because of the upwindbiased Godunov evaluation of the nonlinear terms the method is stable, requires no explicit dissipation, and introduces signicant numerical dissipation only when large gradients occur on the grid scale.The numerical method has been used successfully for studies of large-amplitude internal solitary waves (e.g., Lamb, 2002;White and Helfrich, 2008), and tests of the code show that large internal solitary waves without trapped cores propagate distances of more than 100H with the correct phase speed, amplitude, shape, and minimal loss of energy (< 1%).Two types of calculations are considered.In the first, the stagnant-core DJL solitary waves are tested by using these steady solutions as initial conditions.In the second, trappedcore waves are generated by a lock-release initial condition

Results
To illustrate the solutions to the DJL model for largeamplitude waves a background density given by with is used.This form for s(z) gives bottom-trapped stratification for increasing λ (≥ 0) and N 2 (0) > 0 so that trapped cores are possible. www

Extended-s DJL waves
We first examine waves obtained by extending the ambient N 2 (z) for z < 0 and call those solutions "extended-s DJL" waves.The time-dependent evolution of an extended-s DJL wave for λ=8 and c=0.4 is shown in Fig. 2. Figure 3 shows a larger wave with c=0.58 and the same λ=8.The top panel of each figure shows the theoretical solutions with recirculating (negative vorticity) cores with densities greater than the ambient density along the bottom, s(0)=1 and statically unstable core density distributions.The amplitude η M , core height h c and core length L c of these extended-s DJL solutions for λ = 8 and 4 are shown by the dashed lines in Fig. 8.The wave amplitude η M is defined as the maximum streamline displacement of the flow outside the core at the wave crest.For λ = 8 (4), the linear long-wave phase speed c 0 =0.226 (0.285), the critical overturning speed c * = 0.331 (0.383) and the extended-s DJL conjugate state speed c cs = 0.582 (0.505).
The Richardson number, Ri, of each of these waves falls between 0 and 0.25 in a region in the upper, stably-stratified portion of the core that extends in thinning bands to the two stagnation points as shown in Fig. 4. Recent work by Carr et al. (2008), Fructus et al. (2009), and Barad and Fringer (2010) has shown that internal solitary waves are unstable to Kelvin-Helmholtz instability when the minimum Richardson number is less than about 0.1 and the ratio l Ri /l W ≥0.86.x z (a) Here l Ri is the horizontal length of the Ri≤0.25 region and l W is the solitary wave width (twice the distance from the crest to where η of the maximum amplitude streamline equals 0.5η M ).The latter condition is related to the time (or distance) required for the unstable waves to grow before they are swept out of the unstable region and left behind the solitary wave.However, the extended-s waves are different from the situations considered in Fructus et al. (2009) and Barad and Fringer (2010) in two important ways.The low Ri region is confined to the recirculating core where disturbances could re-enter the unstable region.Furthermore, streamlines in the low Ri region also pass through the gravitationally unstable lower core.
In the time-dependent solution for the c = 0.4 wave in Fig. 2 the trapped core does become unstable.The instability first appears as a wobble, or distortion, of the core center as shown at t = 30.By t = 70 the core begins to breakdown and density inversions are evident.Between t = 70 and t = 200 the flow stabilizes to a new trapped core wave with wave amplitude η M = 0.208 and speed c = 0.394 that are only slightly different from the original wave.However, the trapped core is smaller and the maximum core density has decreased from s = 1.192 to 1.092 and become more homogeneous.The origins of the instability are not clear in this case.The initial wobble could be the result of either gravitational or shear instability.
The instability and its consequences are much more dramatic for the larger, c = 0.58, wave in Fig. 3.This example also provides a clearer indication of the processes involved.At t = 7.5 the first disturbances have developed in the lower, gravitationally unstable part of the core.These disturbances grow and are swept up into the upper half of the core near the forward stagnation point.They then appear to excite Kelvin-Helmholtz billows near the wave crest in the low Ri band (t = 10).An analysis of the energetics of the disturbances at small times would be informative but is beyond the current scope of this work which is focussed on the stable core states.By t = 30 the instability has nearly destroyed the wave.Amazingly, between t = 50 and 100 a new trapped-core solitary wave forms from the disorganized flow.Between t = 100 and 200 the flow completely stabilizes to a slower wave with c = 0.558.
Figure 5 shows a comparison of the vertical profiles of u − c, s, and Ri at the crests of the initial (t = 0) and the equilibrated (t = 200) waves.The final core has a nearly homogeneous density core, s≈1.18, that is much less than the maximum value of s = 1.829 at t = 0.The horizontal velocity along the bottom is larger than the initial wave, indicating an increase in the magnitude of the core vorticity.The wave is stable even though Ri<0.25 in the upper part of the core (0.38 < z < 0.58).This unstable region extends to the two stagnation points similar to the unstable band at t = 0 in Fig. 4b.But unlike the initial wave, the uniform density of the final core inhibits the gravitational instability.It appears that shear instabilities are unable to grow sufficiently without  coupling with the gravitational instability.The low Ri region near z = 0.9 is a result of the vanishing stratification there and remained stable in the calculation.
The final core vorticity is the same sign as the initial core and ambient vorticities.A scatter plot of vorticity versus streamfunction within the core is shown in Fig. 6.A clear relationship, ∇ 2 ψ = f (ψ), has developed in the core.It is distinct from the broad range of core vorticity and streamfunction values at t = 0 (the light gray area in Fig. 6).The increase in the magnitude of the core vorticity following the instability and equilibration is clear from this figure.A similar, although different, relationship develops for the smaller wave in Fig. 2.
It is not possible to predict the vorticity-streamfunction relationships, f (ψ), of the equilibrated waves.They are dependent on the turbulent flow, baroclinic vorticity production, vorticity advection, and slight numerical diffusion that occur in the flows.The important point of these two examples is that these extended-s waves, at least for the ρ(z) considered, are unstable.What emerges from the instability are new trapped-core internal solitary waves that have nearly  homogeneous density cores with distinct relationships between core vorticity and streamfunction.Both of these outcomes support the theoretical model introduced in Sect. 2. However, the current restriction of the theoretical model to stagnant-core solutions is a clear limitation.Regardless, it is useful to explore the stagnant core solution properties, including stability, and ask if this model provides reasonable estimates of trapped-core waves that develop from a general initial-value problem such as a lock-release.

Stagnant-core DJL waves
An example of a stagnant-core solution for λ = 8, s c = 1, and c = 0.4 is shown in Fig. 7.This wave has the same ambient stratification and phase speed as the extended-s wave in Fig. 2. It has slightly larger amplitude, smaller core height, and larger core length than the extended-s wave.Another difference is that the stagnant-core wave has a region of Ri<0.25 above the core (bounded by the heavy in 7).The Ri=0.154 and l /l W 0.37.Recall the extended-s examples had the of Richardson located entirely the core (with exception of values in the very weakly region near upper cf.5).The contour intervals are 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that only the lower half of the domain is shown.
the Turkington et al. (1991) solution method.For λ = 8 (4) the Richardson number first drops below 0.25 at c = 0.368 (0.395) and by the last solution found the minimum Ri above the crest is 0.073 (0.115) and the ratio l Ri /l W =0.8 (0.38).
From Fructus et al. (2009) and Barad and Fringer (2010), these waves should be stable to shear instability.Figure 9 shows the time-dependent evolution of the stagnant-core DJL wave in Fig. 7 (λ = 8 and c = 0.4).In this case, the initial stagnant-core wave is subject to an instability that originates along the core boundary in the low Ri zone.The Kelvin-Helmholtz billows are clear at t = 20.The instability is rather mild and between t = 50 and 100 the flow restabilizes to a wave that is only slightly different from the initial wave.Figure 10a-d shows the wave properties (s, u − c, w, and vorticity ∇ 2 ψ) from the numerical solution averaged between t = 90 and 100.The core density is slightly inhomogeneous.The instability has injected ambient fluid into the core and baroclinic vorticity production and advection produces a weak positive vorticity (the initial and ambient vorticity is everywhere negative) (see Fig. 10d).The flow in the core is very weak (see Fig. 10b and c), but does include a band of u − c > 0 located above the bottom.In contrast, in the extended-s solutions the region of u − c > 0 lies along the bottom (i.e.core vorticity is negative).The phase speed of the equilibrated wave c=0.399 is nearly the same as the stagnant-core model wave.The core boundary from the initial wave (heavy solid line in Fig. 10a-c) also agrees with the apparent core shape in the equilibrated wave.In particular, the region of near-zero vertical velocities in Fig. 10c lies entirely with the initial core boundary.Further agreement between the final wave and the initial stagnant-core model wave is shown in Fig. 10e where the vertical profiles of u − c and s at x − ct = 0 (appropriate c values used) from the theory and equilibrated wave in the numerical model are compared.The agreement is very good.The only significant differences are within the core, where the weak positive vorticity is clear, and s ≈ 0.97.The latter is close to the value s = 0.973 in the bottom cell (0 ≤ z ≤ 1/150) in the ambient flow ahead of the wave.Thus, in the numerical solution the instability replaces the initial core fluid (s = 1) with dense ambient fluid.This is to be expected if the effect of the instability is a turbulent diffusion.The Prandtl-Batchelor theorem then says that the core  should equilibrate to s = 0.973, the value of the density flowing along the bounding streamline.The transport of ambient fluid into the core and the expelling of core fluid out the back of the wave is an example of Lagrangian advection across a hyperbolic trajectory (Wiggins, 2005).The shear instability along the core boundary provides the time-dependent forcing of the hyperbolic trajectory (the bounding streamline approaching the rear stagnation point), subsequent lobe formation and transport.
The equilibrated wave has a region of Ri<0.25 just above the core with minimum Ri≈0.15 and l Ri /l W ≈ 0.37, nearly identical to the initial wave.This suggests that the instability observed in the numerical solution could be a consequence of the initial vorticity discontinuity at the core boundary.One effect of the instability is to eliminate the vorticity discontinuity and reduce the vertical shear at the core boundary (see Fig. 10d).
Another example of the evolution of a stagnant-core DJL wave is shown in Fig. 11.As in the previous example λ = 8 and c = 0.4, but the core of the initial wave has density s c = 1.1.Recall from the boundary condition ( 8) that the jump in density across the core boundary gives a jump in velocity, which should enhance the likelihood of shear instability.The time-dependent solutions show that large overturns develop at the core boundary by t = 20.However, by t = 40 the instability has weakened considerably and by t = 100 the wave has stabilized.This final wave propagates at a speed c = 0.394 and has a core that is nearly stagnant and homogeneous as shown by the dashed lines in Fig. 11e.The mixing from the instability has again replaced the initial core fluid with s = 0.97 fluid.Thus the new wave is close to a stagnantcore DJL wave with s c = 1.The solid lines in Fig. 11e show the vertical profiles of u−c and s from the stagnant-core DJL solution with c = 0.394.With the exception of the weak core circulation and the slightly lower density of the core fluid, the equilibrated wave agrees very well with the theory.
Time-dependent numerical solutions initialized with stagnant-core DJL waves with s c = 1 and other values of c for both λ=4 and 8 show similar behavior.An initial period of shear instability is followed by flow stabilization to a large-amplitude wave with properties very close to those predicted by the stagnant-core DJL model.The amplitude η M , core height h c , and core length L c of these equilibrated waves for λ = 4 and 8 are shown in Fig. 12 (by the open circles).The equilibrated wave from Fig. 11 is indicated by the open diamonds in Fig. 12a-c.The values h c and L c are determined from the vertical and horizontal extent of the region of near-zero vertical velocities (i.e. the w = 0 contours) as shown in Fig. 10c.There is some uncertainty in choosing these values that is reflected in the slight scatter of the data.However, the agreement between the stagnant-core DJL theory and the numerical results is quite good.
These results support the stagnant-core theory despite the fact that the full numerical solutions produce cores with weak circulation with positive vorticity and nearly homogeneous density.However, this is not too strict a test of the theory as it says that if flow starts close enough to the stagnant-core solution, it will remain near it.The calculations shown in Figs. 2 and 3, on the other hand, show that solutions with nearly homogeneous density, but vastly different core circulation are possible.In those cases the wave properties η M , h c , and L c do not agree with the stagnant-core model (not shown).
Additional tests of the stagnant-core theory were conducted with full time-dependent numerical solutions initialized with the lock-release initial condition (21) for λ=4 and 8. Figure 13 shows one example for λ = 8, h d = 1, L d = 0.6 and ρ d = ρ b (s d = 1).One large solitary wave emerges from the release.Between t = 20 and 30 the wave is still evolving, but propagates with a speed c≈0.435, which is greater than the critical speed for core formation c * = 0.331 for λ = 8.Contours of s are shown at the indicated times.The contour interval is 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that part of vertical domain is shown for t > 0.
A close-up of the wave at t = 300 is shown in Fig. 14.This long-time evolution is found by taking the leading wave from the lock-release run at t = 30 as the initial condition for a run with inflow-outflow boundary conditions with a uniform inflow velocity magnitude equal to the propagation speed c = 0.435 from the lock-release run.The wave continues to adjust until about t = 180 after which the phase speed, c = 0.417, is constant.This final wave has a well-defined core region from the w = 0 contour (Fig. 14c).There is a recirculation cell shown by the streamfunction in Fig. 14d and the horizontal velocities u − c > 0 (Fig. 14b).However, the recirculation cell (with positive vorticity) sits above the bottom.Streamlines from the ambient flow ahead of the wave split to flow around this recirculation cell.This is an instantaneous picture of the flow.The flow over the core is still weakly unstable (see the cat's eye of the Kelvin-Helmholtz billows on the trailing side of the wave), inducing a slow, continuous transport of ambient fluid into the core, and expelling core fluid out the back of the wave.This results in a nearly uniform density in the recirculation cell of s ≈ 0.88.The densest ambient fluid flows beneath this cell and is less involved in the cross-boundary transport than the (less dense) ambient streamlines that pass over the cell.
This "core" is quite different than those in either the stagnant-core or extended-s DJL models.There are no stagnation points along the bottom boundary.The density within the w ≈ 0 region (Fig. 14a) is not uniform and closely reflects the streamfunction shown in panel d.We note that a similar circulation pattern was found in numerical solutions of large-amplitude internal solitary wave breaking over www.nonlin-processes-geophys.net/17/303/2010/ Nonlin.Processes Geophys., 17, 303-318, 2010 314 K. R. Helfrich and B. L. White: Trapped-core internal solitary waves slope-shelf topography (Lamb, 2002).Figure 14e shows the vertical profiles of u − c and s at x − ct = 0 from the wave at t = 300 and the stagnant-core theory prediction for a wave with c = 0.409 and λ = 8.The agreement above the core (z > 0.2) is very good.
From the u − c profile it can be seen that the vorticity in the core is positive and opposite sign from the vorticity in the ambient fluid.This positive vorticity is not a direct product of the initial adjustment immediately following the lockrelease.Rather, it is a product of vorticity produced in the unstable shear flow surrounding the core.Figure 15 shows a close-up of the vorticity structure at t = 10 from Fig. 13.At this time the vorticity within the developing core is nearly zero.Patches of positive vorticity from baroclinic produc- .tion are present in the unstable flow around the core.Over time, some of this positive vorticity fluid is entrained into the core.In general, the lock-release runs lead to cores with larger magnitude positive vorticity than the stagnant-core initial conditions.This appears to be a consequence of an initial wave that is more unstable to shear instability along the core boundary, and hence greater baroclinic vorticity production and subsequent entrainment of this vorticity into the evolving core.
A scatterplot of vorticity versus streamfunction from within the core is shown in Fig. 16.For this plot the core region is defined to lie below the arc from x − ct = ±0.48through the center of the chain of cat's-eye vortices (which corresponds closely to the core defined by w = 0 in Fig. 14c).Within the closed recirculation cell where ψ ≤ −0.005, a clear streamfunction-vorticity relation has emerged with the positive vorticity values clustering along a linear trend with ψ.
By varying L c and h d (both ≤ 1) in ( 21), waves with a range of amplitudes and speeds are produced.As shown by the open squares in Fig. 12a and d for λ = 8 and 4, respectively, the amplitudes and speeds of these waves agree quite well with the stagnant-core theory.Furthermore, if h c and L c are defined by the w = 0 contours (see above) the agreement with the theory is also good.If the fluid behind the dam has density s d > 1, the solitary wave that emerges will, after enough time, expel the dense fluid and have a leaky core with density s ≤ 1.Despite the significant differences between the waves produced via a lock-release and the stagnant-core DJL theory, the theory does a surprisingly good job of capturing the overall wave properties.
The lock-release runs produced solitary waves that for λ = 8 were always slower and smaller than the fastest stagnantcore DJL wave found and for λ = 4 were only slightly faster than the maximum theoretical solution.This was the case despite initial conditions that have a larger volume of fluid of density ρ b behind the dam than the core volume of the largest wave.The available potential energy of the initial state is also greater than the total energy of the largest theoretical solitary wave.White and Helfrich (2008) did find conjugate state solutions (indicated by the open triangles) from lock-release runs with L d 1 and h d → 1 when the initial available potential energy per unit length of the region behind the dam was greater than the total energy per unit length of the (infinitely long) conjugate state.The leading face of the conjugate state waves were smooth, but intense Kelvin-Helmholtz billows formed on the uniform region of the conjugate state wave crest (see Fig. 17b in White and Helfrich, 2008).
The wave produced in the lock-release example in Fig. 13 appears to be more unstable than the waves in the laboratory experiments of Grue et al. (2000) and Carr et al. (2008).One possible reason for this difference is that initial dam height h d = 1, which is the maximum possible, is larger than the laboratory experiments (Grue et al., 2000 used h d ≈ 0.5 in experiments that produced large trapped-core waves).Thus the initial state in the model has more available potential energy.The stratification is also different.In the laboratory experiments an equivalent λ ≈ 6. Lastly, the two-dimensionality of the numerical model prohibits the transverse breakdown of the Kelvin-Helmholtz billows.Detailed wave properties will depend on whether the third dimension is included, especially if a careful comparison is to be made with the experiments.However, the calculations still provide insight into wave development and a systematic test of the theory.
One last example of wave evolution is shown in Fig. 17 for a case initiated with the largest stagnant-core wave found with λ = 8 and s c = 1 at c = 0.45 .As in the earlier examples the wave is initially unstable to shear instability localized along the core boundary.Unlike the previous examples, the instability does not weaken.Rather, the instability strengthens between t = 60 and 90 such that the wave continually diminishes in size and speed until t = 200 when the wave reached the downstream end of the domain and the calculation was stopped.If the domain were longer, the flow might eventually equilibrate to a new wave as in the earlier examples, but any final wave will be very much smaller 17.Evolution of the largest stagnant-core DJL solution found for λ = 8 and s c = 1 at c = 0.45.Contours of s are shown at intervals of 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that only the lower portion of the domain is shown.
than the initial wave.Recall that the initial wave has a minimum Richardson number above the core of 0.073 and that l Ri /l W = 0.8.These values put the wave closer to the instability boundary found by Fructus et al. (2009) and Barad and Fringer (2010).It is also consistent with the suggestion that the failure to find stagnant-core DJL solutions for c > 0.45 is related to flow instability and helps explain why the lockrelease runs produced solitary waves within the range of stagnant-core solutions.

Conclusions
In this paper we have developed a model of large-amplitude internal solitary waves with trapped cores that avoids the questionable approach of earlier models that use the ambient stratification-vorticity-streamfunction relation in the recirculating core.The new model consists of a finite size, homogeneous density core beneath a stratified exterior.The exterior flow satisfies the usual Dubreil-Jacotin-Long equation for fully nonlinear internal solitary waves and is matched to the core flow through a pressure continuity condition.The flow within the core must satisfy a relation between streamfunction and vorticity that can be specified arbitrarily in this www.nonlin-processes-geophys.net/17/303/2010/ Nonlin.Processes Geophys., 17, 303-318, 2010 inviscid theory.For to progress, have chosen the possible core flow -a stagnant steady frame of This density, stagnant-core model the incipient solitary wave to the limiting amplitude stagnant-core conjustate of Lamb and Wilkie and White Helfrich for steady, stagnant-core DJL internal waves obtained numerically and the were compared by extending ambient − for negative the stratifications considered, differences in the resulting wave propcould be substantial, especially as the wave amplitude We though, that if the ambient stratification is extended differently, then the disagreement can be reduced.For example, if s is extended with where N 2 0 = −s z(0) from the ambient stratification.This extension gives continuous density and N 2 at z = 0, and avoids numerical problems associated with discontinuous N 2 in the usual DJL solution methods.It asymptotes to s c for z → −∞, and for s c − 1 1, the core density will be nearly homogeneous with maximum density only slightly greater than the ambient environment and the recirculation will be weak.Wave properties computed with s c = 1.02 are nearly the same as the stagnant-core properties shown in Fig. 8.The solution branches even stop at nearly the same wave speed.This extension procedure, though not the form for s above, has been previously proposed by Brown and Christie (1998).While these special extended-s solutions are only slightly different from the new stagnant-core model, the latter provides a theoretical framework for developing homogeneous core solutions with specified circulation.
Time-dependent numerical solutions of the Euler equations initiated with the stagnant-core DJL solutions agree well with the theory, though the solutions exhibit shear instability along the core boundary and develop a weak recirculation (of opposite sign vorticity) within the core.The time-dependent solutions indicate that low Ri regions near the wave crest lead to wave instability that appear to explain the limited range of achievable wave amplitudes.
During the final stages of the submission process this paper, the authors were made aware of a recently submitted manuscript by King et al. (2010) that described a new numerical method for solving the DJL equation.They were able to compute waves with stagnant, uniform density cores using an extension of the ambient density structure similar to (24), but with much more rapid transition than above.Twodimensional, time-dependent numerical calculations initiated with both a stagnant-core wave and a wave found by extending the ambient density profile, were qualitatively similar to the results in this paper.
Solitary waves produced from a lock-release also agree reasonably with the stagnant-core model for wave amplitude, core height and core length.However, the core heights and lengths used in the comparison are defined, somewhat arbitrarily, from the vertical velocity fields.There is a recirculating core, but it has vorticity of opposite sign from the ambient flow and is off the boundary.Furthermore, the cores are not isolated from the ambient flow.Unsteadiness from shear instability gives rise to fluid transport into and out of the nominal core region such that the core density is less than the densest ambient fluid.As in the stagnant-core initiated runs, the positive vorticity in the core must be the result of baroclinic vorticity production as the ambient vorticity is negative.
Very different core circulations from those described above can occur, as illustrated by Figs. 3 and 6.In this example, the core is nearly homogeneous and has vorticity of the same sign as the initial condition, but with very a different vorticity-streamfunction relation.This emphasizes the role of transient wave formation, instability and dissipative processes on the final core structure and circulation.It also raises the question of whether there are general principles that place limits on the range of realizable core flows.Certainly, the presence of shear instability and the tendency of the flow to adjust to marginal stability is one such limit.This is an avenue for future investigation.
As mentioned in the Introduction, Grue et al. (2000) and Carr et al. (2008) carried out experiments on internal solitary waves produced from a lock-exchange and found that large waves could have regions of u > c.We have not made a detailed comparison with their experiments.There are, however, several aspects of the experiments that are in qualitative agreement with the results presented here.One is that in experiments with apparent trapped cores, the vertical profiles of horizontal velocity at the wave trough (the experiments produced waves of depression) show that u/c ≈ 1 in the upper water column as shown in Figs. 14, 15, and 18 in Grue et al. (2000) and Figs. 9 and 20 in Carr et al. (2008).It is not clear that this depth range represents the full height of a core since the authors do not provide independent estimates of the core height (e.g. from the density fields).However, the height of this zone at the wave center is approximately equal to the depth of the upper, linearly stratified layer (h upper ≈ 0.2H ).The wave amplitudes, measured by the vertical displacement of the interface between the two layers, are also approximately equal to the upper layer depth.From Fig. 7 it is clear that the height of the core is always less than the wave amplitude (defined somewhat differently in the theory), except for the largest waves where they are approximately equal.Thus the experimental core heights defined by the u ≈ c depth are at least qualitatively similar to the theory.Of course, the stratifications are different and this should be taken into account in any detailed comparison, but this is suggestive that u ≈ c in the wave cores observed in the experimental.
Another important feature of the experiments is that the cores are not completely stagnant.Small turbulent vortices are present both within the core and along the core boundary (e.g., Fig. 20 in Grue et al., 2000).And there is evidence of small, but finite, vertical shear in the cores.As the free surface is approached from below, the fluid velocity declines, producing a zone of vorticity within the core with opposite sign from the open portion of the flow (see Fig. 20 in Carr et al., 2008, although this may be an artifact of the no-slip lid in this particular experimental run).So while not conclusive, these experiments do lend support to the stagnant-core theory and numerical modeling results.Grue et al. (2000) also report that instability limited the maximum observed wave amplitude to be less than for the theoretical maximum wave (found using an extended-s DJL model).The conclusion that the stagnant-core waves are limited by instability is in qualitative agreement with these experiments.We note that in the experiments done by Carr et al. (2008) the fluid behind the lock was not seeded with particles (for PIV), yet the core regions in their Figs.2-7 are filled with particles (M.Carr, personal communication, 2010), indicating that ambient fluid is incorporated into the core.Unfortunately, neither study made direct measurements of the densities within the core region and these, along with experimental measurements of the trapped-core properties over a broad range of parameters, would be valuable.
Mass transport by large-amplitude internal solitary waves is thought to play an important role in cross-shore larvae transport (Pineda, 1999;Helfrich and Pineda, 2003).If the wave amplitude is large enough for trapped-core formation, the question of core structure and leakiness will be important.
streamline Fig. 14d that larvae near free surface (for waves of depression) incorporated into core, while residing the free surface could brought into a core a finite period of time experience significant transport during that period at the wave phase for possible transport fine suspended sediment elevation propagating the sea The observed Klymak and Moum (2003) show evidence acoustic backscatter from activity and optical backscatter from fine Finally, the time-dependent numerical calculations sented here are all two-dimensional.Clearly, the threedimensional evolution of the solutions needs to be explored.The effects of the gravitational and Kelvin-Helmholtz instabilities on the wave evolution and mixing will certainly be different in three dimensions, possibly leading to different core states.Fluid exchanges between the core and the ambient will likely be modified from the two-dimensional runs.Also, waves with recirculating (or stagnant) cores may themselves be unstable to transverse perturbations.

Fig. 1 .
Fig. 1.Sketch showing the core boundary h(x) and the streamline displacement η.The stagnation points are at x=±L c .
x<L d and z≥h d ρ d , x<L d and z<h d (21) The region of fluid held behind the lock extends from x = 0 to x = L d .Here h d is the depth of a lower uniform-density layer, ρ d ≥ρ b , beneath the ambient stratification.This initial condition mimics the experiments of Grue et al. (2000) and Carr et al. (2008).In both cases the solutions use uniform grids with vertical grid cell size of z = 1/150.The horizontal resolution x = 0.01 in the DJL model initial condition cases and x = 0.01 − 0.02 for the lock-release runs.The time-stepping is controlled internally to keep maximum Courant number less than 0.625.The DJL solution initial condition runs use no normal flow boundary conditions in z and inflow/outflow boundary conditions in x.A uniform flow u = −c with s = s(z) is imposed at the right boundary with an open boundary at the left.The waves are traveling in the positive-x direction in the laboratory frame.The lock-release runs have no normal flow conditions on all boundaries.

Fig. 2 .
Fig. 2. Evolution of an extended-s DJL solution for λ=8 and c = 0.4.Contours of s in intervals of 0.1 for s≤1 (outside the core) and 0.025 for s>1 (inside the core) are plotted.Times are indicated in the panels.Note that only part of the domain is shown.

Fig. 3 .
Fig. 3. Evolution of an extended-s DJL solution for λ = 8 and c = 0.58.Contours of s in intervals of 0.1 are plotted.

Fig. 4 .
Fig. 4. The region in the core of the extended-s DJL solutions in (a) Fig. 2 and (b) Fig. 3 where 0≤Ri≤0.25.The heavy solid line shows the core boundary (s = 1) and the thin lines show Ri=0 (lower line) and Ri=0.25 (upper line) contours.

Fig. 5 .
Fig. 5.The vertical profiles of (a) u − c, (b) s, and (c) Ri through the center of the solitary waves at t = 0 (solid) and t = 200 (dashed) in Fig. 3.

Fig. 6 .
Fig. 6.Scatterplot (black) of the core vorticity versus streamfunction at t = 200 for the case in Fig. 3.The light gray region indicates the values of streamfunction and vorticity in the core at t = 0.

Fig. 7 .
Fig. 7. Stagnant-core DJL for λ = 8, s c and c=0.4.The core is indicated by the shading and the thin lines are s contours in increments of 0.1.The area inside the heavy line and above the core has Ri≤0.25.Note that only the lower half of the water column is shown.
Figure summarizes the properties the stagnant-core solutions.The amplitude, η M core height c and length L are shown as of c λ = and 8.For the extended-s wave are plotted.For values of λ, stagnant-core solutions gin at incipient overturning wave c * end reaching limiting conjugate solutions (open circles) c = and 0.468 λ = 4 8, respectively.was possible continue the beyond what is in the The reason the failure of numerical soluprocedure unclear, the likely is the emergence the regions low Ri.numerical difficulties been encountered Lamb when for solitary in nearly two-layered using The wave η M core h , and length L c wave speed from the DJL solutions are by the lines s =1 and = 8, (a)-(c), λ = (d)-(e).The circles indicate the state solutions.dashed and squares show same quantities from extendeds solutions.Nonlin.Geophys., 17, 303-318, www.nonlin-processes-geophys.net/17/303/2010/ K. R. Helfrich and B. White: Trapped-core solitary waves

Fig. 9 .
Fig. 9. Evolution of a stagnant-core DJL solution for λ = 8, s c = 1, and c = 0.4.The panels show contours of s at the times indicated.The contour intervals are 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that only the lower half of the domain is shown.

Fig. 11 .
Fig. 11.Evolution of a stagnant-core DJL solution for λ = 8, s c = 1.1, and c = 0.4.(a)-(d) Contours of s are shown in at the times indicated.The contour intervals are 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that only the lower half of the domain is shown.(e) Comparison of u−c and s at the wave crest at t = 100 (dashed) with the stagnant-core DJL solution with λ = 8, s c = 1 and c = 0.394 (solid).

Fig. 12 .
Fig. 12.Comparison of wave amplitude η M , core height h c , and core length L c versus wave speed c from the stagnant-core DJL solutions (solid lines), numerical solutions initiated with the model solutions (open circles), and from the lock-release runs (open squares) for λ = 8, (a)-(c), and λ = 4, (d)-(f).The solid circles show the theoretical conjugate state limits.The open diamonds in (a)-(c) show the equilibrated wave initiated with a stagnant-core DJL solution with s c = 1.1 and the open triangles show the conjugate state properties from a lock-release numerical run from White and Helfrich (2008).

Fig. 13 .
Fig. 13.Lock-release run for λ = 8, h d = 1, L d = 0.6 and s d = 1.Contours of s are shown at the indicated times.The contour interval is 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.Note that part of vertical domain is shown for t > 0.

Fig. 14 .
Fig. 14.The leading solitary wave from Fig. 13 at t = 300.(a) s, (b) u − c (c) w and (d) ψ.The wave travels at a constant speed c = 0.417 for t > 100.The s contour interval is 0.1 for s ≤ 0.9 and 0.025 for s > 0.9.The contour interval is 0.05 in (b) and 0.02 in (c), with solid lines for values ≥ 0. The ψ contour intervals in (d) are 0.001 for ψ ≥ 0.01 and 0.04 otherwise.(e) Vertical profiles of u−c and s at x − ct = 0 from the theory (solid) and numerical model (dashed).Note that the horizontal and vertical scales in (d) differ from those in (a)-(c).

Fig. 15 .
Fig. 15.The vorticity structure of the leading wave at t = 10 from the lock-release numerical run in Fig. .

Fig. 16 .
Fig. 16.Scatterplot of the core vorticity versus streamfunction for the wave in Fig. 14.