Articles | Volume 29, issue 2
Nonlin. Processes Geophys., 29, 207–218, 2022
Nonlin. Processes Geophys., 29, 207–218, 2022
Research article
 | Highlight paper
15 Jun 2022
Research article  | Highlight paper | 15 Jun 2022

Effects of rotation and topography on internal solitary waves governed by the rotating Gardner equation

Effects of rotation and topography on internal solitary waves governed by the rotating Gardner equation
Karl R. Helfrich1, and Lev Ostrovsky2,3, Karl R. Helfrich and Lev Ostrovsky
  • 1Department of Physical Oceanography, Woods Hole Oceanographic Institution, Woods Hole, MA, USA
  • 2Department of Applied Mathematics, University of Colorado, Boulder, CO, USA
  • 3Department of Mathematics, University of North Carolina, Chapel Hill, USA
  • These authors contributed equally to this work.

Correspondence: Lev Ostrovsky (


Nonlinear oceanic internal solitary waves are considered under the influence of the combined effects of saturating nonlinearity, Earth's rotation, and horizontal depth inhomogeneity. Here the basic model is the extended Korteweg–de Vries equation that includes both quadratic and cubic nonlinearity (the Gardner equation) with additional terms incorporating slowly varying depth and weak rotation. The complicated interplay between these different factors is explored using an approximate adiabatic approach and then through numerical solutions of the governing variable depth, i.e., the rotating Gardner model. These results are also compared to analysis in the Korteweg–de Vries limit to highlight the effect of the cubic nonlinearity. The study explores several particular cases considered in the literature that included some of these factors to illustrate limitations. Solutions are made to illustrate the relevance of this extended Gardner model for realistic oceanic conditions.

1 Introduction

Oceanic internal waves are an important class of nonlinear wave processes. In particular, the internal solitary waves (ISWs) are the most ubiquitous type of solitons in the natural environment. These waves often propagate for long distances over several inertial periods, and the effect of Earth's background rotation is potentially significant (e.g., Farmer et al.2009; Grimshaw et al.2014; Helfrich2007; Ostrovsky and Helfrich2019). The large ISWs in the South China Sea are prominent examples (Zhao and Alford2006; Alford et al.2010). There are also numerous remote sensing images throughout the coastal oceans that show multiple wave packets (e.g., Jackson2004), indicating that the ISWs persist over periods longer than the local inertial period. It is also known that rotation destroys internal solitons due to resonant radiation of inertia–gravity waves (terminal damping; see Grimshaw et al.1998a). The theoretical modeling of such processes often uses the rotation-modified Korteweg–de Vries equation (rKdV) derived for nonlinear waves in rotating ocean (Ostrovsky1978). In application to the oceanic conditions, this equation may need additional modifications, specifically for the variable depth along the propagation path. This specific case was considered in Grimshaw et al. (2014), Ostrovsky and Helfrich (2019), and Stepanyants (2019).

The rotation-induced solitary wave decay can be suppressed in certain ambient shear flows wherein the sign of the rotation coefficient is changed (Alias et al.2014). In those cases, the rKdV equation supports solitary wave solutions. This interesting situation is not considered in this paper.

Another important feature of oceanic solitons is that, in many cases, they are strongly nonlinear, so that the Korteweg–de Vries approximation (KdV) involving only the quadratic nonlinearity is inapplicable (e.g., Apel et al.2007; Helfrich and Grimshaw2008; Ostrovsky and Grue2003; Ostrovsky and Irisov2017). A better approximation of the real processes can be given by the Gardner equation that extends KdV by adding a cubic nonlinear amplitude term. It is significant when the quadratic and cubic nonlinear terms are of a comparable small amplitude. In that situation, higher-order corrections for dispersion and nonlinear dispersive terms are negligible. For example, this is the case in a two-layer fluid when the layer depths are nearly equal. Gardner solutions preserve the qualitative features of strongly nonlinear waves, such as the existence of a limiting soliton amplitude at which its length infinitely increases, whereas a KdV soliton amplitude can increase unlimitedly, with its width tending to zero. Moreover, in many cases, the Gardner equation gives a good quantitative approximation for strong solitons beyond the formal limits of its applicability (Stanton and Ostrovsky1998; Pelinovsky et al.2000; Grimshaw et al.2004) and thus has become the phenomenological model of choice. There are numerous studies using the Gardner equation in applications to oceanic waves, including its extension for the rotating ocean (e.g., Holloway et al.1999; Obregon et al.2018; Talipova et al.2015).

In this paper, we make the next step by adding a sloping bottom effect to the Gardner model with rotation. Correspondingly, the results of Grimshaw et al. (2014), Ostrovsky and Helfrich (2019), and Stepanyants (2019) can be significantly modified. In particular, we discuss an interplay between the effects of nonlinearity, rotation, and inhomogeneity. Some realistic estimates are also given.

The work by Karczewska and Rozmej (2020) considers the analogs of the KdV and Boussinesq equations for shallow water waves with different relative orders of the following small perturbing factors: nonlinearity, dispersion, and bottom slope. Then various terms, such as nonlinear dispersive terms, can prevail over other perturbations. Here we consider a certain (although rather general) physical problem when all included terms are of the same order and the classical description is sufficient.

2 Rotating Gardner equation

A standard model for the evolution of large-amplitude oceanic internal solitary waves is the rotating Gardner (rG), or extended KdV equation with rotation and variable depth h(x) (Holloway et al.1999), as follows:

(1) x η t + c + α η + ν η 2 η x + β 3 η x 3 + c 2 Q x Q η = γ η .

Here the wave amplitude function η(x,t) depends on the horizontal position x and time t. The linear long wave phase speed c is found from an eigenvalue problem for the structure function Φ(z) of a specified vertical mode. Both c and Φ(z) are slow functions of x. The x-dependent coefficients α, ν, β, Q, and γ of the quadratic nonlinear, cubic nonlinear, non-hydrostatic, inhomogeneous, and rotation terms, respectively, are found as integrals over the depth of Φ or Φ and the background stratification ρ(z) and current u(z). They can be found in numerous publications (e.g., Holloway et al.1999; Grimshaw et al.2004) and are summarized briefly in Appendix A. When considering a spatially inhomogeneous situation, it is advantageous to switch from the (x,t) system to the (s,x) system, where

(2) s = 0 x d x c - t .

If, additionally, ζ=Q1/2η is introduced, then Eq. (1) becomes the following:

(3) s ζ x + α c 2 Q 1 / 2 ζ ζ s + ν c 2 Q ζ 2 ζ s + β c 4 3 ζ s 3 = γ ζ .

We note that here ζ2=Qη2 is the wave action flux. Equation (3) can be shown to have two conserved quantities, as follows:

(4) M = ζ d s and E = ζ 2 d s ,

where the integrals are over the full s domain (infinite or periodic). These are, respectively, mass- and energy-related quantities, as discussed further in the next section.

Additionally, any initial condition to Eq. (3) with γ≠0 must satisfy the zero mass requirement, as follows (Ostrovsky1978):

(5) T / 2 - T / 2 ζ ( s , 0 ) d s = 0 ,

where s=-t at x=0, and T is the length of the time domain.

In the absence of rotation, γ=0, Eq. (3) reduces to the Gardner equation. This permits the solitary wave solution, as follows:

(6) ζ = a 1 + B cosh [ σ ( s - κ x ) ] ,

which is described by the parameter B. Here,

(7) a = Q 1 / 2 α ν ( B 2 - 1 ) , σ 2 = c 2 α a 6 β Q 1 / 2 , κ = β σ 2 c 4 .

The amplitude, A, in terms of η (=ζQ-1/2) is as follows:

(8) A = a Q - 1 / 2 1 + B = α ν ( B - 1 ) .

From the mass constraint (Eq. 5), a solitary wave initial condition requires the addition of a constant pedestal, as follows:


with ζ given by Eq. (6).

There are three families of steady solitary wave solutions given by Eqs. (6)–(8) (Grimshaw et al.1999). When ν<0, solitary wave solutions require 0<B<1 and have polarity αA>0. They approach the sech2 KdV solitary wave as B→1, and as B→0, the solution approaches the maximum amplitude, Amax=-α/ν, or flat-top wave. When ν>0, solitary wave solutions require B2>1, and there are two solution branches. For B>1, αA>0, and there is no limit on the wave amplitude. The sech2 KdV wave is recovered as B→1 from above, and for B≫1, the solution approaches the sech solution of the modified KdV equation (i.e., the Gardner equation with α=0). The third branch occurs for B<-1 with the wave polarity αA<0. In this case, solitary wave amplitude has a minimum amplitude Amin=-2α/ν obtained as B-1 from below. This limiting wave has an algebraic structure, and for B-1, the solution again approaches the sech wave of modified KdV equation. For ν>0, there are also localized pulsating traveling wave solutions (breathers) that have a total negative mass between zero and the mass of the limiting solitary wave at B=-1 (Pelinovsky and Grimshaw1997).

3 Adiabatic evolution of solitary waves

Assuming that ζ→0 for |s|, and integrating Eq. (3), gives the following:

(9) ζ x + α c 2 Q 1 / 2 ζ ζ s + ν c 2 Q ζ 2 ζ s + β c 4 3 ζ s 3 = - γ s ζ d s .

Multiplication by ζ and integration in s from −∞ to gives the following:

(10) d d x - ζ 2 d s = - 2 γ - s ζ d s d s = - γ - ζ d s 2 .

The solution of the rG equation (Eq. 9) implies that the right side of Eq. (10) is zero and E is a constant. However, progress can be made if we start with a solitary wave from Eq. (6), assume that the inhomogeneous and rotational effects are very weak, such that the evolving wave remains a solitary wave but with slowly varying amplitude, and take the limits of integration to contain only the evolving solitary wave. With the solitary wave solution (Eq. 6) written as follows:

(11) ζ = a F ( y ) , F = 1 1 + B cosh ( y ) , and y = σ s .

Equation (9) gives the following:

(12) d d x a 2 σ I 2 = - γ a 2 σ 2 I 1 2 ,


(13) I n = - F n d y , n = 1 , 2 .

The result, Eq. (12), is a statement for the variation in the energy,

(14) E w = a 2 σ I 2 ,

of the slowly evolving solitary wave. In the absence of rotation, Ew is conserved; however, the wave mass,

(15) M w = a σ I 1 ,

is not necessarily conserved. Changes in Mw are compensated by the formation of a trailing shelf. For example, in the typical case with ν<0 of a solitary wave approaching a point of polarity reversal, α=0, the wave mass increases in magnitude (Grimshaw et al.1998b). Thus, the trailing shelf must have the sign opposite to the wave polarity. When ν>0, the situation is more subtle, but again, any variations in Mw are compensated by a shelf (Grimshaw et al.1999; Nakoulima et al.2004). In a homogeneous rotating environment, both Ew and Mw decrease with x and are compensated by the trailing wave radiation (Grimshaw et al.1998a). With both inhomogeneity and rotation, the energy will decrease with x, but the wave mass may increase or decrease, depending on the interplay between these two effects.

In the KdV equation, ρ0Ew is the solitary wave energy, pu dzdt. Here p and u are the first-order, wave-induced pressure and horizontal velocity fields. With the addition of the cubic nonlinear term in the Gardner equation, Ew is not exactly the wave energy but is still a good measure of it. The actual wave mass is Q-1/2Mw, with Q defined as in Eq. (A3d).

3.1 Rotating KdV equation

For later reference, we first consider the adiabatic theory for the rotating KdV equation (ν=0) from Grimshaw et al. (2014). The solitary wave solution is found from Eqs. (6) and (8), with B=1 as follows:

(16) ζ = a F , F = sech 2 ( y ) , σ 2 = c 2 α a 12 β Q 1 / 2 , a = A Q 1 / 2 .

From Eq. (13), 1=2 and I2=4/3, and Eq. (12) gives the following:

(17) d d x 4 3 a 2 σ = - 4 γ a 2 σ 2 .

This can be written as follows:


When integrated, this gives the following:

(18) A A 0 = Q 0 Q 1 / 2 a a 0 = Q 0 Q 1 / 2 w 0 w 1 - ( w 0 a 0 ) - 1 / 2 0 x γ w 2 d x 2 .

The 0 subscript indicates the initial value at x=0.

In the absence of rotation, γ=0, the conservation of wave action gives the following:

(19) A A 0 = Q 0 Q 1 / 2 w 0 w = Q 0 2 β 0 α 0 c 0 2 α c 2 Q 2 β 1 / 3 .

For a homogeneous rotating environment, Eq. (18) gives the following:

(20) A A 0 = 1 - x X eO 2 , X eO = c γ α A 0 12 β 1 / 2 .

The KdV solitary wave is completely extinguished by the radiation of the inertia–gravity waves in the finite distance XeO.

3.2 Rotating Gardner equation

When ν≠0, the Gardner solitary wave solution (Eq. 11) and Eq. (13) give the following:



(23) T ( B ) = tan - 1 B - 1 B 2 - 1 = sgn ( B ) tan - 1 B - 1 B + 1 .

Substituting Eqs. (7), (21), and (22) into Eq. (12) gives the following:

(24) d d x 24 β α 2 Q 2 c 2 ν 3 1 / 2 ( B 2 - 1 ) 1 / 2 - 2 T ( B ) = - 96 γ β Q c 2 ν T 2 ( B ) .

This equation can be integrated to obtain B(x); hence A(x) is derived from Eq. (8). Since B→1 as ν→0, the equation remains regular for cases with ν(x) changing sign. Note also that it remains real for ν<0 and 0<B<1, since tan-1(iy)=itanh-1(y).

Figure 1XeG/XeO as a function of B0 for the three Gardner solitary wave regimes, with ν<0 for 0<B0<1 and ν>0 for B2>1.


Figure 2A/A0 versus x/XeO for homogeneous conditions. (a) ν<0 for B0=[10-1,10-2,10-4] (solid lines from left to right). The dashed line is the rKdV equation solution (Eq. 20). (b) ν>0 and B0=[1.1,5,10] (from left to right). (c) ν>0 and B0=[-2,-5,-10] (from left to right).


Figure 3Comparison of the rG adiabatic radiation decay theory for homogeneous conditions (black lines) and rG numerical solutions for (a) B0=0.3, (b) B0=10-4, (c) B0=5, and (d) B0=-4. In panels (a) and (b), red (blue) indicates the numerical solutions of the rG equation (Eq. 30) for ϵ=2.5×10-5 (2.5×10-3). In panels (c) and (d), red (blue) indicates ϵ=2.5×10-3 (0.25).


While it is not necessary to make the sign of ν explicit in integrating Eq. (24), for ν<0 (0<B<1), it can be written as follows:

(25) 1 - B 2 B d B d x = 1 2 d d x ln β Q 2 α 2 c 2 | ν | 3 2 tanh - 1 1 - B 1 + B - 1 - B 2 + 8 γ 6 β | ν | c 2 α 2 1 / 2 tanh - 1 1 - B 1 + B 2 .

And for ν>0 (B2>1), it can be written as follows:

(26) B 2 - 1 B d B d x = - 1 2 d d x ln β Q 2 α 2 c 2 ν 3 B 2 - 1 - 2 tan - 1 B - 1 B 2 - 1 - 8 γ 6 β ν c 2 α 2 1 / 2 tan - 1 B - 1 B 2 - 1 2 .

As discussed above, in a non-rotating system, the right-hand side of Eq. (24) is zero, and the conservation of wave energy gives the following:

(27) E w 0 = 24 Q 2 α 2 β c 2 ν 3 1 / 2 ( B 2 - 1 ) 1 / 2 - 2 T ( B ) ,

where Ew0 is a constant evaluated at x=0. This can be solved to give B(x).

Radiation decay in a homogeneous environments (where c, α, … are constants) was recently considered by Obregon et al. (2018). Here we briefly redevelop their decay result in our variables for clarity.

In a homogeneous environment, the first term on the right-hand side of Eq. (25) or Eq. (26) is zero. Then the distance, XeG, to the complete radiation decay in the rG equation is found by the integration of Eq. (25) or Eq. (26), with the first term on the right side set to zero, from x=0 to XeG, where B=B0 and Be, respectively. When A0α>0, B0>0, the solitary wave, decays from an initial amplitude A0=αν-1(B0-1) to zero when Be=1, regardless of the sign of ν. While for ν>0 and A0α<0, B0<-1, and the solitary wave decay can be followed only to the limiting amplitude, Alim=-2α/ν, where Be=-1. These considerations give the following:

(28) X eG X eO = 2 8 I B ( B 0 , 1 ) | B 0 - 1 | - 1 / 2 , B 0 > 0 , I B ( B 0 , - 1 ) ( 1 - B 0 ) - 1 / 2 , B 0 < - 1 ,


(29) I B ( B 0 , B e ) = B 0 B e 1 - B 2 B tanh - 1 1 - B 1 + B - 2 d B , ν < 0 , - B 0 B e B 2 - 1 B tan - 1 B - 1 B + 1 - 2 d B , ν > 0 .

Again, XeO is the rKdV equation decay length (Eq. 20) evaluated using |A0α| when B0<-1.

Figure 1 shows XeG/XeO versus B0 for all three wave regimes. For ν<0, where 0<B0<1, XeG/XeO1 for 0.2<B0<1. There is a slight minimum of 0.9924 at B0=0.55. As B0→0, the ratio increases to XeG/XeO=1.1842 at B0=10-14 and appears to approach a finite limit as B0→0. For B0>1 (ν>0), XeG/XeO increases monotonically from one with B0 but remains O(1), even for B0=10. Similar behavior is found for B0<-1 (A0α>0), although it should be remembered that, in this regime, the adiabatic theory gives an amplitude only until the limiting wave at Be=-1 is reached.

Figure 2 shows examples of normalized wave amplitude A/A0 as functions of x/XeO for several values of B0 in each wave regime. For ν<0 and B00.1 (Fig. 2a), the amplitude decay closely follows the solution (Eq. 20) for the rKdV equation, but as B0 decreases, the initial amplitude decay rate slows. As mentioned above, for ν>0 and B0<-1 (Fig. 2c), the decay can only be followed until the limiting wave amplitude is reached.

Figure 4Numerical solution of Eq. (30) with (a) ϵ=2.5×10-5 and B0=10-4 and (b) ϵ=2.5×10-3 and B0=-4. The solid bar indicates the initial wave amplitude A0=sgn(ν)(B0-1). Only a portion of the full ξ domain is shown. The locations of each time series are indicated by the vertical axis.


4 Comparisons with rG numerical solutions

In this section, the adiabatic theory (Eq. 24) is compared with the full numerical solutions of the rG equation (Eq. 3). Example cases employ spatially uniform stratifications, and inhomogeneous effects are introduced by variations in the total water depth. The numerical solutions of Eq. (3) are found using a de-aliased pseudo-spectral scheme in s with a third-order Runge–Kutta integration in x. The relations for the coefficients α(x), ν(x), … are given in Appendix A.

4.1 Rotating homogeneous evolution

In the homogeneous case where the coefficients c, α, … are constants, it is convenient to reduce Eq. (3) to an equation with only one parameter by introducing the following scaling:




This change in variables gives the following:

(30) τ u ξ + u u τ + sgn ( ν ) u 2 u τ + 3 u τ 3 = ϵ u with ϵ = γ β ν 2 α 4 .

The solitary wave solution (Eqs. 68) and the adiabatic radiation decay results carry through after taking c=Q=α=β=1, ν=±1, and γ=ϵ. In these variables, XeO=ϵ-1(|A0|/12)1/2 and A0=sgn(ν)(B0-1).

Figure 3a–d show a comparison of the scaled amplitude A/A0 versus ξ/XeO from the adiabatic theory and full numerical solutions of the rG equation (Eq. 30) with B0= [0.3, 10−4, 5, −4], respectively, and 2.5×10-5ϵ0.25, as indicated. (The small oscillations in the amplitude are due to the periodic boundary conditions used in the numerical solutions that allowed the radiated waves to re-enter the domain upstream of the solitary wave.) For B0=10-4 and 0.3 (ν<0), the agreement between the adiabatic theory and the rG solutions is quite good for ϵ=2.5×10-5. However, for ϵ=2.5×10-3, there is disagreement. Similarly, for B0=-4 and 5, the agreement also declines with increasing ϵ, although, in these examples, the agreement for ϵ=2.5×10-3 is good. Increasing ϵ generally results in slower amplitude decay. The exception is B0=10-4, where the initial decay is more rapid. This rapid decay for near-maximal waves was also found by Obregon et al. (2018), which they attributed to a structural instability of large-amplitude, flat-top solitary waves.

The complicated evolution of the decaying solitary wave and the trailing radiation is illustrated in Fig. 4a and b. Figure 4a is the B0=10-4 and ϵ=2.5×10-5 example in Fig. 3b. As the initial solitary wave decays, the trailing radiation itself steepens to form a group of solitary-like waves that also decay by radiation damping. Over larger distances, this radiation will likely organize into one or more nonlinear wave packets, as found by Grimshaw and Helfrich (2008), for the rotating KdV equation, and Whitfield and Johnson (2015), in the rotating Gardner equation.

The example in Fig. 4b has B0=-4 and ϵ=2.5×10-3 (see Fig. 3d). For these parameters, the evolution is further complicated since the radiation decay ceases at ξ≈50 (=0.21XeO) when the wave amplitude decays to the limiting amplitude Alim=-2α/ν=-2 in these scaled variables. The wave then rapidly forms what appears to be a small solitary wave of reversed polarity (B0>1) and a trailing wave packet that has characteristics similar to the breather solutions of the Gardner equation. However, this packet subsequently disintegrates due to rotational effects. The complicated nature of the wave evolution with rotation in homogeneous conditions suggests even more interesting features when both rotation and inhomogeneous effects are active.

Figure 5Adiabatic theory for wave propagation from deep to shallow in a two-layer system, with h1=50 m, h20=450 m, g=0.005 m s−2, f=10-4 s−1, and A0=-25 m. A/A0, Ew/Ew0, and Mw/|Mw0| are shown as functions of x/L. The solid lines are for L=50 (black), 100 (blue), 200 (red), and 400 km (magenta). The dashed lines are for f=0. The dots are from the solutions of Eq. (3).


4.2 Combined inhomogeneous and rotation effects

Ostrovsky and Helfrich (2019) showed that, for the rKdV equation, the competition between extinction by radiation decay and the point of polarity reversal, α=0, could be characterized by the ratio of the inhomogeneous and rotation terms of Eq. (1) as follows:

(31) C c γ Q x Q η x η c γ L Δ = X eO L ,

where L is the length scale over which inhomogeneous term Q varies, e.g., the distance to the α=0 location. Δ is the solitary wave scale, which is taken above to be the KdV wave scale ΔKdV=(12β/A0α)1/2. For C≫1, inhomogeneous effects dominate and conversely rotational decay dominates for C≪1. Alternatively, Δ might taken to be the Gardner solitary wave scale ΔG=ΔKdV[2/(1+B0)]1/2. However, for B0>0 and not too large, the term in parentheses is O(1), and Eq. (31) is a reasonable scaling estimate. One could simply define C=XeG/L, but since XeG/XeO1 for 0<B0<10 (see Fig. 1), this also gives C from Eq. (31). The exception is for situations with B0<-1, since XeG/XeO can be much less than one.

To illustrate the combined effects of inhomogeneity and rotation, a two-layer Boussinesq stratification with upper layer depth h1, variable lower layer depth h2(x), reduced gravity g=g(ρ2-ρ1)/ρ1, and Coriolis frequency f will be considered. Here ρ1 and ρ2 are the densities of the upper and lower layers, respectively. Relations for the coefficients α(x), ν(x), etc., are given in (A6). Note that ν<0 for two-layer stratifications. Thus, only the 0<B0<1 branch of solitary wave solutions is possible. The wave polarity is given by the sign of α, with α<0 (>0) for h1/h2<1 (>1).

Figure 6rG numerical solutions for h1=50 m, h20=450 m, g=0.005 m s−2, f=10-4 s−1, and A0=-25 m. (a) L=50 km and (b) L=200 km. The wave amplitude η(x,t) is shown at x/L as a function of normalized shifted time, -s(g/H)1/2, where H=h1+h20. The initial wave amplitude is indicated by the scale on the lower right of each panel.


Figure 7η/H at x/L=1.2 for the parameters of Fig. 5. The bottom curve is the non-rotating case with L=50 km. The top four lines include rotation and have L=50, 100, 200, and 400 km from bottom to top. The curves are offset by 0.05.


The bottom slope will be taken constant, and the lower layer depth is given by the following:

(32) h 2 ( x ) h 20 = 1 - 1 - h 1 h 20 x L ,

where h20=h2(0), and L is the distance from the origin to the critical point, where α=0 (i.e., h1=h2). The case h1/h20<1 corresponds to an initial solitary wave, with A0<0 propagating from deep to shallower water. The propagation of a positive wave, A0>0, from shallow to deeper water occurs for h1/h20>1. While the adiabatic solutions can be obtained only up to the critical point, xL, numerical solutions of the rG equation are found beyond the critical depth. In the deep-to-shallow situation, the bottom slope is continued until h2=h1/2, beyond which h2 is constant over a shelf region.

Note that the linear bottom slope in Eq. (32) will not allow for the interesting, but rather special, topographic conditions that Stepanyants (2019) found for the rKdV regime. In those cases, rotation and inhomogeneity effects can be balanced and give a constant soliton amplitude. However, there is no balance for fluid velocity and, more important, soliton energy, which still decreases due to radiation.

Figure 8Adiabatic theory for wave propagation from deep to shallow in a two-layer system with h1=100 m, h20=200 m, g=0.01 m s−2, f=10-4 s−1, and L=40 km. Panels (a, c, e) show A/A0, Ew/Ew0 and Mw/|Mw0|, for the rotating Gardner theory. The solid (dashed) lines are with (without) rotation for A0=-10 m (black) and −45 m (blue). The open (solid) symbols are from the corresponding numerical solutions of Eq. (3) with (without) rotation. Panels (b, d, f) show the equivalent rotating KdV adiabatic solutions.


Figure 9η/H at x/L=1.5 for the parameters of Fig. 8. (a) A0=-10 m and (b) A0=-45 m. The solid (dashed) lines are with (without) rotation.


The evolution of the wave amplitude A(x) is shown in Fig. 5a for a deep to shallow case for a representative oceanographic situation, with h1=50 m, h20=450 m, g=0.005 m s−2, f=10-4 s−1, and A0=-25 m. For this initial wave, B0=0.764 so that effects of the cubic nonlinearity are present but do not dominate initially. Slope lengths L=50,100,200, and 400 km are considered and give C=3.46,1.73,0.86, and 0.43, respectively, from Eq. (31). The solid lines show the rotating adiabatic theory, and the dashed line shows the f=0 solution (equivalent for all L when plotted against x/L). As anticipated from the values of C, rotational decay effects increase significantly as L increases with complete, or nearly complete, extinction for the two longer slopes. However, even for L=50 km, rotation causes a noticeable reduction in wave amplitude compared to the non-rotating solution.

Figure 5b shows the wave energy, Ew(x) from Eq. (14), for the same parameters. The ratio Ew/Ew0 is a measure of the fraction of initial wave energy that remains in the evolving solitary wave, with the difference lost to the trailing inertia–gravity wave radiation. Even for the shortest slope, L=50 km, where the effects of rotation were relatively weak, more than half of the initial wave energy is lost to inertia–gravity wave radiation by x/L0.8.

The variation in wave mass, Mw from Eq. (15), is plotted in Fig. 5c. In all cases, rotation causes an initial reduction in wave mass, which is compensated in the trailing radiation and emerging shelf (with the same mass sign as the initial wave). The mass goes to zero for the two longer slopes prior to the critical point, while for the two shorter slopes the wave mass increases rapidly as the critical point is approached, similar to the non-rotating result. However, increased slope length (i.e., rotational effects) delays this growth with consequences for the magnitude of the trailing shelf and the subsequent wave structure transmitted onto the topographic shelf (see below).

Also shown in Fig. 5a and b are the results from the numerical solutions to the rG equation (Eq. 3). The energy of the solitary wave in the rG model is found by integrating Qη2 in a small region encompassing the solitary wave that does not incorporate an appreciable trailing radiation. The agreement is quite good for all the cases, except for x/L0.8 and L=50 and 100, where the amplitude from the rG numerical calculation does not decay as rapidly as in the adiabatic model. This is consistent with previous studies without rotation (Grimshaw et al.1998b, 1999). As the point of polarity reversal is approached, the wave elongates to form a rarefaction, and the trailing shelf with opposite polarity emerges. Figure 5 shows that the disagreement is associated with rapid changes in the wave amplitude and mass. In this region, that assumption of slow variation of the inhomogeneity fails.

There are two examples of the full rG solutions for L=50 and 200 km shown in Fig. 6, and the rG solutions at x/L=1.2 for all five cases in Fig. 5 are compared in Fig. 7. Note that h2=h1/2 for x/L1.0625 on the shelf. Again, even for L=50 km, the rotation leads to a clear effect on the wave signal transmitted on to the shelf as the leading crest of the (weak) trailing inertia–gravity wave also steepens to form a second transmitted wave packet. The second packet is consistent with the breaking criteria obtained by Grimshaw et al. (2012) and was similarly noted in Grimshaw et al. (2014). Further increases in the rotation effects lead to an additional transmitted packet for L=100 km. For the two longest slopes, L=200 and 400 km, the transmitted signal becomes increasingly disorganized. Figure 6b illustrates the evolution leading to this outcome. In this example, the initial solitary wave is extinguished before the critical point is reached. However, the trailing inertia–gravity wave steepens to produce a solitary wave that is itself scattered through the critical point. The interaction of the scattered signal with the trailing radiation gives rise to the observed disorganization. Recall that integration for calculating the energy of the leading wave, EW=Qη2ds, shown in Fig. 5 by the red solid dots, only includes the leading solitary wave and not the secondary solitary wave that emerges from the trailing radiation.

In the examples above, the cubic nonlinearity was not an essential feature of the evolution. Indeed, the wave evolution is qualitatively similar to the rotating KdV solutions in Fig. 2 of Ostrovsky and Helfrich (2019). In Fig. 8, another example with h1=100 m, h20=200 m, g=0.01 m s−2, f=10-4 s−1, and L=40 km is shown. Initial wave amplitudes of A0=-10 and −45 m correspond to B0=0.788 and 0.0478, respectively. The larger wave is very close to the limiting amplitude Alim=-47.06 m. The competition parameters are C=4.56 and 9.68, respectively. The left column of Fig. 8 shows the rG adiabatic solutions, and for comparison, the right column shows the adiabatic solutions obtained from the rKdV theory (Eq. 18). Now the differences between the rG and rKdV solutions are substantial, especially for the large wave. The rKdV solution has this larger wave decaying more slowly than the smaller wave, while the rG model shows just the opposite. Similar to above, rotation, even for this relatively short slope (and hence weak rotational effect), causes significant energy loss in the primary wave as it climbs the slope. For A0=-45 m, the mass Mw remains finite at the critical point, indicating the generation of a relatively weak trailing shelf. The agreement between the rG adiabatic model and the full rG solutions for A(x) is not very good, although the qualitative prediction that the large wave should decay more rapidly is found in both cases. The origin of the disagreement is likely due to the lack of separation between the wave scales and scale of the inhomogeneity.

The transmitted signals at x/L=1.5 from the full rG numerical solutions with and without rotation are shown in Fig. 9. For both initial wave amplitudes, rotation leads to significant changes in the transmitted signal. This is especially true for A0=-45 m, where the transmitted packet without rotation is replaced by a single, broad wave emerging onto the constant depth shelf with a much weaker trailing signal. However, on the shelf Alim/H=0.0745, so that this leading wave must further adjust and is also subject to continued rotational decay.

5 Concluding remarks

In this paper we have continued a series of studies of nonlinear internal waves of a moderate amplitude in the shallow, stratified areas of the ocean (see the references in the Introduction). Based on the classical Korteweg–de Vries equation, we added the main factors, making the analysis closer to the physical reality, i.e., cubic nonlinearity, Earth's rotation, and sloping bottom. The interplay of these factors makes the problem rather complicated, both physically and mathematically. To better explain the qualitative effect of each of them, we first briefly reproduce the effect of rotation in the medium with quadratic nonlinearity (KdV with rotation), then that with both quadratic and cubic nonlinearities (Gardner with rotation), and, as the main content of this paper, the joint effect of rotation and inhomogeneity in the Gardner equation. The specific qualitative effect of the latter is the limiting soliton amplitude, and the corresponding increase in the wavelength so that the topography effect becomes especially important. Along with the approximate adiabatic approach, a direct numerical study of the rG equation was performed and confirmed the adiabatic theory and highlighted its limitations. This combined approach allowed us to demonstrate a rather complicated behavior of shoaling internal solitons. For example, whereas the soliton energy always decreases due to the radiation losses, the displacement amplitude and mass can increase in a shoaling wave at a finite distance due to the decrease in total depth. In turn, the oscillating tail reveals a complicated behavior that includes, in particular, the formation of nonlinear wave packets as in Grimshaw and Helfrich (2008), followed by an even more complex evolution. Another important point of bifurcation is that of the wave passage, through the condition h1=h2, when the polarity of a soliton is changed. According to the estimations, the interplay between the topography and rotation effects does exist under realistic oceanic conditions.

Future research should include further comparisons of the theory with observational results in different oceanic environments and extend the above results to the strongly nonlinear waves with rotation.

Appendix A: Coefficients of the rotating Gardner equation

The eigenvalue problem for the linear long wave phase speed cn and vertical structure function for the isopycnal displacement Φn(z) of vertical mode n is, for the Boussinesq and rigid-lid limits (after dropping the subscript n), as follows:

(A1) d d z c - u 2 d Φ d z + N 2 ( z ) Φ = 0 , Φ ( - h ) = Φ ( 0 ) = 0 .

The buoyancy frequency is as follows:

(A2) N 2 ( z ) = - g ρ 0 d ρ d z ,

where ρ(z) is the background density profile, u(z) is a background current in the direction of wave propagation, g is the gravitational acceleration, and ρ0 is a reference density. From Grimshaw et al. (2004), the coefficients of the rG equation (Eq. 1) for a particular mode are given by the following:


When u=0, γ=f2/2c, where f is the Coriolis frequency. The general relation for γ with u0 was derived by Alias et al. (2014). The coefficient of the cubic nonlinear term is given by the following:

(A4) ν = I - 1 - h 0 3 ( c - u ) 2 3 T z - 2 Φ 2 Φ 2 - α 2 Φ 2 + α ( c - u ) 5 Φ 2 - 4 T z Φ d z ,

with the nonlinear correction to the vertical structure function, T(z), found from

(A5) d d z c - u 2 d T d z + N 2 T = - α d d z ( c - u ) d Φ d z + 3 2 d d z ( c - u ) 2 d Φ d z 2 ,

with T(-h)=T(0)=0. The solution T(z) is normalized through the addition of bΦ(z), with b chosen so that T(zmax)=0, where Φ(zmax)=1. This gives the isopycnal vertical displacement ξ(x,z,t)=η(x,t)Φ(z)+η2(x,t)T(z).

In a two-layered stratification, with depths h1 and h2 of the upper and lower layers, respectively, and u=0, the following applies:

(A6) c 2 = g h 1 h 2 h 1 + h 2 , α = 3 c 2 h 1 - h 2 h 1 h 2 , ν = - 3 c 8 h 1 2 + 6 h 1 h 2 + h 2 2 ( h 1 h 2 ) 2 , β = c 6 h 1 h 2 , γ = f 2 2 c , Q = 2 g c .

Here g=gΔρ/ρ0, and Δρ=ρ2-ρ1 is the difference in densities between the lower and upper layers.

Code and data availability

This work does not include any externally supplied code, data, or other material. All material in the text and figures was produced by the authors using standard mathematical and numerical analysis by the authors.

Author contributions

KRH and LO equally contributed to formulating the problem, obtaining analytical solutions, and discussing the results. Numerical modeling was performed by KRH.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank Yury Stepanyants, for the helpful discussions.

Financial support

This research has been supported by the Office of Naval Research (grant no. N00014-18-1-2542) and National Science Foundation (grant no. OCE-1736698).

Review statement

This paper was edited by Victor Shrira and reviewed by two anonymous referees.


Alford, M., Lien, R.-C., Simmons, H., Klymak, J., Ramp, S., Yang, Y. J., Tang, D., and Chang, M.-H.: Speed and evolution of nonlinear internal waves transiting the South China Sea, J. Phys. Oceano., 121, 1338–1355,, 2010. a

Alias, A., Grimshaw, R. H. J., and Khusnutdinova, K. R.: Coupled Ostrovsky equations for internal waves in a shear flow, Phys. Fluids, 26, 126603,, 2014. a, b

Apel, J. R., Ostrovsky, L. A., Stepanyants, Y. A., and Lynch, J. F.: Internal solitons in the ocean and their effect on underwater sound, J. Acoust. Soc. Am., 121, 695–722,, 2007. a

Farmer, D., Li, Q., and Park, J.-H.: Internal wave observations in the South China Sea: the role of rotation and nonlinearity, Atmos.-Ocean, 47, 267–280,, 2009. a

Grimshaw, R., He, J.-M., and Ostrovsky, L.: Terminal damping of a solitary wave due to radiation in rotational systems, Stud. Appl. Math, 101, 197–210, 1998a. a, b

Grimshaw, R., Pelinovsky, E., and Talipova, T.: Solitary wave transformation due to a change in polarity, Stud. Appl. Math, 101, 357–388, 1998b. a, b

Grimshaw, R., Pelinovsky, E., and Talipova, T.: Solitary wave transformation in a medium with sign-variable quadratic nonlinearity and cubic nonlinearty, Physica D, 132, 40–62, 1999. a, b, c

Grimshaw, R., Pelinovsky, E., Talipova, T., and Kurkin, A.: Simulation of the transformation of internal solitary waves on oceanic shelves, J. Phys. Oceano., 34, 2774–2791, 2004. a, b, c

Grimshaw, R., Helfrich, K., and Johnson, E. R.: The reduced Ostrovsky equation: Integrability and breaking, Stud. Appl. Math., 129, 414–436, 2012. a

Grimshaw, R., Guo, C., Helfrich, K., and Vlasenko, V.: Combined effect of rotation and topography on shoaling oceanic internal solitary waves, J. Phys. Oceano., 44, 1116–1132, 2014. a, b, c, d, e

Grimshaw, R. H. J. and Helfrich, K. R.: Long-time solutions of the Ostrovsky equation, Studies Appl. Math., 121, 71–88, 2008. a, b

Helfrich, K. R.: Decay and return of internal solitary waves with rotation, Phys. Fluids., 19, 026601, 2007. a

Helfrich, K. R. and Grimshaw, R.: Nonlinear disintegration of the internal tide, J. Phys. Oceano., 38, 686–701,, 2008. a

Holloway, P. E., Pelinovsky, E., and Talipova, T.: A generalized Korteweg-de Vries model of internal tide transformation in the coastal ocean, J. Geophys. Res., 104, 18333–18350, 1999. a, b, c

Jackson, C.: An Atlas of Internal Solitary-like Waves and Their Properties, 2nd ed., 560 pp., Global Ocean Assoc., Alexandria, Va., (last access: 8 June 2022), 2004. a

Karczewska, A. and Rozmej, P.: Can simple KdV-type equations be derived for shallow water problem with bottom bathymetry?, Comm. Nonlin. Sci. Numer. Sim., 82, 105073,, 2020. a

Nakoulima, O., Zahibo, N., Pelinovsky, E., Talipova, T., Sulnyaev, A., and Kurkin, A.: Analytical and numerical studies of the variable-coefficient Gardner equation, Appl. Math. Comp., 152, 449–471, 2004. a

Obregon, M., Raj, N., and Stepanyants, Y.: Adiabatic decay of internal solitons due to Earth's rotation within the framework of the Gardner-Ostrovsky equation, Chaos, 28, 033106,, 2018. a, b, c

Ostrovsky, L.: Nonlinear internal waves in a rotating ocean, Oceanology, 18, 119–1125, 1978. a, b

Ostrovsky, L. and Grue, J.: Evolution equations for strongly nonlinear internal waves, Phys. Fluids, 15, 2934–2948,, 2003. a

Ostrovsky, L. and Helfrich, K. R.: Some new aspects of the joint effect of rotation and topography on internal solitary waves, J. Phys. Oceano., 49, 1639–1649,, 2019. a, b, c, d, e

Ostrovsky, L. and Irisov, V. G.: Strongly nonlinear internal solitons: Models and applications, J. Geophys. Res., 122, 3907–3916,, 2017.  a

Pelinovsky, E. N. and Grimshaw, R.: Structural transformation of eigenvalues for a perturbed algebraic soliton potential, Phys. Lett. A, 229, 165–172, 1997. a

Pelinovsky, E. N., Polukhina, O. E., and Lamb, K.: Nonlinear internal waves in the ocean stratified in density and current, Marine Phys., 40, 757–766, 2000. a

Stanton, T. P. and Ostrovsky, L. A.: Observations of highly nonlinear internal solitons over the continental shelf, Geophys. Res. Let., 25, 2695–2698, 1998. a

Stepanyants, Y.: The effects of interplay between the rotation and shoaling for a solitary wave on variable topography, Stud. Appl. Math., 142, 465–486,, 2019. a, b, c

Talipova, T. G., Kurkina, O. E., Rouvinskaya, E. A., and Pelinovsky, E. N.: Propagation of solitary internal waves in two-layer ocean of variable depth, Izvestiya, Atmos. Ocean. Phys., 51, 89–97, 2015. a

Whitfield, A. J. and Johnson, E. R.: Wave-packet formation at the zero-dispersion point in the Gardner-Ostrovsky equation, Phys. Rev., 91, 051201,, 2015. a

Zhao, Z. and Alford, M.: Source and propagation of nonlinear internal waves in the northeastern South China Sea, J. Geophys. Res., 111, C11012,, 2006. a

Short summary
Internal solitons are an important class of nonlinear waves commonly observed in coastal oceans. Their propagation is affected by the Earth's rotation and the variation in the water depth. We consider an interplay of these factors using the corresponding extension of the Gardner equation. This model allows a limiting soliton amplitude and the corresponding increase in wavelength, making the effects of rotation and topography on a shoaling wave especially significant.