2.1 The case $\unicode[STIX]{x1D6FC}\ll 1$

Working on the assumption that $\unicode[STIX]{x1D6FC}$ is small, we seek an asymptotic solution to (2.3) in the form

Substituting (2.5) into (2.3) and solving at successive orders of approximation using the first of the conditions in (2.4), we find

where the $\overline{u}_{i}$ are constants of integration, and $\unicode[STIX]{x1D707}(\unicode[STIX]{x1D70F})=\int _{0}^{\unicode[STIX]{x1D70F}}\text{exp}(-{\unicode[STIX]{x1D70F}^{\prime }}^{2})\hspace{2.22198pt}\text{d}\unicode[STIX]{x1D70F}^{\prime }$ . In the far field we find that

as $\unicode[STIX]{x1D709}\rightarrow \infty$ , which indicates that the expansion will fail when $\unicode[STIX]{x1D709}\sim O(\unicode[STIX]{x1D6FC}^{-1/3})$ . Accordingly we introduce an outer region in which the new variable $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6FC}^{1/3}\unicode[STIX]{x1D709}=O(1)$ and posit the expansion $u=\unicode[STIX]{x1D6FC}^{2/3}\tilde{u} _{0}(\unicode[STIX]{x1D701})+\cdots \,$ inside this outer region. Substituting into (2.3) we find at leading order and at first order,

The solutions which satisfy the second condition in (2.4) are

where $d_{1}$ , $d_{2}$ are constants of integration. Matching to the inner region we find that

as well as $d_{2}=0$ and $\overline{u}_{1}=1/2$ . We conclude that for small $\unicode[STIX]{x1D6FC}$ , there exists a solution with negative elevation at the point of symmetry with $u(0)<0$ .

2.2 The case $\unicode[STIX]{x1D6FC}\gg 1$

In this case it will be convenient to rescale the problem by introducing the new dependent variable $u=\unicode[STIX]{x1D6FC}^{1/2}U$ so that (2.3) takes the form

where $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FC}^{-1/2}$ and $\unicode[STIX]{x1D6FF}\ll 1$ . Referring to figure 2 we see that the asymptotic solution for small $\unicode[STIX]{x1D6FF}$ has four distinct regions. In the core region, labelled region ii in the figure, we expand the solution by writing

Substituting into (2.11) we obtain algebraic equations for the $U_{i}$ at successive orders of approximation which are solved easily. At leading order we find

Discussion of the choice of the sign in (2.13) will be given below. At $O(\unicode[STIX]{x1D6FF})$ we have

which is valid regardless of the sign taken for $U_{0}$ in (2.13). From the next two orders, we determine the forms

Although each of the $U_{i}$ in (2.12) satisfies the symmetry condition $U_{i\unicode[STIX]{x1D709}}(0)=0$ required by (2.4), the expansion must nonetheless break down as $\unicode[STIX]{x1D709}\rightarrow 0$ since the first term on the left-hand side of (2.11) will become comparable to the second in this limit. Moreover, (2.12) evidently breaks down as $\unicode[STIX]{x1D709}$ increases since successive terms in the expansion grow to become comparable in size to their predecessors. We address the first of these issues first, and consider the solution in the boundary layer at $\unicode[STIX]{x1D709}=0$ , namely region i in figure 2.

Figure 2.

(a) Sketch of the four regions used in the matched asymptotics for large $\unicode[STIX]{x1D6FC}$ . The decay is exponential in region ii and algebraic in region iv. (b) Numerical solution to the region iii problem (2.27) subject to (2.28) with the minus sign.

In region i, we introduce the boundary-layer variable $y=\unicode[STIX]{x1D6FF}^{-1/2}\unicode[STIX]{x1D709}$ and expand the solution by writing

Substituting into (2.11), at leading order we find

subject to the symmetry condition $\tilde{U} _{0y}(0)=0$ and the condition that $\tilde{U} _{0}\sim \pm 1$ as $y\rightarrow \infty$ to match with region ii, the choice of sign in the latter condition being dictated by the choice of sign taken in (2.13). This problem has three possible solutions given by

Referring to the $(\tilde{U} _{0},\tilde{U} _{0y})$ phase plane, these correspond respectively to fixed points at $\tilde{U} _{0}=\pm 1$ , and a homoclinic orbit encircling $\tilde{U} _{0}=1$ which connects $\tilde{U} _{0}=-1$ to itself.

We examine cases (i) and (ii) first. At $O(\unicode[STIX]{x1D6FF})$ we have

The solution which satisfies $\tilde{U} _{1y}(0)=0$ , and which does not grow exponentially as $y\rightarrow \infty$ , is given by

Since the expressions in (2.20) are forced directly by the solution in the core region ii, the matching with that region is immediate. For case (iii), we have

By making the change of variables $t=\tanh (y/\sqrt{2})$ , this can be transformed into the associated Legendre equation

The solution which satisfies the required symmetry condition $\tilde{U} _{1t}(0)=0$ is given by

where

are the homogenous solutions to (2.22). We note that $\tilde{U} _{1}(t)\sim \operatorname{arctanh}^{2}t$ as $t\rightarrow 1$ and so $\tilde{U} _{1}\sim y^{2}/2$ ensuring a match with region ii.

We now turn our attention to region iii. Scrutinising the terms $U_{0}$ and $U_{1}$ in (2.13) and (2.14), we see that the expansion (2.12) in region ii breaks down when $\exp (-\unicode[STIX]{x1D709}^{2}/2)\sim 1/2\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D709}^{2}-1)$ . This occurs when $\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D709}_{0}(\unicode[STIX]{x1D6FF})$ , where

with $\unicode[STIX]{x1D706}=(-2\log \unicode[STIX]{x1D6FF})^{1/2}$ . Inside this narrow region, we introduce the internal layer variable $z=\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}_{0})$ , and expand by writing

Substituting into equation (2.11), we obtain at leading order,

The boundary conditions are

Here the $+/-$ sign is required to match to the solution in region ii, namely (2.13) with the $+/-$ sign selected respectively. It is possible to compute a solution to this problem numerically for either choice of sign; however, great care is needed when attempting to match the solution with region ii, as will be discussed below. We discuss first the case of the minus sign in (2.28). Figure 2 shows the result of a calculation using fourth-order Runge–Kutta integration to solve the problem (2.27) with the minus sign chosen in the boundary condition (2.28). An analysis for large negative $z$ reveals that the solution takes the form

as $z\rightarrow -\infty$ , where $\text{K}_{2}$ is a modified Bessel function of the second kind, and $c_{1}$ is a constant which must be determined globally but is not required here. Since it is known that $\text{K}_{2}(w)\sim \sqrt{\unicode[STIX]{x03C0}/2z}\,\text{exp}(-2z)\cdots \,$ (e.g. Abramowitz & Stegun Reference Abramowitz and Stegun1964), the final term in (2.29) decays superexponentially and matching to region ii is possible. For large positive $z$ , the forcing term on the right-hand side of (2.27) becomes negligible and the dominant balance of the two terms on the left-hand side requires the behaviour $M_{0}\sim -6/(z+c_{2})^{2}$ , for constant $c_{2}$ , as $z\rightarrow \infty$ . In region iv, the expansion proceeds as

where $w=\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}_{0}=O(1)$ . At leading order we have $\hat{U} _{0ww}+U_{0}^{2}=0$ and the solution satisfying the far-field condition in (2.4), namely $\hat{U} _{0}\rightarrow 0$ as $w\rightarrow \infty$ , is

for constant $c_{3}$ which can in principle be found by matching to region iii.

As noted above, it is possible to numerically compute a solution to the problem (2.27) subject to (2.28) with the plus sign selected in the boundary condition. However, the numerical solution exhibits small decaying oscillations for large negative $z$ . Indeed, in the limit $z\rightarrow -\infty$ it can be shown that

for constants $\unicode[STIX]{x1D6FC}_{1}$ , $\unicode[STIX]{x1D6FC}_{2}$ , where $\text{J}_{2}$ and $\text{Y}_{2}$ are Bessel functions of the first and second kind respectively. As $\unicode[STIX]{x1D709}\rightarrow \infty$ , $\text{J}_{2}(\unicode[STIX]{x1D709})\sim \sqrt{2/\unicode[STIX]{x03C0}\unicode[STIX]{x1D709}}\,\cos (\unicode[STIX]{x1D709}-5\unicode[STIX]{x03C0}/4)+\cdots \,$ and $\text{Y}_{2}(\unicode[STIX]{x1D709})\sim \sqrt{2/\unicode[STIX]{x03C0}\unicode[STIX]{x1D709}}\,\sin (\unicode[STIX]{x1D709}-5\unicode[STIX]{x03C0}/4)+\cdots \,$ (Abramowitz & Stegun Reference Abramowitz and Stegun1964) and so $M_{0}$ has decaying oscillations as $z\rightarrow -\infty$ in agreement with the numerical solution. Since the solution in region ii is non-oscillatory, it is not possible to complete a match between regions ii and iii, and on this basis we must reject the choice of the plus sign in (2.13) since it does not lead to a consistent asymptotic solution.

We conclude that consistent asymptotic solutions may be constructed for cases (ii) and (iii) in (2.18), but not for case (i). This means that there are two asymptotic solutions consistent with the choice of a minus sign in (2.13), namely cases (ii) and (iii), but there are no solutions consistent with the choice of a plus sign in (2.13), namely case (i).

2.3 Numerical results

In summary, we have found in § 2.2 that it is possible to construct a consistent asymptotic solution using boundary-layer theory for cases (ii) and (iii) in (2.18), but not for case (i). This together with the results of the small $\unicode[STIX]{x1D6FC}$ asymptotic analysis in § 2.1 suggests that there exist at least two distinct branches of steady solutions. To carry out a full investigation of the solution space, we solved the steady problem (2.3) numerically employing a shooting technique based on the fourth-order Runge–Kutta method. The problem was first reformulated as a first-order system in the form $\boldsymbol{x}_{\unicode[STIX]{x1D709}}=\boldsymbol{F}(\boldsymbol{x})$ , where $\boldsymbol{x}=(u,u_{\unicode[STIX]{x1D709}})^{\text{T}}$ and $\boldsymbol{F}=(u_{\unicode[STIX]{x1D709}},\unicode[STIX]{x1D6FC}s-u^{2})^{\text{T}}$ . Following the assumed symmetry of the solution, the system was integrated forward from the starting point $\boldsymbol{x}(\unicode[STIX]{x1D709}=0)=(u_{0},0)^{\text{T}}$ for a chosen $u_{0}$ . The correct value of $u_{0}$ was determined iteratively until the trajectory of the solution in the $(u,u_{\unicode[STIX]{x1D709}})$ phase plane approaches the origin. Formally this should happen in infinite time as $\unicode[STIX]{x1D709}\rightarrow \infty$ but in computational practice the value of $u_{0}$ was refined so that an increasingly large value of $\unicode[STIX]{x1D709}$ was required to reach the origin. There are two qualitatively different ways in which this can happen as follows:

as $\unicode[STIX]{x1D709}\rightarrow \infty$ . Option (a) represents a balance between the two terms on the left-hand side of (2.3) and constitutes algebraic decay of the solution as $-6\unicode[STIX]{x1D709}^{-2}$ for large $\unicode[STIX]{x1D709}$ . In this case the slope of the trajectory tends to $0-$ as the origin is approached. Option (b) represents a balance between the first term on the left-hand side of (2.3) and the topographic forcing term on the right-hand side. In this case the trajectory enters the origin from below with an infinitely negative slope. The computed solution was checked against an independent calculation in which the boundary value problem was solved using a finite difference method on a uniform grid over the truncated domain $[-L,L]$ with $L$ taken to be suitably large. Approximating the second derivative in (2.3) with a centred difference furnishes a set of nonlinear algebraic equations at each grid point, which were solved using Newton’s method.

Figure 3(a–c) shows numerical computations of the two solution types predicted by the asymptotic analysis. We refer to the depression wave in figure 3(a) as a type i solution and the elevation waves in figure 3(b,c) as type  ii solutions. The solution space, characterised by the value $u(0)$ , is shown in figures 3(h) and 3(i). The predictions of the asymptotic analysis are included in figures 3(h) and 3(i). Numerically, we find that the far-field behaviour (a) in (2.33) is obtained on the lower branch of type i solutions and also on almost all of the upper branch except at the point $\unicode[STIX]{x1D6FC}\approx 35.14$ , where this branch terminates and where the numerical evidence appears to suggest that behaviour (b) occurs. Further experimentation revealed the existence of a seemingly infinite sequence of new solutions iia,b,c,d etc., which are not predicted by our asymptotic analysis, and which have oscillations on the free surface localised around the origin. Examples are shown in figure 3(d–g). We note the similarity with the solution space over a rectangular trench (Keeler Reference Keeler2017).

Figure 3.

(a–g) Free-surface profiles: (a) type i at $\unicode[STIX]{x1D6FC}=50$ ; (b) type ii at $\unicode[STIX]{x1D6FC}=100$ ; (c) type ii at $\unicode[STIX]{x1D6FC}=35.2$ , close to the branch termination point; (d) type iia at $\unicode[STIX]{x1D6FC}=300$ ; (e) type iib at $\unicode[STIX]{x1D6FC}=1500$ ; (f) type iic at $\unicode[STIX]{x1D6FC}=3000$ ; (g) type iid at $\unicode[STIX]{x1D6FC}=5000$ . (h) Solution space. The large $\unicode[STIX]{x1D6FC}$ behaviour $u(0)\sim -\unicode[STIX]{x1D6FC}^{1/2}+1/2$ , case (ii) in (2.18), is shown with a broken curve overlaid on the lower branch. The broken curve overlaid on the upper branch is the large $\unicode[STIX]{x1D6FC}$ behaviour for the type ii solution $u(0)\sim 2\unicode[STIX]{x1D6FC}^{1/2}-4\log (2)$ , case (iii) in (2.18). The inset shows behaviour for small $\unicode[STIX]{x1D6FC}$ . The broken curve illustrates the small $\unicode[STIX]{x1D6FC}$ behaviour $u(0)\sim -(6/d_{1}^{2})\unicode[STIX]{x1D6FC}^{2/3}+(1/2)\unicode[STIX]{x1D6FC}$ , with $d_{1}$ as given in (2.10). (i) Extended view of the solution space including type iia, b, c and d solutions. The latter solutions terminate at $\unicode[STIX]{x1D6FC}\approx 214.0,\hspace{2.22198pt}699.1,\hspace{2.22198pt}1709.0$ and $3508.0$ respectively.

Figure 4.

Fully nonlinear results for the Gaussian topography (solid curves), $b=0.3$ . The broken curves are for the weakly nonlinear solutions. (a) Type i solution, $a=-0.0017$ . (b) Type ii solution, $a=-0.039$ . (c) Type iia solution, $a=-0.134$ . (d) Solution space for the type i and ii solutions. (e) Solution space for the type iia solution. (f) Reduced surface angle $\unicode[STIX]{x1D703}/\unicode[STIX]{x1D703}_{c}$ , where $\unicode[STIX]{x1D703}_{c}=\unicode[STIX]{x03C0}/6$ , for the highest computed type ii wave at $a=-0.06484$ . In (f) the inset shows a close-up near to $x=0$ .

2.4 Fully nonlinear calculations

We now compare the results obtained for the weakly nonlinear model to those found for the fully nonlinear problem. In the non-dimensionalised steady water-wave equations the velocity potential $\unicode[STIX]{x1D719}(x,y)$ satisfies Laplace’s equation throughout the flow domain,

On the free surface the dynamic and kinematic boundary conditions need to be satisfied,

along with the no-penetration condition on the bottom,

We solved this system of equations for the unknown functions $\unicode[STIX]{x1D719}(x,y)$ and $\unicode[STIX]{x1D702}(x)$ at $F=1$ numerically for the Gaussian forcing $f(x)=a\exp (-b^{2}x^{2})$ with the aim of computing solitary wave solutions. Based on the weakly nonlinear analysis we would expect the tails of the solitary waves to decay algebraically, and indeed we establish analytically in appendix A that the full nonlinear solutions decay like $-(4/3)x^{-2}$ as $x\rightarrow \pm \infty$ , which is precisely the far-field behaviour found for the fKdV equation (2.3). Consequently a large truncation length for the computational domain is required for the computations. The numerical computations were performed first using the boundary-integral method described by Tam et al. (Reference Tam, Yu, Kelso and Binder2015) over the truncated domain $[-L,L]$ for chosen $L$ . Anticipating an inverse square algebraic decay in the far field, the boundary conditions $\unicode[STIX]{x1D702}_{x}/(1+\unicode[STIX]{x1D702})=\mp \,2/L$ were applied at $x=\pm L$ respectively. Second, the fully nonlinear system is first reformulated as a single partial differential equation with a non-local term in terms of surface variables alone, as described in appendix A. The resulting equation (A 11) together with (A 5) is solved with spectral accuracy over the truncated domain $[-L^{\prime },L^{\prime }]$ , for large $L^{\prime }$ in the mapped domain, assuming periodic boundary conditions, by employing the fast Fourier transform. In the results presented here the computations were done with the first of the numerical approaches, except for the type II branch solutions, where the second approach proved more expeditious particularly in obtaining the steep wave profiles.

The comparison with the weakly nonlinear solutions for the steady fKdV equation (2.3) is shown in figure 4(a–e). The type i solution demonstrates a very good agreement with a small discrepancy at the wave minimum. Note that only a portion of the wave is shown around the origin since the two solutions are visually indistinguishable over the remainder of the range. There is very good agreement on the type i solution branch for small $|a|$ as expected. The fully nonlinear type ii solution branch exhibits similar behaviour to the weakly nonlinear prediction but while the latter branch continues indefinitely through negative values of $a$ , the former terminates at a finite value of $a$ . We note that the fully nonlinear branch cannot possibly follow the weakly nonlinear branch for large $a$ since it is clear from the Bernoulli equation in (2.35) that the surface elevation is bounded above with $\unicode[STIX]{x1D702}\leqslant 1/2$ in the case $F=1$ ; equality is obtained at a stagnation point where the first term in the Bernoulli equation vanishes. Such an eventuality is expected at the crest of a Stokes-limiting configuration, as is found for unforced solitary waves (e.g. Hunter & Vanden-Broeck Reference Hunter and Vanden-Broeck1983). We computed the branch up to $a=-0.06484$ (see figure 4 d) at which point a very large number of mesh points are required to obtain an accurate solution as the wave profile does indeed seem to be approaching a Stokes-limiting configuration with an interior angle of $120^{\circ }$ , and a stagnation point, at the crest. The surface angle $\unicode[STIX]{x1D703}$ , defined as the angle between the tangent to the surface and the horizontal, is shown in figure 4(f) for the type ii profile obtained at $a=-0.06484$ . This particular calculation was performed using the numerical method discussed above with $2^{14}=16\,384$ points equally spaced in the mapped domain with $L^{\prime }=100$ . Given the very high resolution required, we did not attempt to move further along the branch toward the Stokes-limiting configuration. Wade et al. (Reference Wade, Binder, Mattner and Denier2014) experienced similar numerical difficulties close to the Stokes-limiting configuration for an imposed distribution of pressure on the free surface. In the current work the algebraic far-field decay presents an added difficulty.

Figure 5.

Far-field decay of fully nonlinear solutions for a Gaussian topography, $a=-0.03$ and $b=0.3$ , computed on domain $[-L^{\prime },L^{\prime }]$ with $2048$ collocation points using the method described in appendix A. Profiles for type i (a) and type ii (b) shown on a $\log$ - $\log$ plot. In both (a) and (b) the asymptote $\log (-\unicode[STIX]{x1D702})=\log (4/3)-2\log x$ is shown with a dashed line.

As was discussed in § 2.3, the weakly nonlinear branch terminates at its right-hand end at a negative value of $a$ , as can be seen in figure 4(d). We find that the fully nonlinear branch also appears to terminate at its right-hand end. Generally solutions along this fully nonlinear branch decay algebraically as $-(4/3)x^{-2}$ and this is confirmed in figure 5, which demonstrates excellent agreement between the analytic prediction and the numerical computation for a sample case along this branch. It is numerically extremely difficult to identify the shift from algebraic decay to exponential decay (behaviour (b) in (2.33)) which we anticipate will herald the end of the branch as it does for the weakly nonlinear theory.

The solution type iia also exists in the fully nonlinear solution space but there is a large difference between the fKdV and fully nonlinear solution branches (see figure 4 c,e) as might be expected since the amplitude $a$ is here rather large. It seems reasonable to suppose that the type iib,c, $\ldots$ branches are also present in the fully nonlinear regime but since they would likely be quantitatively quite different to their weakly nonlinear counterparts it may not be straightforward to find them, and we have not attempted to do so. We also note in common with the weakly nonlinear theory, the type iia branch terminates at its right-hand end; however as with the type ii solutions it appears also to terminate at its other end.

3.1 Stability analysis

The aim is to classify the steady solutions we have computed as being neutrally, spectrally, linearly or formally stable. For a definition of these terms, see for example, Holm et al. (Reference Holm, Marsden, Ratiu and Weinstein1985); however we will remind the reader of their precise meaning as the analysis proceeds.

Our starting point is the time-dependent fKdV equation (1.2) at the critical condition $F=1$ , which is repeated here for convenience,

To investigate the stability we write $\unicode[STIX]{x1D702}(x,t)=\unicode[STIX]{x1D702}_{s}(x)+\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}(x,t)$ , where $\unicode[STIX]{x1D702}_{s}$ is a steady solution of (3.1) and $\unicode[STIX]{x1D701}$ is a time-dependent perturbation, and where $\unicode[STIX]{x1D716}\ll 1$ . We demand that the perturbation $\unicode[STIX]{x1D701}$ satisfies the same far-field conditions as $\unicode[STIX]{x1D702}$ , namely equation (1.3). Substituting into (3.1) and retaining terms at $O(\unicode[STIX]{x1D716})$ we obtain

The general solution to (3.2) may be written in the form (e.g. Chang & Demekhin Reference Chang and Demekhin2002, § 7.3)

where the $a_{k}$ and $a(s)$ are expansion coefficients corresponding respectively to the $N_{p}$ values of $s_{k}$ in the point spectrum and to the continuous spectrum of the operator ${\mathcal{M}}$ . These satisfy the relation

(with $k$ subscripts appended in the case of the point spectrum). The integral in the second term in (3.3) is taken over the contour $\unicode[STIX]{x1D6E4}$ which delineates the continuous spectrum in the complex $s$ plane. The eigenfunctions in the point spectrum, $g_{k}(x)$ , reside in $L_{2}(-\infty ,\infty )$ . Furthermore, they satisfy the integral constraint,

which is obtained by integrating (3.4) over the real line and using the fact that $g_{k}\rightarrow 0$ as $|x|\rightarrow \infty$ (Camassa & Wu Reference Camassa and Wu1991). Consequently the eigenfunctions $g_{k}$ carry zero excess mass. The functions $g(x,s)$ associated with the continuous spectrum $\unicode[STIX]{x1D6E4}$ are bounded as $|x|\rightarrow \infty$ . The continuous spectrum itself is determined by examining the behaviour as $|x|\rightarrow \infty$ (e.g. Sandstede Reference Sandstede2002). In this limit, the second term in ${\mathcal{K}}$ defined in (3.2) can be neglected and (3.2) reduces to $\unicode[STIX]{x1D701}_{t}-(1/6)\unicode[STIX]{x1D701}_{xxx}=0$ ; adopting the form $\unicode[STIX]{x1D701}=\text{exp}(st)g(x,s)$ , we find $g=\text{exp}(\text{i}kx)$ and $s=-\text{i}k^{3}/6$ for $k\in \mathbb{R}$ so that the continuous spectrum occupies the imaginary axis: $\unicode[STIX]{x1D6E4}=\{\text{i}q:q\in \mathbb{R}\}$ .

In the case of type i steady solutions, which have the property that $\unicode[STIX]{x1D702}_{s}(x)<0$ for all $x$ , it is straightforward to show that the point spectrum of ${\mathcal{M}}$ does not include the origin (Keeler Reference Keeler2017). In the case of type ii steady solutions, numerically computed spectra are found to not include the origin, as will be discussed below.

If there is no eigenvalue in the spectrum in the right-half-complex-plane, that is with a positive real part, then $\unicode[STIX]{x1D702}_{s}(x)$ is spectrally stable. If the spectrum is contained along the imaginary axis, then the solution is neutrally stable. The solution is linearly stable if for every $\unicode[STIX]{x1D716}>0$ there exists a $\unicode[STIX]{x1D6FF}>0$ such that if in a suitable norm $\Vert \unicode[STIX]{x1D701}(x,0)\Vert <\unicode[STIX]{x1D6FF}$ then $\Vert \unicode[STIX]{x1D701}(x,t)\Vert <\unicode[STIX]{x1D716}$ for $t>0$ . While linear stability clearly implies spectral stability, the converse is not necessarily true (see Holm et al. (Reference Holm, Marsden, Ratiu and Weinstein1985) for some examples). As was pointed out by Camassa & Wu (Reference Camassa and Wu1991), for a steady solution $\unicode[STIX]{x1D702}_{s}(x)$ which is even in $x$ (resulting from an even topography), the operator ${\mathcal{M}}$ is odd in $x$ , i.e. ${\mathcal{M}}g(-x)=-{\mathcal{M}}g(x)$ , which means that for each eigenvalue $s$ satisfying (3.4), $-s$ , $s^{\ast }$ and $-s^{\ast }$ are also eigenvalues, so that the spectrum has a fourfold symmetry in the complex plane. As a consequence, spectral stability can only occur if all of the eigenvalues lie on the imaginary axis, in which case the solution is neutrally stable. Therefore for the fKdV equation with even topography, neutral stability for even solutions is a necessary and sufficient condition for spectral stability.

We have computed the spectrum of the ${\mathcal{M}}$ operator numerically on a truncated domain $[-L,L]$ with periodic boundary conditions using a spectral method. Sandstede & Scheel (Reference Sandstede and Scheel2000) showed that calculations performed on a periodic domain faithfully reproduce the point and continuous spectra provided the size of the domain is large enough. Bridges, Derks & Gottwald (Reference Bridges, Derks and Gottwald2002) and Barashenkov & Zemlyanaya (Reference Barashenkov and Zemlyanaya2000) have noted that numerical algorithms for computing eigenvalues in KdV-type problems can lead to spurious eigenvalues leaving/entering the imaginary axis and posing as unstable/stable modes. With this in mind, we have performed careful convergence tests on our results and, furthermore, we confirm our results by comparison with time-dependent simulations to be discussed in the next subsection. For type i solutions the spectrum is found to be contained on the imaginary axis indicating that the flow is spectrally stable in this case. This conclusion will be confirmed via a formal argument to be presented below. In contrast type ii solutions are found numerically to be linearly unstable. In all cases we examined we found that the operator ${\mathcal{M}}$ has just one unstable eigenvalue $s$ in the right-half-plane (excluding $s^{\ast }$ ). In fact, this was found to be true for all of the type ii and type iia solutions shown in figure 3(h). Table 1 illustrates the convergence of the unstable eigenvalue for a sample type ii case as the discretisation level $N$ and the truncation length of the computational domain $L$ are varied (similar results are found for a typical iia case). Evidently convergence is achieved with a rather large number of discretisation points. The spectrum for this case is shown in figure 6(a). The unstable eigenfunction is shown in figure 6(b). The large $x$ behaviour of the eigenfunction is controlled by the balance $(1/6)g_{k}^{\prime \prime \prime }\sim sg_{k}$ in (3.4). Therefore the decay is exponential with $g_{k}\sim \text{exp}(\unicode[STIX]{x1D707}x)$ as $|x|\rightarrow \infty$ , where $\unicode[STIX]{x1D707}=(6s)^{1/3}+2\text{i}k\unicode[STIX]{x03C0}/3$ for $k=0,1,2$ . For the unstable eigenvalue given in table 1, we obtain the dominant decay rate $\unicode[STIX]{x1D707}=0.73-0.41\text{i}$ for $x\rightarrow -\infty$ and $\unicode[STIX]{x1D707}=-0.01+0.84\text{i}$ as $x\rightarrow \infty$ . The latter value indicates a very slow decay downstream requiring a large computational domain. This explains the large value of $L$ required to achieve numerical convergence (see table 1).

Table 1.

Type ii solution at $a=-0.1357$ , $b=0.3$ ( $\unicode[STIX]{x1D6FC}=226.1$ ). Convergence of the eigenvalues in the discrete ${\mathcal{M}}$ spectrum computed on a computational domain $[-L,L]$ with $N$ equally spaced points. One value of $s$ in the ( $s$ , $-s$ , $s^{\ast }$ , $-s^{\ast }$ ) quartet is shown for various $N$ and $L$ (and grid size $\text{d}x=2L/N$ ).

Figure 6.

Stability properties of the type ii solution at $a=-0.1357$ , $b=0.3$ ( $\unicode[STIX]{x1D6FC}=226.1$ ): (a) Spectrum of the ${\mathcal{M}}$ operator and (b) the real and imaginary parts of the unstable eigenfunction for $s=0.0036+0.0997\text{i}$ . In (b) the inset shows a close up of the eigenfunction near to the origin.

As noted above, in the case of a neutrally stable steady solution to the present problem, all of the eigenvalues in the spectrum are necessarily located on the imaginary axis. Linear instability remains a possibility, however. For example, it may be the case that there exist non-simple point-spectrum eigenvalues embedded in the continuous spectrum on the imaginary axis which may lead to algebraic growth in time (Holm et al. Reference Holm, Marsden, Ratiu and Weinstein1985). To proceed we adopt an alternative approach suggested by Camassa & Wu (Reference Camassa and Wu1991), Chardard et al. (Reference Chardard, Dias, Nyguyen and Vanden-Broeck2011) and others, which exploits the fact that (3.1) has a Hamiltonian structure. In fact we can write it as

where ${\mathcal{H}}$ is the Hamiltonian. Consequently a steady solution $\unicode[STIX]{x1D702}_{s}$ of (3.1) corresponds to a solution of the equation $\unicode[STIX]{x1D6FF}{\mathcal{H}}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D702}=0$ . With a general perturbation $\unicode[STIX]{x1D701}$ as defined above, i.e. as $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{s}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}(x,t)$ , the second variation of ${\mathcal{H}}$ at $\unicode[STIX]{x1D702}_{s}$ is given by

If the second variation $\unicode[STIX]{x1D6FF}^{2}{\mathcal{H}}$ is positive (or negative) definite then the steady solution $\unicode[STIX]{x1D702}_{s}$ is said to be formally stable. Formal stability implies linear stability, although the converse is not necessarily true (see Holm et al. Reference Holm, Marsden, Ratiu and Weinstein1985). Linear instability can therefore be ruled out if $\unicode[STIX]{x1D6FF}^{2}{\mathcal{H}}$ is definite.

Since ${\mathcal{K}}$ is a self-adjoint Schrödinger operator, computing its eigenvalue spectrum constitutes a singular Sturm–Louiville problem on the real line (e.g. Titchmarsh Reference Titchmarsh1962; Coddington & Levinson Reference Coddington and Levinson1955) and is a well-studied problem in quantum mechanics and inverse scattering theory (e.g. Landau & Lifshitz Reference Landau and Lifshitz1977; Drazin & Johnson Reference Drazin and Johnson1989). Restricting the admissible perturbations $\unicode[STIX]{x1D701}$ to the class of $L_{2}(-\infty ,\infty )$ square-integrable functions, we obtain the lower bound (Camassa & Wu Reference Camassa and Wu1991; Keeler Reference Keeler2017),

Since $\unicode[STIX]{x1D708}_{0}>0$ , we can conclude from (3.8) that $\unicode[STIX]{x1D6FF}^{2}{\mathcal{H}}$ is positive definite. This means type i solutions (for the trench and Gaussian dip), which are negative definite, are formally stable and hence linearly stable. For our other solution types, $\unicode[STIX]{x1D6FF}^{2}{\mathcal{H}}$ is difficult to calculate analytically and instead we used the numerical pseudo-spectral method to determine the ${\mathcal{K}}$ spectrum and, hence, to calculate $\unicode[STIX]{x1D708}_{0}$ . (It should be emphasised that if for a given solution it is found that $\unicode[STIX]{x1D708}_{0}<0$ , so that the solution is formally unstable, this does not imply that it is linearly unstable and further investigation is needed.) For the type i solutions we confirmed numerically that $\unicode[STIX]{x1D708}_{0}>0$ implying linear stability in agreement with the calculations of the ${\mathcal{M}}$ spectrum discussed above. For the type ii and iia solutions we found that $\unicode[STIX]{x1D708}_{0}<0$ ; as we have already discussed, these solutions are found to be unstable when considering the spectrum of the ${\mathcal{M}}$ operator.

In the next subsection we perform time-dependent simulations of the unsteady fKdV equation with the goals of confirming the stability properties of the type i and ii solutions, and of following the development of the free surface in each case to large time.

3.2 Time-dependent simulations

We carried out time-dependent simulations of (1.2) on a truncated spatial domain $[-L,L]$ , for a specified value of $L$ , assuming periodic boundary conditions and using a method similar to that adopted by Chardard et al. (Reference Chardard, Dias, Nyguyen and Vanden-Broeck2011) and Wade et al. (Reference Wade, Binder, Mattner and Denier2014). The method computes spatial derivatives spectrally using a discrete Fourier transform and integrates forward in time using the fourth-order Runge–Kutta algorithm. Further details of the method can be found in (Trefethen Reference Trefethen2000). In numerical practice, we may obtain solutions that satisfy the boundary conditions in the far field by taking $L$ sufficiently large. In practice we found that taking $\unicode[STIX]{x0394}t\sim 10^{-2}$ and $512$ Fourier modes is sufficient to ensure the accuracy of the computations.

Since we assume periodicity in the computations, any disturbance leaving the domain will appear at the opposite end, and such behaviour is clearly unphysical. To prevent this in our computations we introduce a so-called sponge layer (see, for example, Boyd Reference Boyd2000; Alias, Grimshaw & Khusnutdinova Reference Alias, Grimshaw and Khusnutdinova2013; Grimshaw & Maleewong Reference Grimshaw and Maleewong2016) and in fact solve the equation

where $A(x)=D\{1+\tanh (b_{0}[x-b_{1}])\}$ for positive parameters $D$ , $b_{0}$ and $b_{1}$ . If we choose $b_{1}$ to be close to $L$ , then $A$ is non-zero over a narrow region at the left-hand end of the domain. This has the effect of preventing disturbances leaving the upstream end from entering the downstream end (due to the assumed periodicity) and vice versa.

Multiplying (1.2) by $\unicode[STIX]{x1D702}$ and integrating, after making use of (1.3) we obtain

where $P=\int _{-\infty }^{\infty }\unicode[STIX]{x1D702}^{2}\,\text{d}x$ is the momentum (see Camassa & Wu (Reference Camassa and Wu1991, Appendix C)). The expression $W$ is termed the wave resistance coefficient by Camassa & Wu (Reference Camassa and Wu1991) and is interpreted physically as a measure of the power being supplied by the forcing (see also Wu Reference Wu1987). In the absence of forcing, $P$ is a conserved quantity. In general $P=P(t)$ and so $W\neq 0$ ; it follows from the integral on the right-hand side of (3.10) that for a topography which is even in $x$ , as will be assumed in this paper, $\unicode[STIX]{x1D702}(x,t)$ will have a left–right asymmetry about $x=0$ . In the presence of forcing it is straightforward to show that the excess mass $Q=\int _{-\infty }^{\infty }\unicode[STIX]{x1D702}\,\text{d}x$ is a conserved quantity. The Hamiltonian defined in (3.6) provides a second conserved quantity, corresponding to the energy, and both of these can be used as a check on the fidelity of numerical time-dependent simulations of (1.2). A further check can be made by computing $W$ using either of the two alternative expressions given in (3.10) and confirming that the same result is obtained in either case. In the numerical results to be discussed in the next section, both $Q$ and ${\mathcal{H}}$ were performed in the absence of the sponge layer and it was confirmed that they are conserved.

In the results to be presented below a time-dependent calculation is started using the initial condition $\unicode[STIX]{x1D702}(x,0)=\unicode[STIX]{x1D702}_{s}+\unicode[STIX]{x1D716}\hat{\unicode[STIX]{x1D701}}$ where $\unicode[STIX]{x1D702}_{s}$ is a chosen steady solution, $\hat{\unicode[STIX]{x1D701}}$ is a perturbation and $\unicode[STIX]{x1D716}$ is taken to be a small number. In all of the calculations to be presented below we took $\unicode[STIX]{x1D716}=0.025$ . If the numerically computed eigenspectrum of the operator ${\mathcal{M}}$ revealed a positive largest growth rate, the initial perturbation $\hat{\unicode[STIX]{x1D701}}=(g_{k}+g_{k}^{\ast })/2$ , with $g_{k}$ the corresponding eigenfunction, was used. Otherwise if the numerically computed spectrum was found to contain only purely imaginary eigenvalues, the initial perturbation $\hat{\unicode[STIX]{x1D701}}=\unicode[STIX]{x1D702}_{s}$ was used.

We note that for a given initial condition and a particular choice of the topographic forcing $f(x)$ , the solution to the fKdV equation (1.2) is known to be unique (Wu Reference Wu1987). Also, we remark that when presenting unsteady simulations of the fKdV equation for topographic forcing similar to that considered here Wu (Reference Wu1987) and Camassa & Wu (Reference Camassa and Wu1991) both solved the so-called regularised form of the fKdV equation (see Benjamin, Bona & Mahony Reference Benjamin, Bona and Mahony1972), whose equivalence is assured only in the limit of long waves.

We focus on the case of Gaussian topography with $f(x)=a\exp (-b^{2}x^{2})$ . Our first goal is to confirm the stability results of figure 6. Figure 7 shows the results of a numerical simulation starting from a small perturbation of the unique type i steady solution $\unicode[STIX]{x1D702}_{s}$ at $\unicode[STIX]{x1D6FC}=141$ ( $a=-0.085$ , $b=0.3$ ). The free-surface profile $\unicode[STIX]{x1D702}(x,t)$ is shown at different times over a portion of the computational domain $[-L,L]$ . We note that the solution in the left side of the domain, which is not shown, remains uniformly zero over the duration of the simulation and that no disturbances are observed to propagate upstream. Also shown in the figure are the time signals of the momentum $P$ and the wave resistance coefficient $W$ . The initial condition $\unicode[STIX]{x1D702}(x,0)=1.025\unicode[STIX]{x1D702}_{s}$ has momentum $P(0)=5.53\times 10^{-3}$ which is in excess of the momentum of the steady state type i solution momentum $P_{s}=5.26\times 10^{-3}$ . The simulation shows that ultimately the system returns to the steady state $\unicode[STIX]{x1D702}_{s}$ as time increases. The wave resistance coefficient $W$ approaches zero at large time, which is a necessary condition for the steady state to be recovered. Also $W$ is negative during the simulation and so, according to (3.10), the momentum $P$ must decrease monotonically. In particular we see that $P\rightarrow P_{s}$ as $t\rightarrow \infty$ . A $\log$ - $\log$ plot of the $W$ signal indicates that $W\sim O(t^{-2})$ for large $t$ , and so $P$ is expected to approach $P_{s}$ algebraically in time; this explains the very long integration time required to recover the steady state.

We emphasise again that momentum cannot be released from the free surface as long as it is symmetric about $x=0$ since $W$ is zero in this case. Since the initial condition is symmetric, and moreover has a greater momentum than the steady state $\unicode[STIX]{x1D702}_{s}$ , it follows that an asymmetric free-surface disturbance must develop in $t>0$ if the steady state is to be recovered in the long term. As can be seen in figure 7(a–d), small amplitude waves appear on the leeward side of the free surface. These waves propagate downstream carrying momentum with them. Eventually these waves evanesce and the steady state is recovered. Close inspection of these figures reveals that the free surface, shown with a solid line, overlaps with the steady state $\unicode[STIX]{x1D702}_{s}$ , shown with a broken line, over an increasingly wide region as time goes on.

Figure 7.

Numerical simulation of (3.9) with Gaussian forcing for $\unicode[STIX]{x1D6FC}=141$ , ( $a=-0.085,b=0.3$ ) with $L=200$ and $N=512$ . The initial condition is $\unicode[STIX]{x1D702}(0,t)=1.025\unicode[STIX]{x1D702}_{s}$ , where $\unicode[STIX]{x1D702}_{s}$ is the unique type i solution at this $\unicode[STIX]{x1D6FC}$ (see figure 3 h). (a–d) Show the solution over a portion of the domain at times $t=40$ , $200$ , $2000$ and $10\,000$ . (e,f) Show the momentum $P$ and $W$ respectively against time. The broken line in (e) is the momentum of the type i steady state, $P=5.26\times 10^{-3}$ .

Figure 8.

Numerical simulation of (3.9) with Gaussian forcing for $\unicode[STIX]{x1D6FC}=226.1$ , ( $a=-0.1357,b=0.3$ ) with $L=400$ and $N=2048$ (note not all of the domain is displayed). The initial condition is $\unicode[STIX]{x1D702}(0,t)=1.025\unicode[STIX]{x1D702}_{s}$ , where $\unicode[STIX]{x1D702}_{s}$ is the unique type ii solution at this $\unicode[STIX]{x1D6FC}$ (see 3 h). (a–d) Show the solution over a portion of the domain at times $t=100$ , $t=500$ , $t=1000$ , $1100$ and $1500$ . (f) Shows $W$ respectively against time. The dotted curve is the steady solution $\unicode[STIX]{x1D702}_{s}$ .

Figure 8 shows the results of a numerical simulation starting from a small perturbation of the type ii steady solutions $\unicode[STIX]{x1D702}_{s}$ at $\unicode[STIX]{x1D6FC}=226.1$ ( $a=-0.1357$ , $b=0.3$ ) respectively. Very similar results are found for a simulation started from a small perturbation of a type iia steady solution. The free-surface profile $\unicode[STIX]{x1D702}(x,t)$ is shown at different times over a portion of the computational domain $[-L,L]$ . The calculated unstable eigenmode is used as an initial condition for the linearised time-dependent equation (3.2) and excellent agreement is found with the calculated unstable eigenvalue (Keeler Reference Keeler2017). The agreement for the nonlinear time-dependent equation is also good for small times when the initial condition is given by $\unicode[STIX]{x1D702}(x,0)=\unicode[STIX]{x1D702}_{s}+1/2\unicode[STIX]{x1D716}(g_{k}+g_{k}^{\ast })$ , with $g$ the unstable eigenmode. The initial condition for the calculations in figure 8 is $\unicode[STIX]{x1D702}(x,0)=1.025\unicode[STIX]{x1D702}_{s}$ and the different panels show the free-surface profile at different times together with the steady solution, which is shown as a dashed curve. These snapshots show waves travelling downstream in an irregular manner. Additionally, the solution in $x<0$ does not remain uniformly zero over the duration of the simulation and after a finite interval soliton-like disturbances emerge and propagate upstream. To conserve $Q$ , the excess mass created by these soliton-like disturbances is balanced by a non- uniform depression region which appears directly downstream of the topography. The time signal of the wave resistance coefficient $W$ is also shown in the figure and is seen to oscillate with increasing amplitude as time increases. We conclude the solution does not return to the steady state as $t$ increases.

We also carried out simulations with the initial condition $\unicode[STIX]{x1D702}(x,0)=0$ . Although a stable steady solution exists in this case, this is not approached at large time. Instead unsteady behaviour is observed whose features are in line with those described by Wu (Reference Wu1987), Camassa & Wu (Reference Camassa and Wu1991), Grimshaw & Smyth (Reference Grimshaw and Smyth1986) (see in particular § 6 and figure 11) and Grimshaw, Zhang & Chow (Reference Grimshaw, Zhang and Chow2007).