Nonlinear Processes in Geophysics Strongly nonlinear , simple internal waves in continuously-stratified , shallow fluids

Strongly nonlinear internal waves in a layer with arbitrary stratification are considered in the hydrostatic approximation. It is shown that “simple waves” having a variable vertical structure can emerge from a wide class of initial conditions. The equations describing such waves have been obtained using the isopycnal coordinate as a variable. Emergence of simple waves from an initial Gaussian impulse is numerically investigated for different density profiles, from twoand three-layer structure to the continuous one. Besides the first mode, examples of secondand third-mode simple waves are given.


Introduction
It is well known that in strongly nonlinear, non-dispersive media, such as those considered in water waves, gas dynamics, and magnetic hydrodynamics, a smooth initial perturbation, however strong it would be, disintegrates into a set of "simple" waves, each evolutioning in such a way that each point of the wave profile propagates at a constant velocity depending on the local perturbation.This speed is in general a smooth function of the local perturbation.As a result such a wave either stretches in space or steepens; in the latter case its profile eventually overlaps, preventing the further use of non-dispersive approximation and leading to the formation of a shock wave or, in a dispersive medium (such as the plasma) to oscillations and possible generation of a series of solitons.Here the existence of simple waves is discussed for strongly nonlinear waves with a vertical structure, in application to internal gravity waves in an incompressible stratified fluid, considered in the hydrostatic approximation.An important and interesting complication in our case, in contrast Correspondence to: L. A. Ostrovsky (lev.a.ostrovsky@noaa.gov)with the well-studied water waves and classical gas dynamics problems, is that in a stratified system the wave has a vertical dimension, and the local speed of a simple wave must be independent of the vertical coordinate.
The question of existence of simple waves in strongly nonlinear, continuously stratified systems seems to be of significant interest both from the heuristic and practical viewpoints.Indeed, it is now a common knowledge that in many cases the tide-generated internal waves can be "genuinely" strongly nonlinear in the sense that isopycnal displacement in the wave can be comparable with and even exceed the characteristic vertical scale of the stratification (e.g., pycnocline depth).A bulk of observational data confirms this statement; a review of the problem for dispersive, solitary waves can be found in, e.g., Apel et al. (2007).However, except for a weakly nonlinear case considered, for example, by Yermakov and Pelinovsky (1975), theoretical consideration of such processes was limited to the case of a two-layer fluid where only one internal mode exists.Even in this case the main attention was concentrated on the formation and parameters of strong solitons.At the same time, as follows from the relevant models and observations both in deep and shallow ocean areas, in many cases strong internal waves are long enough to remain practically non-dispersive and deform only due to the dependence of the local wave velocity on the displacement, i.e., as simple waves.Even when the groups of solitons are observed, these groups are separated by stretching intervals of weaker, long waves.These cases were modeled much less thoroughly than the solitons, and the aim of the present paper is to address this gap.
In the non-dispersive, two-layer case the exact expression for the simple wave velocity c(η), where η is the local displacement of the interface (pycnocline), has been obtained (Sandström and Quon, 1993;Baines, 1995;Zahibo et al., 2007) and studied for general initial conditions by numerous authors.See, for example, Lyapidevsky (2000) for an application to roll waves.Recent work Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.L. A. Ostrovsky and K. R. Helfrich: Simple waves in continuously-stratified, shallow fluids by Chumakova et al. (2009) discussed the possibility of simple waves in layered and continuously stratified hydrostatic systems.However, they did not show that simple waves do indeed develop from more general initial conditions.Here we consider long-wave evolution in a general case of strongly nonlinear, long waves in a smoothly stratified fluid.One of the interesting questions to be answered is whether a progressive, distorting wave (similar to a simple or Riemann wave in gas dynamics) can emerge from an arbitrary localized initial perturbation.
The governing equations for fully-nonlinear, continuously stratified, hydrostatic flows are derived in Sect. 2. In Sect. 3 these equations are then used to obtain a set of nonlinear partial differential equations that govern the simple wave dynamics in this system.These equations can be reduced to the known results for weakly nonlinear waves with continuous stratification and fully-nonlinear waves in a two-layer system.However, analytical solutions for simple waves for general initial conditions and general stratification profiles do not seem possible at this point.In Sect. 4 initial value numerical solutions are used to demonstrate that simple wave behavior in systems with many layers (approaching a continuous stratification) does indeed exist and are consistent with the theory of Sect. 3. We note that in hyperbolic systems such as studied here, shocks are expected to form and our solutions do develop shocks.However, our focus is not on the shock dynamics, but rather on the rarefying portions of the solutions where simple wave behavior can be unambiguously identified.In Sect. 5 we give an example of a simple wave with complex vertical structure and conclude in Sect.6.

Basic equations
Hydrodynamic equations for a inviscid, non-diffusive, stratified fluid (for simplicity, a 2-D problem is considered; generalization to the 3-D case can readily be made) are Here u and w are horizontal and vertical components of velocity vector, p is pressure, ρ is total density, ρ 0 is a constant reference density, and g is gravity acceleration.These equations assume the Boussinesq approximation which is almost always acceptable in the case of the ocean where the density variations are small (ρ/ρ 0 ≈ 1).
For an arbitrary stratification, we apply the approach similar to that suggested in Ostrovsky (1978) for weakly nonlinear motions, where the Boussinesq equations for internal modes have been derived.Namely, the basic hydrodynamic equations will be written in the variables x, t, and β, where β is the vertical coordinate of a resting isopycnal.In these variables the density ρ depends only on β.
The local height z is denoted as h(β,x,t) which is the level of an isopycnal characterized by its initial level β, so that ζ = h − β is the isopycnal displacement in the wave, and at infinity where the perturbation is absent, h = β.At an isopycnal, vertical velocity is w = ζ t + uζ x where the derivatives are taken at a constant β.Also, for a function f (t,x,z) ⇐⇒ f (t,x,β) we have (here, by definition, h z = 1) In particular, Here ∂/∂z is the derivative at constant x as in the above system Eqs.(1) − (4).Also Here we distinguish between x-derivatives at constant z as above and at constant β.Similarly, Consider first the continuity Eq. ( 1) which now has the form In the new variables it reads From the dynamic Eqs.(3) and (4) it follows, upon differentiating them over x and z, respectively, and then subtracting the results gives As mentioned, we consider long, non-dispersive waves which means that w u, and ∂/∂x ∂/∂z.This can be formalized by introducing the corresponding spatial scales from which it readily follows that all terms in the right-hand side of Eq. (6) are small compared to each term on the left-hand side and can be neglected.Thus, or, in new variables, Here N 2 (β) = (g/ρ 0 )∂ρ/∂β is the given square of buoyancy frequency (same as its distribution N 2 (z) at infinity).We note that in the final Eqs. ( 5) and ( 7) the x-derivatives are at constant β, while in the original variables (e.g., ( 6)) the x-derivatives are at constant z.Equations (5) and (7) are the exact long-wave (nondispersive) equations which are constructive for analyzing the propagation and formation of simple waves.In what follows we consider a flat bottom at z = 0 and the rigid-lid condition on the surface z = H (which is typical of the consideration of internal waves in the Boussinesq limit): w = 0 at z = 0,H .

Simple waves
Assume now that a non-steady progressive wave does exist which at each level β propagates as a simple wave depending on two variables, z and ξ(x,t) = x −c(ξ )t.Let us look for the solution in the form h = h(ξ,β) and u = u(ξ,β), where ξ and β are considered independent.Thus, we have Substituting this into Eqs.(7) and (5) we obtain, respectively, h ξ ξ t β + ξ x h β u ξ + uh ξβ = 0.
Taking into account Eq. (8) and independence of ξ on β, we have Hence, the problem is reduced to a system of lower order as the x and t dependence is reduced to just ξ .
To relate these equations to the known results, we consider two particular cases.Let us first apply equations in the form Eq. (10) to small-amplitude waves.In the linear approximation, after letting c = c 0 and representing the solution in the form h = β + h (ξ,β), u = u(ξ,β), where h and u are small, (10) becomes After integrating these equations over ξ we have where F 1,2 are arbitrary functions.Seeking a solution in the form we have to let F 1,2 = 0, so that A(ξ ) is an arbitrary function and With the conditions f (0) = f (H ) = 0, this is a usual linear eigenvalue problem defining the vertical mode structure f (β) and long-wave propagation speed, c 0 .Now consider a weakly nonlinear wave.At small nonlinearity, letting c = c 0 + c (ξ ), we rewrite Eq. (10) in the form Here the nonlinear terms on the right-hand sides of this system are assumed small as compared with each linear term on the left-hand side.Now we will use the asymptotic perturbation method.Let h = h 1 + h 2 and u = u 1 + u 2 , where the subscripts 1 and 2 denote the first-and second-order perturbations, respectively.Then the second-order correction to Eq. (13) yields From here, after substituting u 1 = c 0 h 1β from Eq. ( 11) into the right-hand side of these equations and eliminating u 2 , we obtain with h 2ξ (0,H ) = 0.Here Eq. (11) has been used to eliminate h 1 .
According to the Fredholm theorem, for the existence of a finite solution for this inhomogeneous boundary problem for the function h 2ξ , the right-hand side of Eq. ( 14) must be orthogonal to the eigenfunctions of the left-hand-side linear operator, thus From this the known expression for c is found to be (e.g., Yermakov and Pelinovsky, 1975;Ostrovsky, 1978), This expression determines the nonlinear correction to wave velocity which is, as expected, proportional to the displacement A(ξ ).As expected, when N 2 is constant and f ∝ sinkβ from Eq. ( 12), c = 0, and nonlinearity reveals itself only at higher order.Examples of the application of weakly nonlinear theory to single and multi-modal cases can be found in Yermakov and Pelinovsky (1975 As a further test of the simple wave Eq. ( 10) it is shown that they give the known simple wave results (i.e., characteristic, or phase speeds and the corresponding Riemann invariants) for the fully-nonlinear two-layer system with a rigid lid (Sandström and Quon, 1993;Baines, 1995;Zahibo et al., 2007).In what follows the lower layer is labeled 1 and the upper layer 2, with a density jump of ρ 2 −ρ 1 at β = β 0 .In the first approximation the horizontal velocities u 1 (ξ ) and u 2 (ξ ) in each layer are independent of β.The vertical location h of the interface is approximated as a linear function of β, so that for each layer Here After that the second Eq. ( 10) yields From where it follows that In this two-layer case, N 2 = g δ(β −β 0 ), where g = (ρ 1 − ρ 2 )/ρ 0 .Thus, integrating the first Eq.( 10) over β in a small vicinity of β = β 0 , gives Substituting Eq. ( 18) into this equation gives from which the two-layer characteristic speeds, c ± , are Equations ( 18) and ( 19) can be manipulated to give . Integration then gives the Riemann invariants on the characteristics c ± from Eq. ( 19).Both Eqs. ( 19) and ( 20) are the same as found in Sandström and Quon (1993), Baines (1995), and Zahibo et al. (2007).

Evolution of strongly nonlinear perturbations
Deriving a general analytical solution to the continuouslystratified, simple-wave theory Eq. ( 10) does not appear possible at this point.However, we can obtain numerical solutions of the continuously-stratified, shallow-water Eqs. ( 5) and ( 7) and analyze the results for evidence of simple wave behavior.The integration is done by taking the fluid to be composed of M layers.The thickness of each layer is d j , j = 1,2,...,M (labeled from the bottom up) and the interface between layer j and j + 1 is at h j +1/2 .Far away from the disturbance the depth of interface h j +1/2 → β j +1/2 .Note that the fluid is contained between a flat bottom at h 1/2 = β 1/2 = 0 and a rigid lid at h M+1/2 = β M+1/2 = H .Each layer has uniform velocity u j and density ρ j .The reduced gravity between layers g j +1/2 = g(ρ j − ρ j +1 )/ρ 0 .With these definitions, where D j = β j +1/2 − β j −1/2 is the resting thickness of layer j .Similarly Here β j is the resting position of mid-point of layer j .
Thus the layered version of the continuity Eq. ( 5) is and Eq. ( 7) becomes an equation for the velocity jump u j +1 − u j between layers for j = 1,2,...,M − 1, and These give 2M − 1 equations in 2M unknowns, u j and d j , j = 1,2,...,M.The final equation comes from summing ( 21) over all the layers and noting that from the rigid-lid assumption.This gives where we have assumed that u j = 0 far from any disturbance.
We also eliminate solving one of the individual continuity Eq. ( 21) by using Eq. ( 23).
The layered Eqs. ( 21)-( 24) are solved with the nonoscillatory, shock-capturing method of Jiang and Tadmor (1998).A shock-capturing method was employed since these nonlinear, hydrostatic equations can be expected to lead to wave steepening and breaking.This numerical method is advantageous since it naturally allows the numerical solutions to proceed smoothly when shocks form as demonstrated in the results below.We note that the equations above, while in a flux form suitable for shock capturing, do not preserve momentum flux across a shock.Indeed, even the two-layered version of these equations do not possess this property.The fundamental issue concerns the distribution between the layers of the energy loss across the shock.A discussion of this issue can be found in Klemp et al. (1997) for two-layer flows.Jiang and Smith (2001) highlight the role of viscosity in resolving the problem.Since our focus is not on the shock dynamics, we allow shocks to form, but restrict our analysis of the numerical results to shock-free regions where we search for simple-wave behavior.Note that if either non-hydrostatic effects or significant (turbulent) viscosity were included in the model the shocks would have a finite length or not form at all.
In all the results reported the numerical solutions use a grid step of x = 0.02.The time step adjusts automatically to keep the Courant number max(|u|) t/ x ≈ 0.5.In some case t is fixed for convenience of analysis.The lateral boundary conditions are simply that ∂/∂x = 0.The continuous, background density profile (non-dimensional) which the M-layer discretization approximates is given by where with The normalization in (26) gives s(1) = 0 and s(0) = 1 so that the dimensional density varies from ρ 0 + ρ at β = 0 to ρ 0 at β = 1.The tanh density profile in Eq. ( 27) was chosen for simplicity and to provide a definite example.It allows both the location of the interface, given by z 0 , and the thickness, set by λ, to be varied.The M-layer discretization is done so that the density difference is approximately the same across all layers, although in some cases additional layers are added to keep individual layer thickness below a maximum.Figure 1 shows s(β) profiles with z 0 = 0.6 and λ = 5 and 15.Also shown by the solid circles are the densities at the layer mid-points for a discretization for M = 20.
Most of the time-dependent numerical solutions of Eqs. ( 21) -( 24) use a Gaussian-shaped initial condition for the interface displacements, with amplitude a 0 and width b.Here φ(β) is the vertical structure function for the linear vertical mode of interest (max(φ) = 1).The vertical mode structure is found numerically from the corresponding M-layer linear eigenvalue problem that also gives the linear long wave phase speed c 0 .
The initial condition for the velocity field is found from the linear approximation to Eq. ( 21) Before proceeding to the continuously stratified cases (i.e., large M) the numerical code was tested for a two-layer case (M = 2) where simple wave solutions can be found analytically from the characteristic speeds Eq. ( 19) and Riemann invariants Eq. ( 20). Figure 2 shows the numerical solution for the interface at the times indicated for a two-layer case with the resting layer depths D 1 = 0.6 and D 2 = 0.4 (z 0 = 0.6 and λ → ∞ in Eq. ( 27)) and an initial Gaussian disturbance Eq. ( 28) with a 0 = 0.25 and b = 1/3.The initial velocity field is a simple wave propagating to the right with R − = sin −1 (D 2 − D 1 ) in Eq. ( 20).The leading face of the initial disturbance rarefies and a shock quickly forms on the trailing face.Both are expected from the simple-wave characteristic speed, c + , that decreases monotonically from c 0 = (D 1 D 2 ) 1/2 as the interface displacement ζ = d 1 − D 1 increases above zero (for this example with a 0 > 0).The dashed lines in the figure are the theoretical simple wave so- lution.Prior to the onset of breaking at t ≈ 4.76, the numerical and theoretical solutions are indistinguishable.After this time the agreement between the theoretical and numerical solutions is excellent on the leading face of the disturbance ahead of the shock.For simplicity, theoretical solution does not include the trailing shock that could be obtained using a shock-joining analysis.
The comparison in Fig. 2 indicates that the numerical solution procedure is very accurate, especially in the smooth parts of the flow that are of interest here.We do not show, but do note, that the two-layer model also accurately captures the initiation of multiple shocks due to the non-monotonic, simple-wave relation between c and displacement that can occur for certain layer depths and initial interfacial displacements (Smyth and Holloway, 1988;Zahibo et al., 2007).
In this two-layer example, and the numerical solutions below we will focus on stratifications with z 0 > 0.5 and initial conditions with a 0 > 0 and c 0 > 0. As seen above, this initial condition produces a leading rarefaction propagating in the positive-x direction.This isolates the simple wave behavior ahead of any trailing shocks and simplifies the analysis for simple waves.If a 0 < 0 and z 0 > 0.5 a shock would form on  the leading face and for general stratifications the shock will generate slower, higher-mode disturbances that would then interfere with the trailing rarefaction (see below).While the shocks are certainly interesting and merit study, our focus is on simple wave dynamics.
When the stratification consists of more than two layers there is no analytical simple wave solution.Figure 3 shows a numerical solution for a three-layer stratification with a thin middle layer that closely approximates the two-layer case in Fig. 2. In this example D 1 = 0.54, D 2 = 0.12 and D 3 = 0.34, with scaled density jumps 3/2 = 5/2 = 0.5.The initial condition is given by Eq. ( 28) and the linear wave velocity structure from Eq. ( 29) with a 0 = 0.25, b = 1/3 and c 0 = 0.458.This case also produces a leading rarefaction with a first vertical mode structure.At t = 6 the upper interface is approaching breaking (i.e., it is nearly vertical).As the breaking proceeds (t = 12 − 20) a second-mode disturbance develops behind the trailing shock of the leading disturbance.Interestingly, this second-mode wave develops shocks on both the leading and trailing ends.However, it is not clear that this behavior is physically correct as there are questions about the conservation properties of the governing equations across discontinuities.Small first and secondmode disturbances propagating to the left can also be seen.
If the leading rarefaction is evolving as a simple wave, the speed c of any point along either interface should be independent of the interface (i.e., β) and depend only on x.This was tested by taking the numerical solution at any   15) and the initial condition.(c) Close-up of the interface positions from the numerical solution in Fig. 6 at t = 60 (blue solid lines) and the prediction from assuming simple wave behavior from (30) with c(ξ ) from (b) (red dashed lines).t = t * and extracting the propagation speed of points on each interface.This is done by differencing the x-positions (at two closely separated times around t * ) of the interface point with displacement ζ j +1/2 (x,t * ).For example, at t * = 30, and x = 15 in Fig. 3, the upper interface has an amplitude ζ 5/2 (15,30) = 0.149.At t * − dt, the upper interfacial point with ζ 5/2 = 0.149 is at x = 15 − dx 1 and at t * + dt, it is at x = 15 + dx 2 .The speed of the interface point with amplitude ζ 5/2 = 0.149 is then c = (dx 2 + dx 1 )/(2dt).This procedure gives the speed c j +1/2 for each interface as a function of ξ (= x at t = t * ).
The computation of c j +1/2 (ξ ) for x = 15 − 34 at t * = 30 from Fig. 3 is shown in Fig. 4a.The speeds along both interfaces are shown, the lower is dashed and the upper is solid, and they are almost exactly the same and thus not distinguishable on the figure as it is plotted.The largest difference is 0.001 at x = 15 and is < 10 −5 for x > 18.The solution is evolving as a simple wave since the speeds are independent of the interface.A further test of simple-wave behavior is shown in Fig. 4b where a close-up of the numerical solution at t = 60 is plotted.Also plotted is the estimate of the interface positions (dashed lines) found by assuming simple wave evolution of the t = 30 solution from Here c(ξ ) is found at t * = 30 in Fig. 4a.The interface positions at t = 30 and 15 < x < 34 propagate with speeds c(ξ ) independent of the vertical position (i.e., β or j ).With the exception of the shock, that propagates faster than trailing portion of the rarefaction, the agreement is excellent.
Another test for the presence of a simple wave is to solve the simple wave Eqs.(10) for the c(ξ ) derived above.In the layered system, the non-dimensional versions of Eq. ( 10) are and If d j and u j are known at one end of the domain where c(ξ ) is given, say ξ 0 , then Eqs. ( 31) and (32), along with the constraints Eqs. ( 23) and ( 24), can be integrated to obtain d j (ξ ) and u j (ξ ).
Figure 5 shows a comparison of the full numerical solution in Fig. 3 at t = 30 (solid lines) and the simple wave solution (dashed lines) found from integrating Eqs. ( 31) and (32) starting at ξ 0 = x = 15 using a Runge-Kutta method and c(ξ ) from Fig. 4a.The agreement is very good, again indicating that the leading disturbance is evolving as a simple wave.
The role of continuous stratification is explored for the density profiles from Eq. ( 27) shown in Fig. 1.Both cases use M = 20 layers as shown in the figure.The time-dependent numerical solution for the stratification with z 0 = 0.6 and λ = 15, and c 0 = 0.363 is shown in Fig. 6.The initial condition is Gaussian-shaped from ( 28) and ( 29) with a 0 = 0.25 and b = 1/3 as in Figs. 2 and 3.The evolution is qualitatively the same as above with a leading rarefaction followed by a shock.Slower second and higher vertical mode disturbance are also present.The very high wavenumber oscillations (at near the grid spacing) behind the shock at t = 15 and later appear to be the expression of Kelvin-Helmholtz instability in analogy with the ill-posed nature of the two-layer hydrostatic system for large vertical shear.However, in this case the instability does not grow to overwhelm the solution.
Analysis of the propagation speeds of points along the interfaces of the leading rarefaction are shown in Figs.7a and  b. Figure 7a shows the speed c as a function of β for the x-locations indicated by the dashed vertical lines in Fig. 6 at t = 30.Just ahead of the shock at x = 14 , the speeds of the interfaces in the lower part of the water column are greater than the upper part, with a two-layer-like depen-dence of c on β.The speed c becomes essentially independent of β for x ≥ 16.The minimum, mean, and maximum values in β of c(x) for 14 ≥ x ≥ 34 are shown in Fig. 7b.Both figures testify to the formation of a simple wave.Also shown by the dash-dot line is c(x) from the weakly nonlinear approximation (15) with the same initial condition and stratification.Only the rarefying portion of the solution ahead of the crest is shown.Clearly, this numerical example demonstrates simple-wave behavior well beyond the weakly nonlinear regime.A comparison of the full numerical solution at t = 60 with the prediction from assuming simple wave evolution from (30) is given in Fig. 7c.This calculation used the mean c(ξ ) in Fig. 7b (t * = 30).With the exception of the shock, the agreement is excellent.
Unlike the three-layer example shown in Fig. 5, attempts to directly solve the layered simple wave Eqs. ( 31) and (32) using c(ξ ) from the time-dependent numerical solution were not successful.This appears to be a result of the sensitivity of the problem to small errors in c(ξ ).However, the full numerical solution at t = 30 was found to satisfy ( 31) and (32) in 14 ≥ x ≥ 34 with maximum residuals of < 10 −5 .
Figures 8 and 9 show the same plots for the second case, with a more diffuse interface, λ = 5 and c 0 = 0.339.The other parameters are all unchanged.The dependence of c on β is slightly more pronounced (Fig. 9a and b) but again c does become independent of β as the leading portion of the rarefaction is approached.The application of (30) in Fig. 9c indicates that the leading disturbance is evolving as a simple wave.The comparison with the weakly nonlinear theory in Fig. 9b again indicates that simple waves exist in the fully nonlinear regime.
The vertical structure of the interface displacements, ζ (β), from the example in Fig. 8 at t = 30 between x = 16 and 24 are shown in Fig. 10a. Figure 10b shows comparison of the vertical structure at x = 24, where the displacements are small, and the vertical structure from the linear eigenvalue problem scaled to have the same maximum displacement as the full numerical solution.A similar comparison between the linear vertical structure and the vertical structure in the large-displacement region (x = 16 − 20) shows significant differences.
Finally, the development of simple wave behavior in a mode-two disturbance is shown in Fig. 11 for a three-layer case with D 1 = D 2 = 0.4, and D 3 = 0.2, and 3/2 = 5/2 = 0.5.The initial condition is again from Eqs. ( 28) and ( 29) with a 0 = 0.23, b = 1/3 and c 0 = 0.235.The unbalanced initial condition produces both first-and second-mode disturbances that propagate to the right.By t ≈ 50 the faster firstmode disturbance (with the frontal shock) has moved ahead of the slower mode-two rarefaction (with the rear shock).The speeds of points on each interface at t = 75 in the modetwo rarefaction (15 ≤ x ≤ 32) are shown in Fig. 12a.The computed speeds are noisy, with the upper interface generally propagating slightly faster than the lower one.The noise probably reflects the difficulty of making finite difference www.nonlin-processes-geophys.net/18/91/2011/ Nonlin.Processes Geophys., 18, 91-102, 2011   estimates in this slower moving mode-two disturbance that may not be free of small amplitude mode one waves.The different speeds on each interface would presumably come together if the simple wave behavior had more time become established.Despite these issues, the comparison of the numerical solution at t = 100 with an application of the sim-ple wave assumption for interface evolution (30) shown in Fig. 12c is quite good.The simple wave estimate was found using t * = 75 and c(ξ ) given by the mean speed on the two interfaces.28) and ( 29) with a 0 = 0.25 and b = 1/3 and the vertical structure from the second vertical mode.

Construction of a prescribed simple wave
Although the solutions of the Eqs.(5) and (7) at the corresponding stage do satisfy the lower-order simple wave equations (10), it is often difficult to solve the latter directly since the velocity c(ξ ) is not known in advance, and the vertical (mode) structure is not fixed.However, Eq. (10) can be effectively used for constructing a wide spectrum of simple waves with any prescribed velocity c(ξ ) and different initial conditions.Sometimes the corresponding wave structure can be rather complex.Here we give an example of a solution with the following parameters: h(0,β) = β − gsin(3πβ), u(0,β) = −u 0 cos(3πβ), h(ξ,0) = 0, h(ξ,1) = 1, u(ξ,0) = −vcosξ, u(ξ,1) = vcosξ.
In the linear approximation this corresponds to a third-mode vertical structure which is harmonic in time.
Then equations (10) for the simple wave were solved with the initial conditions (33) using Mathematica.As an example, we have taken the following parameters: q = 0.2,g = u 0 = v = 0.1.Figures 13 and 14   time of breaking the hydrostatic approximation becomes inapplicable, and either dispersion (typically resulting in soliton formation) or dissipation (forming a shock wave -internal bore), or both, have to be taken into account.components are present, the vertical profiles at even and odd phases are different: the latter is significantly distorted with respect to the initial profile.Nevertheless, this solution is a simple wave as it is discussed above.

Conclusions
At present there exist many recorded observations of strongly nonlinear internal waves, mostly of those generated by tides in coastal zones.These waves often exist as solitary waves and their groups ("solibores"), presumably formed as a result of nonlinear steepening of long baroclinic waves generated by interaction of barotropic tidal waves with bottom features such as shelf breaks.However, the long-wave stage of the process has not been studied in any detail.Besides, the rarefaction stage may be of interest itself, since it may contain most of the wave energy and defines the total length of the internal tidal cycle.The models analyzed above allow a simplification of the description of strongly nonlinear quasihydrostatic processes while remaining exact and, in particular, expand the concept of simple waves known in many areas of fluid dynamics, to the processes with a complex (and not fixed) vertical structures.
Fig. 1. s(β) profile from Eq. (27) for z 0 = 0.6 and λ = 5 and 15.The solid circles indicate the values of s at the layer mid-points for M = 20.

Fig. 2 .
Fig.2.The solid lines show the interface positions from a two-layer numerical solution for D 1 = 0.6 and D 2 = 0.4.The initial condition is the Gaussian interface displacement Eq. (28) with a 0 = 0.25, and b = 1/3.The initial velocity field is for simple wave with R − = sin −1 (D 2 − D 1 ) in Eq. (20).The dashed lines shown for t ≥ 6 are the theoretical simple wave solutions that have become multi-valued after the onset of breaking at t ≈ 4.76.The numerical solution and the simple wave solution are nearly indistinguishable on the leading face of the wave.Note that only the upper half of the domain is shown.
Fig. 4. (a) Speed of propagation of points on each interface, c, at t = 30 from Fig. 3. Speeds of the lower and upper interfaces are indicated by solid and dashed lines, respectively.At the magnification of the plot these two curves are indistinguishable.(b) Close-up of the interface positions (solid lines) from the numerical solution in Fig. 3 at t = 60.Also shown (dashed lines) are the interface positions assuming simple wave behavior from Eq. (30) with c(ξ ) from (a).The numerical and simple wave lines are nearly indistinguishable on the the leading part of the wave (x > 25).

Fig. 5 .
Fig. 5. Interface positions at t = 30 from the numerical solution in Fig. 3 (blue lines) and the solution of the layered simple wave Eqs.(31) and (32) (red dashed lines).
L. A. Ostrovsky and K. R. Helfrich: Simple waves in continuously-stratified, shallow fluids

Fig. 7 .
Fig. 7. (a) Speed c of interfacial points as a function of β at x = 14, 16, 20, 22, and 24 computed from the numerical solution in Fig. 6 at t = 30.(b) Interfacial speed c as a function of x in the rarefaction region of the numerical solution at t = 30.The solid line shows the average and the dashed lines the maximum and minimum in β.The dash-dot line shows c(x) predicted from the weakly non-linear model Eq.(15) and the initial condition.(c) Close-up of the interface positions from the numerical solution in Fig. 6 at t = 60 (blue solid lines) and the prediction from assuming simple wave behavior from (30) with c(ξ ) from (b) (red dashed lines).

Fig. 9 .
Fig. 9. (a) Speed c of interfacial points as a function of β at x = 15, 17, 19, 21, 23, and 25 computed from the numerical solution in Fig. 8 at t = 30.(b) Interfacial speed c as a function of x in the rarefaction region of the numerical solution at t = 30.The solid line shows the average and the dashed lines the maximum and minimum in β.The dash-dot line shows c(x) predicted from the weakly non-linear model (15) and the initial condition.(c) Close-up of the interface positions from the numerical solution in Fig. 8 at t = 60 (blue solid lines) and the prediction from assuming simple wave behavior from (30) with c(ξ ) from (b) (red dashed lines).
Fig. 10.(a) The vertical structure of the interface displacements, ζ (β), in the leading rarefaction at the indicated x locations from the numerical solution at t = 30 in Fig. 8.(b) ζ (β) at x = 24 (solid with dots) and the linear eigenfunction scaled to have the same maximum displacement (dashed line).
show the resulting wave evolution.The wave front becomes vertical (breaks) at t ≈ 17, then it becomes multi-valued as shown for t = 20.Past the

Fig. 12 .
Fig. 12.(a) Speed of propagation of points on each interface, c, at t = 75 from Fig. 11.The lower (upper) interface is indicated by the solid (dashed) line.(b) Close-up of the interface positions from the numerical solution in Fig. 3 at t = 100 (blue solid lines) and the prediction from assuming simple wave behavior from Eq. (30) with c(ξ ) from (a) (red dashed lines).

Fig. 14 .
Fig. 14.Vertical structure of the isopycnal displacement h − β as a function of β at three values of ξ indicated next to each curve.The solution is 2π periodic in ξ .