Shallow water models for stratified equatorial flows

Our aim is to study the effect of a continuous prescribed density variation on the propagation of ocean waves. More precisely, we derive KdV-type shallow water model equations for unidirectional flows along the Equator from the full governing equations by taking into account a prescribed but arbitrary depth-dependent density distribution. In contrast to the case of constant density, we obtain for each fixed water depth a different model equation for the horizontal component of the velocity field. We derive explicit formulas for traveling wave solutions of these model equations and perform a detailed analysis of the effect of a given density distribution on the depth-structure of the corresponding traveling waves.


Introduction
The dynamics of inviscid fluids are governed by Euler's equations of motion. In the context of water wave models one usually assumes the density function to be constant throughout the entire fluid body, which is a reasonable simplification since the compressibility of water is negligible: even very high amounts of hydrostatic or hydrodynamic pressure have almost no effect on the water density. In fact, the pressure acts like a Lagrange multiplier, enforcing the incompressibility condition (see the discussion in [3]). However, other factors like temperature and salinity have a considerably higher effect on the density and give rise to stratified flows which become relevant e.g. in the context of ocean dynamics. A significant increase of density due to changes in temperature or salinity gives rise to a pycnocline-a layer that separates surface water of lower density from water of higher density in greater depths. The propagation of disturbances of a pycnocline gives rise to internal waves, which are much larger than surface waves but hard to detect since there is only weak interaction between internal waves and surface waves, cf. [19,11] and the references therein.
In recent studies of stratified equatorial ocean flows, cf. [4,11,12], the pycnocline has been considered as an infinitely thin interface between two layers of different but constant densities. This approach opened up rich developments in the theoretical investigations of equatorial flows, including the derivation of exact solutions [5,6,7], a generalization to three-dimensional fluid flows [14], discovery of the Hamiltonian structure of the model [8,9,22,23] and qualitative analysis of the thermocline see e.g. [32]; see also the overview articles [26,27] and the reference therein. A comprehensive account of phenomena pertaining to physical oceanography in general and to equatorial ocean dynamics in particular are presented in [13].
While the above approach assumes a predefined discontinuous piecewise constant density function with a jump at the pycnocline, the aim of the present paper is to study the effect of a continuous change of density on the propagation of water waves. To this end we derive the weakly nonlinear KdV-type shallow water equation for prescribed but arbitrary depth-dependent density distributions. Here the unknown u, depending on time τ and space ξ, describes the horizontal velocity component of the flow field beneath a unidirectional free surface shallow water wave of small amplitude at a specific water depth; the coefficient ν of the nonlinear term uu ξ , to be specified in the sequel, depends on the prescribed density distribution and the water depth. The major novelty compared to the case of constant density is the depth dependence of the horizontal velocity component of the flow at leading order. In particular, the vorticity of these flows does not vanish, cf. Remark 3.3. This is expected in view of the general principle that stratified flows can never be irrotational, see the discussion in [3] and [33]. The constant Ω 0 incorporates the Coriolis effect along the Equator. We obtain a similar equation for the elevation of the free surface. By a thorough study of the explicit traveling wave solutions of these KdV type equations we acquire a solid theoretical understanding of the effect of the density distribution on the size and shape of these traveling waves at various depths. This analysis shall serve as a basis for a better understanding on how a continuous change of density affects the dynamics of more general types of waves.

Governing equations in dimensionless form
In this section we introduce the governing equations for unidirectional gravity water waves; we provide the equations in physical variables and bring them into dimensionless form in the sequel. The main difference compared to the classical water waves problem is that we allow for a non-constant "background" density distribution depending smoothly on the vertical direction caused by e.g. increasing salinity or decreasing temperature with depth. Such stratification phenomena occur in particular in equatorial regions, where waves typically propagate in purely azimuthal direction, so that a unidirectional setting is a reasonable simplification, cf. [31] and the discussion in [11]. Motivated by this geophysical point of view, in which the effects of the earth's rotation on large waves become relevant, we consider the f -plane approximation for the inviscid Euler equations for two-dimensional unidirectional flows near the Equator. The motion of the water between the flat bed and the free surface is therefore governed by the following set of equations: where the bars distinguish these physical variables from dimensionless variables appearing later. The independent variablesx andz denote the directions of increasing azimuth and vertical elevation respectively; the horizontal and vertical fluid velocity component in the direction of increasing azimuth and elevation are given byū =ū(x,z,t) andw =w(x,z,t). The pressure is denoted byP =P (x,z,t). We consider a prescribed but arbitrary continuously differentiable density distributionρ of the water which may vary with depth only, that isρ The gravitational acceleration is denoted byḡ ≈ 9.81m s −2 and the constantΩ ≈ 7.29 · 10 −5 rad s −1 Figure 1. Fig. 1a illustrates the fluid domain in the physical (x,z)-plane between the flat bed atz = 0 and the free surfacez =h 0 +η(·,t) at a certain instant of timē t. The average water levelh 0 is indicated by a dashed line,λ shows the distance between two consecutive crests andā is the vertical deviation of a typical crest from h 0 . Fig. 1b shows a prescribed depth dependent density distributionρ(z) with a significant density increase in the region between the two dotted horizontal lines close to the surface giving rise to a pycnocline.
denotes the rotational speed of the Earth around the polar axis. Thus, the twoΩ-terms in (1) capture the effects of the so-called Coriolis force. Let us emphasize at this point that all considerations we make in this paper equally apply to the corresponding non-geophysical settings, i.e. to unidirectional inviscid gravity water waves with a prescribed smooth depth-dependent density distribution. In this case we simply setΩ = 0. As already mentioned, we assume that the fluid body is bounded from below by a flat bed atz = 0; it is bounded from above by the free water surface which is assumed to be described by the graph of the functionh 0 +η(·,t) at any instant of timet, whereη(x,t) measures the deviation of the free surface from the average water depthh 0 at (x,t). We impose the usual dynamic boundary condition for the pressure at the water surface and the kinematic boundary condition for the vertical velocity at the surface and the flat bed:P the constantP atm denotes the atmospheric pressure.
To bring the set of equations (1) and boundary conditions (3) into dimensionless form, we use the following standard reference length scales: a typical amplitude of the surface waveā and the average undisturbed water depthh 0 as the vertical scales, and a typical wave lengthλ as the horizontal scale, cf. [10], [24] and [25] for more details. We introduce the following set of non-dimensional variables: where the dimensionless pressure P acts as a measure for the deviation ofP from the hydrostatic pressure distribution. The densityρ is scaled according tō By employing (4)-(5) and introducing the two fundamental dimensionless parameters referred to as amplitude and shallowness parameter, the set of equations (1) and boundary conditions (3) transform into the following dimensionless form: In order to transfer the free surface to a fixed boundary, we rewrite the boundary conditions at the free surface by means of Taylor expansions of the involved variables u, w and P about z = 1. This yields the expressions where we have used the fact that ρ(1) = 1 in view of the nondimensionalization (5). Finally, following [10], we perform the usual scaling on the system (7) and set Ω = εΩ 0 , where Ω 0 is an appropriate constant. The latter scaling is a valid manoeuvre since the dimensionless parameters Ω and ε are of the same order of magnitude when assuming typical values of ocean depth and wave amplitude for offshore waves near the Equator, see the discussion in [21]. Additionally we exploit that the shallowness parameter δ in system (7) can be scaled out in favor of ε via the transformation where the scaling of w is required to ensure conservation of mass in the resulting system. We note that this transformation remains non-singular for ε and δ tending to zero only if these parameters satisfy a certain asymptotic relation. That is, by imposing the transformation (11) one considers the restriction to shallow water waves of small amplitude: The application of (9), (10) and (11) on the system (7) yields the following system of governing equations for unidirectional equatorial waves in scaled and dimensionless form: 3. Derivation of a shallow water model The system (13) of governing equations for equatorial waves serves as a starting point for our derivations. We follow the techniques presented in [24,25] to obtain model equations for the surface elevation η and the horizontal velocity u. By retaining only the leading order terms of (13) (i.e. by setting ε = 0) we obtain the linear system from which we infer that the leading order approximation of η (again denoted by η) satisfies the linear wave equation: To see this, we first note that the second equation in (14) and the dynamic boundary condition imply that η = ρP for all z ∈ [0, 1]. In view of the first equation this yields that ρu t = −η x . Next we invoke the equation of mass conservation which is equivalent to (ρw) z = −ρu x , hence Differentiating both sides with respect to t, using the boundary condition at the top and the fact that ρ(1) = 1, we obtain that η tt = − 1 0 ρu xt ds on z = 1. Hence, we infer that η satisfies (15).
The general solution of (15) is of the form η(x, t) = F (x − t) + F (x + t) for arbitrary functions F and F . This motivates the introduction of the far-field variables to follow right-propagating waves governed by (13). Rewriting (13) in terms of the far-field variables (16) yields the system To obtain an asymptotic solution of system (17), we formally expand the respective variables η, u, w and p in the form q ∼ ∞ n=0 q n ε n . At leading order we obtain the linear system Integrating the second equation and using the boundary condition on the pressure we obtain that η 0 = ρP 0 for all z ∈ [0, 1]. In view of the first equation we may therefore deduce that where we have invoked the assumption that any change in u can occur only due to a perturbation of η, see [24]. From the equation of mass conservation, the bottom boundary condition and (19) we deduce that w 0 = −zρ −1 η 0ξ . Next we consider the first order system where we have written all leading order expressions in terms of η 0 . Viewing the second equation as a linear ODE in P 1 with z being the independent variable, we can solve it using the pressure boundary condition and obtain that This results in the following representation of P 1 in terms of η 0 and η 1 : Next we solve the equation of mass conservation as a linear ODE in w 1 . By expressing the inhomogeneity using the first equation, and incorporating the boundary condition on z = 0 and P 1ξ obtained in the previous step, we find that Finally, we calculate the integral in (21), evaluate the resulting expression for w 1 at z = 1 and compare the outcome with the corresponding boundary condition in (20). This yields the weakly nonlinear third order partial differential equation where which governs the leading order approximation of the free surface η, where η ∼ η 0 + O(ε). The equation for the leading-order approximation of the horizontal velocity u at any depth z ∈ [0, 1] has a depth-dependent coefficient of the nonlinear term: Remark 3.1. For constant density ρ ≡ 1 and vanishing Coriolis force, both equations (22) and (24) reduce to the standard KdV [29] 2u τ + 3uu ξ + 1 3 which approximates right-propagating shallow water waves of small amplitude at leading order. We observe that for fixed z ∈ [0, 1] and for any given density function ρ the coefficient of the nonlinear term in both equations is a fixed real number. Therefore, these equations can always be brought to the form (25) via the transformation where α = ν(ρ) and u = η for (22) and α = ρ(z)ν(ρ) for (24). This transformation is valid only as long as α > 0. The sign of α can be determined for any arbitrary given density distribution ρ; we note that α is strictly positive for realistic density distributions in the context of (geophysical) water waves and provide examples for different densities at the end of Section 4.  (25), like soliton interaction, bi-Hamiltonian structure, the existence of a Lax pair, integrability, global Birkhoff coordinates, see e.g. [20,30,18,28], and stability of traveling wave solutions (cf. Section 4.3 for some references regarding stability) remain valid for equations (22) and (24) by means of the transformation (26).
Remark 3.3. The depth dependent coefficient in the velocity equation (24) due to non-constant density ρ represents a major difference compared to the standard KdV, which models both the free surface η and the horizontal velocity u at any depth. In the next section we will consider explicit traveling wave solutions of (22) and (24) for certain prescribed density distributions to illustrate the change of the corresponding wave profiles with depth. Moreover, a non-constant density distribution ρ entails a non-vanishing vorticity 1 of the flow field at leading order, cf. (18) and (19). This is to be expected in view of the rotational nature of stratified flows, see [3,33].

Explicit traveling wave solutions
In this section we present explicit traveling wave solutions of equations (22) and (24). We study the effects of the non-constant density on the z-structure of solutions in general and for two concrete examples of density distributions. Traveling wave solutions of (22) of the form η(ξ, τ ) = ϕ(ξ − cτ ) with wave speed c > 0 satisfy the ordinary differential equation with ν(ρ) defined in (23).
Remark 4.1. We observe that ρ appears only as a multiplicative factor in front of the solution and hence any changes in ρ affect only the amplitude, cf. Fig. 2a. This is different from the dependence on the Coriolis parameter and the wave speed: as c or Ω 0 increase, the profiles become taller as well as narrower, see [21]. Moreover, observe that we have made no assumption on the sign of ν(ρ). Therefore, in principle our approach allows us to obtain solitary waves of elevation when this number is positive and waves of depression when it is negative, see the remark at the end of Example 2.
For the horizontal velocity component u traveling wave solutions of the equation (24) of the form u(ξ, z, τ ) = φ(ξ − cτ, z) can be obtained directly from the explicit traveling wave solutions of (22) in view of the relation between the surface and the velocity component given in (19). This relation provides an explicit measure for the attenuation of the amplitude with depth which is inversely proportional to ρ(z), see also Fig. 2b. Therefore, we obtain a traveling wave solution of (24) for each depth z ∈ [0, 1] given by where ϕ is a solution of (27). We observe that the amplitude of the wave does not vanish at the bottom.
Example 1. Consider the linear density function where A > 0, i.e. ρ(z) is increasing with depth. In this case the coefficient of the nonlinear term in the equation (22) is ν(ρ) = 2 ln(1 + A)/A + 1 ∈ (1, 3), which is a strictly monotonically decreasing function of A. Therefore, for every fixed z ∈ [0, 1], the profiles of the solitary traveling wave solutions become taller and also wider as A increases, cf. (29) and Fig. 2a. The amplitude of the velocity component φ(ξ − cτ, z) decreases with depth like (Az) −1 in view of (19) and (31), see Fig. 2b.
Example 2. Next we consider the density function for suitable parameters a i ∈ R, to obtain a more realistic density distribution as described in the introduction and depicted in Fig. 1b. We choose the parameters to match the actual physical situation in the mid-equatorial Pacific, with an average water depth of about 4000 m and a pycnocline located at a depth of roughly 200 m. In this oceanic region, the density of the water at the surface is about 1027 kg/m 3 with a density increase of 1% from the surface to the seabed, cf. [11]. Nondimensionalizing and scaling these quantities to depths z ∈ [0, 1] and ρ(1) = 1 leads us to the following choice of parameter values: a 3 = 0.95 determines where the region of biggest density gradient is located, i.e. at the center of the pycnocline; a 2 = 300 regulates the thickness of the pycnocline; with the last two parameters a 0 = 1.0047 and a 1 = 0.0034 we ensure the correct density increase. This choice results in the density distribution depicted in Fig. 3a.  (22) with linear density function ρ(z) = 1 + A(1 − z). In Fig. 2a we see that for larger values of the parameter A > 0 the profile becomes taller and also wider. In Fig. 2b we see how the amplitude of solutions decreases with depth z ∈ [0, 1] like (Az) −1 according to (30).
It is interesting to see how this density profile is reflected in the decay of the amplitude of the solitary wave profile (29) with depth. This is because the solitary waves φ(ξ − cτ, z) of the velocity equation (24) decay with depth inversely proportional to ρ(z) according to relation (30), see Fig. 3b 2 .
Finally, we observe that for the choice of parameters described above the coefficient of the nonlinear term ν(ρ) = 2.97 defined in (23) is positive. We note that in this setting negative values of ν(ρ), which give rise to solitary waves of depression, would arise only for a density increase of almost 300% which is not a physically relevant scenario.

4.2.
Periodic traveling waves. Periodic traveling waves of (22) can be obtained in a similar fashion, cf. [17] or [16]. To this end we integrate (27) twice to obtain (ϕ ) 2 = F (ϕ) where F (ϕ) := −ν(ρ)ϕ 3 + 6(c + Ω 0 )ϕ 2 + Aϕ + B, where A, B are constants of integration. One can show that this equation has real and bounded solutions only if F (ϕ) has three simple roots satisfying ϕ 3 < ϕ 2 < 0 < ϕ 1 when ν(ρ) > 0 (the other situation being similar), in which case we can write The solution oscillates between ϕ 2 ≤ ϕ ≤ ϕ 1 with period 2 . Therefore we can express the solution implicitly as  Figure 3. In Fig. 3a we see a plot of the density profile ρ(z) = a 0 − a 1 arctan (a 2 (z − a 3 )) defined in (32) for suitable choices of parameter values a i ∈ R to model a density increase of 1% from surface to bed. Fig. 3b shows a schematic representation of the fact that the amplitude of solitary wave solutions φ(ξ − cτ, z) decays with depth inversely proportional to ρ(z), cf. (30).

Stability.
Having established existence and some qualitative properties of traveling waves, we may ask about the stability of these solutions. The orbital stability of periodic and solitary traveling waves of the standard KdV have been studied by many authors, see for instance [1,2,15]. Recall the observation presented in Remark 3.1 that for every fixed depth z ∈ [0, 1] the KdV type equations (22) and (24) derived in Section 3 can be transformed into the standard KdV. For this reason, we infer that the explicit solitary and periodic traveling wave solutions (29) and (30) of equations (22) and (24), respectively, are orbitally stable as well.