2. The KdV–Burgers equation

The KdV equation for water waves is very well known. In standard notation for a surface displacement $\eta$ above an undisturbed depth $h$ it is the left-hand side of (2.1). The forcing/friction term is derived in MG24 and reduced here to the shallow-water limit in a unidirectional setting as described above and in more detail below

(2.1)\begin{gather} \eta _t + C \eta_x + \frac{3C}{2h}\eta \eta_x + \frac{Ch^2}{6}\eta_{xxx} = F(\eta), \quad C = (gh)^{1/2}, \end{gather}

(2.2)\begin{gather}F(\eta )= \frac{h}{2C} \left\{ \left[-\frac{\rho_{a}}{\rho_w}\beta W_{r}^2 + 4\kappa \frac{C}{h } \right]\eta_{xx} - \frac{1}{h} \frac{\rho_a}{\rho_w} c_d U_{a}^2 (|\eta| \eta)_{xx} - \frac{g}{h^2} C_D |\eta | \eta \right\} . \end{gather}

Here, $g$ is the acceleration due to gravity, $\rho _a$ is the air density, $\rho _w$ is the water density, $\kappa$ is the kinematic viscosity, $c_d$ is the wind surface drag coefficient, $W_r$ is a reference wind velocity introduced by Miles (Reference Miles1957), $U_a$ is a scaling velocity used in MG24 for the surface wave stress and $C_D$ is the bottom drag coefficient. The first term has two parts. The first part describes an unstable KdV–Burgers equation and the second part describes a stable KdV–Burgers equation. The second term describes a nonlinear unstable KdV–Burgers equation. The third term is a decay term modelling turbulent bottom stress. Hereafter, we refer to (2.1) as the KdV–Burgers equation, although it is a much modified version.

The first part of the first term in (2.2) arises from a representation of air pressure $P_a$ at the air–water interface, due to Miles (Reference Miles1957) and used in his pioneering paper on wind waves

(2.3)\begin{equation} \frac{P_{a}}{\rho_a} = (\alpha |k|\eta+ \beta\eta_x ) W_{r}^2 . \end{equation}

Here, $k, \omega$ are the wavenumber and frequency for a sinusoidal wave, $\eta \propto \exp {({\rm i} kx - {\rm i} \omega t )}$ in a linear reduction and $\alpha, \beta$ are non-dimensional real parameters, found by solving the linearised air flow equations. Although $\alpha, \beta$ depend on the depth $h$ this emerges only for very small depths too small to be of practical interest, so here we set $\alpha = 0.872, \beta = 1.745$ and $W_r =0.9\ {\rm m}\ {\rm s}^{-1}$, as in MG24. Further, since only the term in $\beta$ is related to the critical level instability mechanism of Miles (Reference Miles1957), the term in $\alpha$ is omitted to derive the first term in (2.2); see Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010). Finally, we put ${\rm i} k \to \partial /\partial x$. It is pertinent to note here that Miles (Reference Miles1957) assumed that $k> 0$, as we did in MG24 and was done in many works, which is allowed in the application to linear, sinusoidal waves. When $k< 0$, then $k\eta$ in Miles (Reference Miles1957) has been replaced by $|k|\eta$ in (2.3), noting that for real-valued functions $\eta$, both $|k|, {\rm i} k$ are valid Fourier transforms, unlike just $k$. Also, the solutions of the modal equation for the linearised air flow depend on $|k|$ rather than just $k$. The second part of the first term in (2.2) is due to a laminar frictional boundary layer at the water surface (see Dutykh & Dias (Reference Dutykh and Dias2007a,Reference Dutykh and Diasb), Dutykh (Reference Dutykh2009) and Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010)); the kinematic viscosity $\kappa = 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$.

The second term in (2.2) arises from a representation of wave stress in the turbulent near-surface air flow. The wave stress is $\rho _a c_d |W_a|W_a$, where $W_a$ is a near-surface wind speed; see Tang & Grimshaw (Reference Tang and Grimshaw1995) and Tang et al. (Reference Tang, Grimshaw, Sanderson and Holland1996) for instance. In MG24 we set $W_a = u_a$, the horizontal air flow velocity at the free surface found from the linearised air flow equations, and given by $\rho _a c u_a = P_{a}$, where $P_a$ is given by (2.3) and $c = \omega /k$ is the phase velocity of the underlying carrier wave with frequency $\omega$ and wavenumber $k$. However, in MG24 the resulting expression is too sensitive to the depth in the shallow-water limit to be useful here. Instead, we retain the expression $u_a = |k|\eta U_a$ from MG24 but here fix the scaling velocity $U_a = 5\ {\rm m} {\rm s}^{-1}$ as a typical value for the parameters we use here. Further, the asymptotic KdV reduction has lead to an expression $|\eta |\eta$ in (2.2) which is difficult to handle, both numerically in a Fourier spectral code and analytically in the modulation theory which follows. Hence, we take advantage of the fact that $c_d$ increases with the wind speed; see Tang & Grimshaw (Reference Tang and Grimshaw1995) and Tang et al. (Reference Tang, Grimshaw, Sanderson and Holland1996) for instance. The outcome is that we replace $c_d$ with $c_{d}|\eta | /A_0$, where $A_0$ is a scaling wave amplitude and $c_d$ is now the constant value at $\eta = A_0$. Hence, in (2.2), the term $c_d U_{a}^2 (|\eta | \eta )_{xx}$ is replaced with $c_d U_{a}^2 ( \eta ^3)_{xx}/A_0$. We set $c_d = 0.0023, A_0 = 0.5\ {\rm m}$ as typical values.

The third term is a conventional representation of decay due to turbulent bottom stress $C_D |U| U$ (see Tang & Grimshaw (Reference Tang and Grimshaw1995) and Tang et al. (Reference Tang, Grimshaw, Sanderson and Holland1996) for instance) where $U$ is the bottom velocity. In the full shallow-water limit $U = C \eta /h$, but in the linearised periodic wave formulation of MG24, $U = (g \eta /c) \cosh {(kh)} = c k\eta /\sinh {(kh)}$, where $c = \omega /k$ is the linear phase speed for a carrier wave of frequency $\omega$ and wavenumber $k$. In the shallow-water limit $kh \to 0, c \to C$ and we recover $U =C \eta /h$ but here we retain the dependence on $kh$, even though $kh$ is specified for (2.1) only through (2.2) and the initial conditions. In effect, $C_D$ acquires a factor $f_r = c^2 (kh)^2 / (C^2 \sinh ^2{(kh)})$. This depends only on $kh$, $f_r = (kh) \tanh {(kh)}/\sinh ^2{(kh)} = 2kh/\sinh {(2kh)}$ and $f_r \to 1$ as $kh \to 0$ but becomes quite small as $kh$ increases. The device used in the second term to replace $|\eta |\eta$ with $\eta ^3$ is not available here as $C_D$ does not depend directly on the wind speed. Instead, we use the strategy in MG24 and assume that as $|\eta |$ decreases the nonlinear term $|\eta |\eta$ can be replaced by the linear term $A_0 \eta$. As the KdV equation is a small-amplitude approximation, we replace $C_D | \eta | \eta$ in (2.2) with $C_D A_0 \eta$ where the drag coefficient remains at $C_D = 0.0015$ and we again put $A_0 = 0.5\ {\rm m}$, the same value as in the second term. In effect the quadratic bottom stress is replaced by a linear term. If this is modelled by a laminar bottom boundary layer we find from MG24 (see also Dutykh & Dias (Reference Dutykh and Dias2007a,Reference Dutykh and Diasb)) that $C_D A_0$ is then replaced with $kh [\kappa /2 \omega ]^{1/2}$. This gives a connection between $C_D$ and the kinematic viscosity $\kappa$ but we caution that the connection depends on $k$ and $h$. Nevertheless, for the plausible values $\omega =1\ {\rm s}^{-1}$, $kh =1$ and $A =0.5\ {\rm m}$, $\kappa = 10^{-6} \ {\rm m}^2 \ {\rm s}^{-1}$ yields $C_D = 0.0014$ quite close to the accepted value of $C_D = 0.0015$.

Next, we express (2.1) in a reference frame moving with speed $C$, so that $x$ is replaced by $x^{\prime } = x - C t$; henceforth, the superscript on $x^{\prime }$ is omitted. Then, we put (2.1) into canonical non-dimensional form using amplitude, length and time scales $E_S, L_S, T_S$; that is, $\eta _d = E_S \eta _n , x_d = L_S\, x_n, t_d = T_S\, t_n$, where the subscript $d, n$ denotes the dimensional, non-dimensional variables, respectively. We set $E_S = 2h(3K^2)^{-1} m, T_S = 6K^3h C^{-1} s, L_S = Khm$, where $K$ is a non-dimensional scaling parameter chosen to ensure that the waves have small amplitudes and are long compared with the depth. The outcome is, replacing ((2.1) and (2.2)), omitting the subscript $n$ and incorporating the changes to $c_d, C_D$ specified above

(2.4)\begin{gather} \eta _t + 6\eta \eta_x + \eta_{xxx} = F(\eta), \end{gather}

(2.5)\begin{gather}F(\eta ) ={-}\mu_1 \eta_{xx} - \mu_2 (\eta^3)_{xx} -\mu_3 \eta , \end{gather}

(2.6)\begin{gather}\mu_{1,2,3} = \left\{K \left[3\frac{\rho_{a}}{\rho_w} \beta \frac{W_{r}^2}{C^2} - \frac{12 \kappa}{hC} \right] , \frac{2 }{K}\frac{\rho_a}{\rho_w} c_d \frac{hU_{a}^2 }{A_0 C^2 }, 2~K f_r C_D \frac{A_0 }{h} \right\}. \end{gather}

Here, $\mu _{1, 2, 3}$ are dimensionless constants, adjusted slightly as described above from the equivalent values in (2.2); $\mu _1$ is the sum of a term due to critical level instability and a term due to laminar friction at the water surface, and will be positive or negative indicating wave growth or decay, depending on which term is the larger (see Miles (Reference Miles1957) and Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010) for a discussion of this in the deep water context); $\mu _1$ varies with water depth $h$ due to the dependence of $\beta$ on $h$ and the factor $W_{r}^2 /C^2$ in (2.6), see figure 1(a), where the plot of $\beta$ is based on the expression in MG24 for a logarithmic wind shear profile and a $5$-second carrier wave. In general, there is a critical depth $h_{cr}$ such that $\mu _1 \gtrless 0$ according as $h \gtrless h_{cr}$. For our parameter settings, $h_{cr} = 2.5 \ {\rm m}$ is quite shallow. The term with coefficient $\mu _2$ describes wave growth due to a representation of wave stress in the air at the water surface induced by a turbulent wind. In our derivation, $\mu _2$ is independent of $h$ as the term in (2.6) $hU_{a}^2 /A_0 C^2 =U_{a}^2 /gA_0$ does not depend on $h$. The term with coefficient $\mu _3$ describes wave decay due to a representation of a turbulent bottom boundary layer. The term $\mu _3$ increases rapidly when the depth is very small and decays to zero as the depth increases, see figure 1(b). Based on these considerations we restrict attention to $h > 5\ {\rm m}$ when $C > 7\ {\rm m}\ {\rm s}^{-1}$, as the variability at smaller depths is unrealistic in practice, and then each $\mu _{i}$ is a small parameter, so facilitating the asymptotic analysis of § 3. In the numerical simulations of § 4, we choose a specific $h = 10\ {\rm m}$ when $C = 9.9 \ {\rm m}\ {\rm s}^{-1}$. Although it was shown in MG24 that $\beta$ depends on depth the variation is slight for $h \ge 10 \,m$ and so here we put $\beta = 1.745$, constant, the value found by MG24 in water of sufficient depth for a logarithmic wind shear profile and a $5$-second carrier wave.

Figure 1. Plots of $\mu _{1}$ (a) and $\mu _3$ (b) vs depth $h$ in shallow water for a $5$-second carrier wave.

The KdV–Burgers equation (2.4) has two conservation laws, for mass and wave action. For modulated periodic waves, there is an additional conservation law, conservation of waves, see § 3. We assume that $\eta$ is spatially periodic in a domain $[-L, L]$ and then the mass conservation law comes directly from (2.4)

(2.7a,b)\begin{equation} \boldsymbol{M} = \int_{{-}L}^{L} \eta \, {{\rm d}\kern0.07em x}, \quad \boldsymbol{M}_t = \int_{{-}L}^{L} F(\eta )\, {{\rm d}\kern0.07em x} ={-} \mu_3 \boldsymbol{M} . \end{equation}

Note that the bottom stress term with coefficient $\mu _3 >0$ causes $\boldsymbol {M} \to 0$, an artefact of the one-way KdV asymptotic approximation, which requires that the horizontal velocity field is proportional to $\eta$. The energy equation obtained from (2.4) is

(2.8)\begin{equation} \left( \frac{\eta ^2 }{2} \right)_t + (2\eta^3 )_x + \left(\eta \eta_{xx} - \frac{\eta_x^2}{2}\right)_x = \eta F(\eta) . \end{equation}

Wave action, the same as wave energy here, is then defined by

(2.9)\begin{equation} \boldsymbol{E} = \int_{{-}L}^{L} \frac{1}{2} \eta^2\, {{\rm d}\kern0.07em x} . \end{equation}

In the absence of forcing/friction $\boldsymbol {E}$ is conserved. Otherwise

(2.10)\begin{equation} \boldsymbol{E}_t = \int_{{-}L}^{L} \eta F(\eta)\, {{\rm d}\kern0.07em x} = \int_{{-}L}^{L} \left[\mu_1\eta_{x}^2 + 3\mu_2 \eta_{x}^2 \eta^2 \right]\, {{\rm d}\kern0.07em x} - 2\mu_3 \boldsymbol{E} . \end{equation}

The first term is a growth or decay term according as $\mu _1 \gtrless 0$, the second term is a growth term as $\mu _2 > 0$ and the third term is a decay term as $\mu _3 > 0$. There are analogous expressions for solitary wave solutions where $\eta \to 0$ as $x \to \pm \infty$ and $[-L, L ]$ is replaced by $(-\infty, \infty )$.

The stable KdV–Burgers equation with $\mu _1<0, \mu _2 =\mu _3 =0,$ was first proposed by Su & Gardner (Reference Su and Gardner1969) as an extension of the KdV equation which combined nonlinearity, dispersion and canonical dissipation. Independently, in a study of waves in a liquid-filled elastic tube, Johnson (Reference Johnson1970) proposed the same equation, and examined its steady-state solution as a potential model for a bore, either undular or monotonic. In the wind-wave context, Zdyrski & Feddersen (Reference Zdyrski and Feddersen2021) and Manna & Latifi (Reference Manna and Latifi2023) have recently proposed the unstable ($\mu <0$) KdV–Burgers equation as a model for wind waves in shallow water. Both used the Jeffreys pressure asymmetry theory as the driving mechanism rather than the Miles critical level instability theory. Zdyrski & Feddersen (Reference Zdyrski and Feddersen2021) used numerical simulation to examine wave shape, while Manna & Latifi (Reference Manna and Latifi2023) used an asymptotic analysis, similar to that we present in § 3 to show that a solitary wave could grow in amplitude and develop a singularity.

3. Modulated solitary wave train

3.1. Cnoidal waves

When unforced and in the absence of friction the KdV–Burgers equation (2.4) reduces to the KdV equation when it has a well-known periodic wave solution, the cnoidal wave, described in many places, including the original 1895 KdV paper. In the presence of perturbation effects such as variable depth, or as here the forcing/friction term $F(\eta )$ (2.5), this becomes a modulated cnoidal wave, introduced by Gurevich & PitaevskiI (Reference Gurevich and PitaevskiI1974) and Whitham (Reference Whitham1974). It has been heavily used since, see for instance Grimshaw & Smyth (Reference Grimshaw and Smyth1986), Smyth (Reference Smyth1987, Reference Smyth1988), Myint & Grimshaw (Reference Myint and Grimshaw1995), Grimshaw (Reference Grimshaw2010), El, Grimshaw & Tiong (Reference El, Grimshaw and Tiong2012) and Grimshaw & Yuan (Reference Grimshaw and Yuan2016, Reference Grimshaw and Yuan2018) in a fluid mechanics context. Here, we follow a similar approach as in these references, described briefly here in the following text, and based on the premise that $F(\eta )$ is a small perturbation.

For a modulated periodic wave with wavenumber $k$ and frequency $\omega =kV$, where $V$ is the phase speed, we introduce a phase $\theta (x, t)$ such that

(3.1a,b)\begin{equation} \theta_x = k , \quad \theta_t ={-} kV . \end{equation}

The equation for conservation of waves comes from eliminating $\theta$ from (3.1a,b)

(3.2)\begin{equation} k_t + (kV)_x =0 . \end{equation}

We seek a modulated wave of the form $\eta =\eta (\theta, x,t)$ which is periodic in $\theta$ with period $2P$ and where the dependence on $x,t$ is slowly varying and expressed through the available wave parameters. Without loss of generality we can put $P = {\rm \pi}$. The spatial period for a single wave is $2P/k$ and the full spatial domain is $[-L, L]$ containing $N$ waves. Next, we define a wave average for a single wave

(3.3)\begin{equation} \langle [\cdots ]\rangle = \int_{{-}P}^{P} [\cdots] \frac{{\rm d}\theta}{2P} . \end{equation}

We suppose that

(3.4)\begin{equation} \eta = \hat{\eta}(\theta, x, t) + D (x, t), \quad \langle \hat{\eta}\rangle =0 , \end{equation}

where $D( x,t )= \langle \eta \rangle$ is the mean level and is a slowly varying function of $x, t$. The averaged wave action density, just wave energy in this context, is

(3.5)\begin{equation} E = \left\langle \frac{\hat{\eta}^2 }{2}\right\rangle , \end{equation}

and is likewise a slowly varying function of $x, t$.

On expressing the fast dependence on $\theta$ and the slow dependence on $x, t$, the KdV–Burgers equation (2.4) becomes at leading order

(3.6)\begin{equation} -k(V- 6D)\hat{\eta}_{\theta} + 3k(\hat{\eta}^2)_{\theta } +k^3 \hat{\eta}_{\theta \theta \theta } = \cdots , \end{equation}

where the omitted terms $[\cdots ]$ contain the slowly varying $x, t$ derivatives and the forcing/friction terms. One integration with respect to $\theta$ yields

(3.7)\begin{equation} {-}(V- 6D)\hat{\eta} + 3\hat{\eta}^2 +k^2 \hat{\eta}_{\theta \theta } = 6E + \cdots . \end{equation}

The wave equation (3.7) has the modulated cnoidal wave solution given by

(3.8)\begin{gather} \hat{\eta} = A\{b + cn^2(\gamma \theta; m)\} , \quad b ={-}\langle cn^2 \rangle = \frac{1-m}{m}-\frac{E(m)}{mK(m)}, \end{gather}

(3.9a,b)\begin{gather}A = 2m \varGamma^2, \quad V = 6D + 4\varGamma^2 \left\{2-m -\frac{3E(m)}{K(m)}\right\},\quad \varGamma = \gamma k . \end{gather}

Here, we have used the description in Grimshaw & Smyth (Reference Grimshaw and Smyth1986), El (Reference El2007), Grimshaw (Reference Grimshaw2010) and Grimshaw & Yuan (Reference Grimshaw and Yuan2016). In ((3.8) and (3.9a,b)) $cn(x; m)$ is the Jacobi elliptic function of modulus $m$ ($0< m<1$) and $K(m), E(m)$ are the elliptic integrals of the first and second kind, respectively, defined by

(3.10)\begin{gather} cn(x;m) = \cos{(\phi)} , \quad x = \int_{0}^{\phi } \frac{{\rm d}\phi^{\prime }}{(1-m\sin^2 {\phi^{\prime }})^{1/2}} , \quad 0 \le \phi \le \frac{\rm \pi}{2} , \end{gather}

(3.11a,b)\begin{gather}K(m) = \int_{0}^{{\rm \pi} /2 } \frac{{\rm d} \phi }{(1-m\sin^2 {\phi })^{1/2}} , \quad E(m) = \int_{0}^{{\rm \pi} /2 }\, (1-m\sin^2 {\phi })^{1/2}\, {\rm d}\phi . \end{gather}

The cnoidal wave (3.8) is periodic in $\theta$ with a period $2P = 2K(m)/\gamma$, which for fixed $P = {\rm \pi}$ defines $\gamma = K(m)/{\rm \pi}$. The corresponding spatial period is $2{\rm \pi} /k =2K(m)/\varGamma$. Here, $V$ is the wave phase speed, $A$ is the wave amplitude and the maximum, minimum values of $\hat {\eta }$ are $A_M = A(b + 1), A_m = Ab$. The wave action density is given by

(3.12) \begin{equation} \left.\begin{array}{c@{}} \displaystyle E = \left\langle \dfrac{\hat{\eta}^2 }{2} \right\rangle = \dfrac{A^2}{2}(C_4 - b^2), \quad C_4 = \langle cn^4\rangle , \\ \displaystyle C_4 = \dfrac{1}{3m^2 K(m) } \{3m^2 K(m) - 5mK(m) + 4mE(m) +2K(m)-2E(m) \} , \end{array}\right\} \end{equation}

where $C_4 , b$ are functions of the modulus $m$ only.

As the modulus $m \to 0$, $cn(x; m) \to \cos {(x)}$, $\gamma \to 1/2$, $b \to -1/2$ and then the cnoidal wave (3.8) collapses to a linear sinusoidal wave

(3.13a–c)\begin{equation} \hat{\eta} = \frac{A}{2}\,\cos{( \theta)}, \quad V = 6D - k^2 , \quad E = \frac{ A^2}{8}. \end{equation}

In this limit $A \to 0$, but $A/m = 2\varGamma ^2$ stays finite.

As the modulus $m \to 1$, $cn(x; m) \to \hbox {sech}(x)$, $K(m), E(m) \to \infty, 1$ and $b \to 0$. The cnoidal wave (3.8) becomes a solitary wave train, defined on a periodic lattice, and riding on a background level $D$

(3.14a–c)\begin{equation} \hat{\eta}= A \,\hbox{sech}^2 (\gamma \theta ), \quad A = 2\varGamma^2 , \quad V = 6D + 4\varGamma^2 . \end{equation}

In this limit, $\gamma \to \infty, k \to 0$ but $\varGamma =\gamma k$ stays finite, and $\gamma \theta$ can be replaced by $\varGamma (x- V t )$ in integral expressions; $C_4 , E \to 0$ but $\mathcal {E} = E/k$ remains finite and is given by

(3.15)\begin{equation} \mathcal{E} = \frac{2A^2}{3\varGamma } = \frac{8\varGamma^3 }{3}= \frac{(2A)^{3/2}}{3} . \end{equation}

For pure solitary waves, where it is required that $\eta \to 0$ at infinity, formally $D$ must be zero. However, under perturbation a trailing shelf is generated; see Grimshaw (Reference Grimshaw1979) and Grimshaw & Mitsudera (Reference Grimshaw and Mitsudera1993) for instance. It is important here to distinguish between a pure solitary wave and a solitary wave train. The latter is essentially a useful approximation to a modulated cnoidal wave train, often used to describe undular bores; see Grimshaw & Smyth (Reference Grimshaw and Smyth1986), Smyth (Reference Smyth1987), El (Reference El2007) and Grimshaw (Reference Grimshaw2010); perturbative effects such as friction and variable depth have been examined by Smyth (Reference Smyth1988), Myint & Grimshaw (Reference Myint and Grimshaw1995), El et al. (Reference El, Grimshaw and Tiong2012) and Grimshaw & Yuan (Reference Grimshaw and Yuan2016, Reference Grimshaw and Yuan2018), amongst many others.

3.2. Modulation analysis

The parameters $[A, \varGamma, k, m, V, D]$ are reduced to three by the relations in (3.9a,b). They are slowly varying functions of $(x, t )$ found from a multi-scale asymptotic expansion, or more directly using the conservation laws for mass (2.4) and wave action (2.8) combined with the equation for conservation of waves (3.2). The wave average of the conservation of mass equation (2.4) yields

(3.16)\begin{equation} D_t + 6\mathcal{M}_x = \langle F(\eta ) \rangle ={-}\mu_3 D , \quad \mathcal{M} = \frac{\langle \eta^2 \rangle }{2 } =\frac{D^2}{2} + E. \end{equation}

A slowly varying term $D_{xxx}$ on the left-hand side has been omitted because it is much smaller than the retained terms due to the slow variation in $x$. Likewise, two small terms $-\mu _1 D_{xx}$ and $-\mu _2 \langle \eta ^3 \rangle _{xx}$ on the right-hand side of (3.16) have been omitted. Further, when $D^2 \gg 2E$, (3.16) reduces to a Hopf equation with implicit solution

(3.17) \begin{equation} \left.\begin{array}{c@{}} \hat{D}_{\tau } +6\hat{D}\hat{D}_x = 0, \quad D = \hat{D}\exp{(-\mu_3 t )}, \quad \tau = \dfrac{1- \exp{(-\mu_3 t )}}{\mu_3 } , \\ \hat{D} = D_{0} (x - 6\hat{D}\tau ) , \quad D(x, t =0) = D_0 (x) . \end{array}\right\} \end{equation}

For $0 < t < \infty$, $0 < \tau < 1/\mu _3$. If $D_{0}(x) = D_0$ is a constant the solution is the constant $\hat {D} = D_0$; $D =D(t)$ is independent of $x$ and decays exponentially to zero as $t \to \infty$. Otherwise, an infinite slope is reached if $1 + 6D_{0x} \tau =0$ which needs $D_{0x} < 0$ somewhere. For a sufficiently small initial condition such that $\mu _3 + 6D_{0x} >0$ for all $x$ this does not occur since $\mu _3 \tau < 1$. However, otherwise it is inevitable and can only be prevented by the restoration of a dispersive term such as $D_{xxx}$ to (3.16).

One more modulation equation is needed and this comes from taking the wave average of the conservation of energy equation (2.8)

(3.18) \begin{equation} \left.\begin{array}{c@{}} \mathcal{M}_t +\mathcal{Q}_x = \langle \eta F(\eta)\rangle, \quad \mathcal{Q} = \left\langle 2\eta^3 - \dfrac{3\hat{\eta}_{x}^2 }{2} \right\rangle , \\ \mathcal{Q} = 2D^3 + 12DE + \left\langle 2\hat{\eta}^3 - \dfrac{3\hat{\eta}_{x}^2 }{2} \right\rangle . \end{array}\right\} \end{equation}

Next, $D_t$ is eliminated from (3.18) using (3.16) and the wave equation (3.7) to generate an equation for $E$

(3.19)\begin{gather} E_t +(VE + \mathcal{U})_x +6ED_x = \langle \hat{\eta} \hat{F} \rangle , \quad \mathcal{U} = \left\langle \frac{\hat{\eta}^3 }{2} - \hat{\eta}_{x }^2 \right\rangle ,\quad \hat{F} = F - \langle F\rangle , \end{gather}

(3.20)\begin{gather}\text{or if } E = k\mathcal{E} , \quad k[\mathcal{E}_t + V\mathcal{E}_x + 6\mathcal{E}D_x ] +\mathcal{U}_x = \langle \hat{\eta} \hat{F} \rangle , \end{gather}

(3.21)\begin{gather}\langle \hat{\eta} \hat{F} \rangle = \mu_1 \langle \hat{\eta}_{x}^2 \rangle + 3\mu_2 \langle \hat{\eta}_{x}^2 \eta^2\rangle - \mu_3 2E . \end{gather}

Here, $\mathcal {E}$ is the spatially averaged wave action density, and ((3.18) and (3.20)) agree with analogous derivations in the literature; for instance, see Grimshaw & Yuan (Reference Grimshaw and Yuan2016) for a detailed review of modulated cnoidal waves for a variable-coefficient KdV equation. The term $\mathcal {U}$ is needed to ensure that the modulated wave propagates with the nonlinear group velocity and not in general with the phase speed $V$

(3.22) \begin{align} \mathcal{U} \!=\! A^3 \left[- b^3 + \frac{3bC_4}{2} +\frac{C_6}{2} \right] - 4A^2 \varGamma^2 [(1\!-\!m)(C_2- C_4) +m (C_4- C_6) ], \quad C_n =\langle cn^n\rangle . \end{align}

When $m = 0$, $\mathcal {U} = -2k^2 E$ and as then $V =6D-k^2$, we obtain the linear group velocity of $6D-3k^2$ corresponding to a linear frequency of $6Dk -k^3$. However, when $m \to 1$ we find that $\mathcal {U} \to 0$.

The three ((3.2), (3.16), (3.19)) form a nonlinear hyperbolic system for the chosen three parameters from the set $[A, \varGamma, m, V, D]$, say $A, m, D$. In the absence of the forcing/friction term, that is $F(\eta )=0$, the elegant Whitham modulation theory exploits the underlying integrability of the KdV equation and yields solutions using Riemann invariants; see Whitham (Reference Whitham1974) and Gurevich & PitaevskiI (Reference Gurevich and PitaevskiI1974) and the comprehensive review by El (Reference El2007). Of particular interest is the similarity solution

(3.23) \begin{equation} \left.\begin{array}{c@{}} \displaystyle \dfrac{x}{t} = 2 A_0 \left[1 +m - \dfrac{2m(1-m)K(m)}{E(m) -(1-m)K(m)} \right] , \quad -6A_0 < \dfrac{x}{t} < 4A_0 , \\ \displaystyle A = 2A_0 m , \quad D =A_0 \left[ m-1 + \dfrac{2E(m)}{K(m)} \right] , \quad \varGamma = \gamma k = A_{0}^{1/2} . \end{array}\right\} \end{equation}

This describes an undular bore, an expanding wave train connecting a zero level at the front $x/t = 4 A_0 , D=0$, where $m \to 1$, to the rear $x/t = -6 A_0 D =A_{0} > 0$, where $m \to 0$. At the front the leading wave is a solitary wave of amplitude $2A_0$ moving in the positive $x$-direction, and at the rear the waves are linear sinusoidal waves moving in the negative $x$-direction. Note that both $x, t$ can be translated to different origins, leading to additional flexibility in fitting to initial conditions. This solution has been widely used as a model for nonlinear wave trains in shallow water and in density stratified flows; see Grimshaw & Smyth (Reference Grimshaw and Smyth1986), Smyth (Reference Smyth1987) and Grimshaw & Yuan (Reference Grimshaw and Yuan2016, Reference Grimshaw and Yuan2018) for instance. With the forcing/friction term $F(\eta ) \ne 0$ the system ((3.2), (3.16), (3.19)) is difficult to handle analytically, and so we simplify it using the solitary wave train approximation.

In the solitary wave train approximation (3.13a–c) when $m \to 1$, $\gamma \to \infty, k \to 0$ while $\varGamma =\gamma k$ stays finite, $E \to 0$ but $\mathcal {E} = E/k$ remains finite and is given by

(3.24)\begin{equation} \mathcal{E} = \frac{2A^2}{3\varGamma } = \frac{8\varGamma^3 }{3}= \frac{(2A)^{3/2}}{3} . \end{equation}

Since in this limit $E \to 0$ the mean level (3.16) reduces to (3.17) and yields a solution $D(x,t)$ independently of other parameters. Hence, $D$ can be regarded as known. In order to simplify the subsequent analysis of a solitary wave train, we shall set $D = 0$ in the detailed calculations to follow in the remainder of this section.

In this solitary wave train limit, $\mathcal {U} \to 0$, and the wave action (3.20) becomes

(3.25)\begin{equation} \mathcal{E}_t + V\mathcal{E}_x + 6\mathcal{E} D_x = \mathcal{F}_S . \quad V = 6D + 4\varGamma^2 . \end{equation}

Here, $\mathcal {F}_S$ is the solitary wave train limit of $\langle \hat {\eta } \hat {F} \rangle /k$, given by

(3.26)\begin{equation} \mathcal{F}_S =\mu_1 \frac{16A^2 \varGamma}{15} +\mu_2 \left[D^2\frac{16 A^2 \varGamma}{5} + D\frac{256 A^3 \varGamma}{35} +\frac{128A^4 \varGamma}{105 } \right] - \mu_3 \frac{4A^2}{3\varGamma } . \end{equation}

The solitary wave expressions (3.13a–c) have been used to express the $\mathcal {F}_S$ in terms of $\varGamma$ noting that $A = 2\varGamma ^2$. Since (3.25) is an equation for $E$ or $A$, or $\varGamma$ alone the system ((3.2), (3.16), (3.19)) has been reduced to a single equation. The wavenumber $k$ is then determined by the conservation of waves equation (3.2) with $V$ given in (3.25). The solitary wave train modulation equation (3.25) is a single equation for either $\mathcal {E}$, $A$ or $\varGamma$. Here, we choose $\varGamma$ and then, after putting $D=0$, (3.25) becomes

(3.27)\begin{gather} \varGamma_t + 4\varGamma^2 \varGamma_x = \mathcal{G} , \end{gather}

(3.28)\begin{gather}\displaystyle \mathcal{G} = \frac{\mathcal{F}_S}{8 \varGamma^2 } =\mu_1 \frac{8}{15} \varGamma^3 + \mu_2 \frac{256}{105} \varGamma^7 -\mu_3 \frac{2\varGamma }{3} . \end{gather}

Equation (3.27) is a nonlinear hyperbolic equation which nonetheless in general is quite difficult to solve explicitly. When $\mathcal {G}=0$ it reduces to a Hopf equation and solutions can be found using standard techniques. In particular, we note the similarity solution

(3.29)\begin{equation} A = 8\varGamma^2 = 2\frac{x}{t} . \end{equation}

This solution exhibits a linear progression in amplitude of a solitary wave train. Along the track of a specific soliton with amplitude $A_n$ and speed $V_n = A_n /2$, $x/t \to V_n$ as ${t \to \infty}$. There is a superficial resemblance to the undular bore solution (3.23) and it can be interpreted as related to that when $m \approx 1$. Like that solution, $x/t$ can be replaced by $(x-x_0 )/(t- t_0 )$ for arbitrary constants $x_0, t_0$ provided that it stays positive.

When $\mathcal {G} \ne 0$ we first consider solutions which depend on $t$ alone, when the solution is given by

(3.30)\begin{equation} \mathcal{H}(\varGamma ) = \int_{\varGamma_0 }^{\varGamma } \frac{{\rm d} \varGamma }{\mathcal{G}} = t, \end{equation}

where $\varGamma _0 \ne 0$ is the initial constant value. On its own, the first term in $\mathcal {G}$ proportional to $\mu _1$ generates the solution

(3.31a,b)\begin{equation} \frac{1}{2\varGamma_0^2} -\frac{1}{2\varGamma^2} = \frac{8 \mu_1 t}{15} , \quad A =2\varGamma^2 = \frac{A_0}{1 - A_0 8 \mu_1 t/15} . \end{equation}

When $\mu _1 < 0$, $A, \varGamma \to 0$ as $t \to \infty$, but when $\mu _1 > 0$, $A, \varGamma \to \infty$, albeit in a very long but finite time proportional to $\mu _{1}^{-1}$. This singularity in a modulated solitary wave in the unstable KdV–Burgers equation was noted by Manna & Latifi (Reference Manna and Latifi2023) but we emphasise it is a singularity in the modulation asymptotic theory, and not in the full equation. Again on its own, the second term proportional to $\mu _2 > 0$ generates the solution

(3.32a,b)\begin{equation} \frac{1}{6\varGamma_0^6} -\frac{1}{6\varGamma^6} = \frac{256 \mu_2 t}{105} , \quad A^3 = 8\varGamma^6 = \frac{A_{0}^3}{1 - A_{0}^3 64 \mu_2 t/35} , \end{equation}

with unlimited growth to infinity in a very long time proportional to $\mu _{2}^{-1}$. Also on its own the third term proportional to $\mu _3$ generates the solution

(3.33a,b)\begin{equation} \varGamma = \varGamma_0 \exp{({-}2\mu_3 t /3)}, \quad A =2\varGamma^2 = A_0 \exp{({-}4\mu_3 t/3)}, \end{equation}

with exponential decay to zero as $t \to \infty$. In general, $\mathcal {J}(\varGamma ) = \mathcal {G}(\varGamma )/\varGamma$ is initially negative at $\varGamma =0$, then increases or decreases according as $\mu _1 \gtrless 0$ and becomes zero at $\varGamma = \varGamma _a$, after which it remains positive and tends to infinity as $\varGamma \to \infty$. The expressions ((3.27) and (3.30)) then imply that, if $\varGamma _0 < \varGamma _a$, $\varGamma \to 0$, while if $\varGamma _0 > \varGamma _a$ then $\varGamma \to \infty$.

In general, (3.28) can be reduced to an inhomogeneous first-order nonlinear hyperbolic equation which can, in principle, be solved implicitly, using characteristics

(3.34)\begin{equation} \mathcal{H}_t + V\mathcal{H}_x = 1 , \quad V = V(\mathcal{H}) = 4\varGamma^2 . \end{equation}

Solutions with zero $x$-dependence were discussed in the previous paragraph. Here, $\mathcal {H}$ is a function of $\varGamma$, which can in principle be inverted to express $\varGamma$ as a function of $\mathcal {H}$. Although this cannot be done explicitly, we have expressed $V = V(\mathcal {H})$. The inhomogeneous term on the right-hand side is removed by the transformation

(3.35)\begin{equation} \mathcal{H} = \hat{\mathcal{H}} + t , \quad \hat{\mathcal{H}}_t + V\hat{\mathcal{H}}_x = 0 , \quad V =V( \hat{\mathcal{H}} + t ) . \end{equation}

This is a Hopf type of equation and we expect the general solution to exhibit wave steepening. It can be formally solved using characteristics, that is, $\hat {\mathcal {H}}$ is constant on the curves ${\textrm {d}\kern0.07em x}/\textrm {d}t = V$. This inevitably leads to wave steepening on a time scale determined by the slope of the initial condition for $A =2\varGamma ^2$, which may be shorter than the time scales in the $x$-independent solutions ((3.31a,b)–(3.33a,b)).

Finally, we consider the linear sinusoidal wave limit (3.13a–c) when $m \to 0$ and the left-hand side of (2.4) is linearised. In this limit, only the linear part proportional to $\mu _1$ in the forcing/friction terms is retained. The parameter set is reduced to three $[A, k, D]$ noting that the phase velocity $V =6D - k^2$. The mean (3.16) reduces to $D_t =0$ and so we put $D = D_0$, taken here as a constant. In (3.2) for conservation of waves $\omega = kV = 6kD_0 - k^ 3$ and the linear group velocity is $c_g = \omega _k$, so that

(3.36)\begin{equation} k_t + c_g (k) \,k_x = 0 , \quad c_g = 6 D_0 - 3 k^2 . \end{equation}

Equation (3.36) is a Hopf equation, which can be solved implicitly. In general, $c_g (k)= c_g (x,t)$ and we define a new variable $X$ by

(3.37)\begin{equation} \frac{{\rm d}\kern0.07em x}{{\rm d}t} = c_g (x,t), \quad x = X \ \hbox{at} \ t =0 , \end{equation}

with solution $x = x(X,t)$. Using (3.36), we find that $c_g = c_g(X)$ and then

(3.38)\begin{equation} x = X + c_g(X)t . \end{equation}

The solution of (3.36) is

(3.39)\begin{equation} k = k_{0} (X), \quad k (x, t=0) = k_{0} (x) , \end{equation}

with a corresponding solution for $c_g = c_g(k)$. Although (3.38) determines $x =x(X, t)$ explicitly, the inverse relation $X=X(x, t)$ is implicit. Since $X_x = (1 + c_{gX} t )^{-1} ,$ $X_t + c_g X_x =0$ a singularity develops when $1 + c_{gX} t =0$. This only occurs if $c_{gX} < 0$, or since $c_{gX} = -6kk_X$, if $k_X > 0$.

The wave energy expressions (3.12) reduces to

(3.40a,b)\begin{equation} E = \frac{ A^2 }{16} , \quad \mathcal{U} ={-}2k^2 E , \end{equation}

and the wave energy equation (3.19) becomes

(3.41)\begin{equation} E_t + ( c_g (k) E )_x = \mu_1 2k^2 E . \end{equation}

Here, $E$ propagates with the group velocity $c_g$ modulated by the KdV–Burgers term with coefficient $\mu _1$, which can be either positive or negative. This is a linear hyperbolic equation for $E$ and, using the transformation $x=(X, t)$, $E (X, t)$ satisfies

(3.42) \begin{equation} \left.\begin{array}{c@{}} E_t + \dfrac{ c_{gX}}{1 +c_{gX}t}E = \mu_1 2k^2 E , \quad E(X, t=0) = E_0 (X) , \\ \text{so that} \quad E = \dfrac{E_0(X)}{1 +c_{gX}t }\exp{( \mu_1 2k^2 t )} . \end{array}\right\} \end{equation}

The singularity at $1 + c_{gX} t =0$ is inherited from those in $k, c_g$ and otherwise $E$ grows or decays exponentially according as $\mu _1 \gtrless 0$.

Finally, we note that the linear wave modulation system ((3.36), (3.41)) has a steady solution in which $k, c_g$ are constants and

(3.43)\begin{equation} E= E (x) =E_0 \, \exp{( \nu x )}, \quad \nu = \frac{\mu_1 2k^2 }{c_g }. \end{equation}

This is exponentially increasing or decreasing with $x$ according to $\mu _1c_g \gtrless 0$. It could be a partial, small-amplitude representation of the steady-state bore solutions analysed by Johnson (Reference Johnson1970) for the KdV–Burgers equation when $\mu _1 < 0$. However, we will not pursue that here as it is beyond the scope of this article.

4. Numerical simulations

In this section we describe numerical simulations of the non-dimensional KdV–Burgers equation (2.4), repeated here,

(4.1)\begin{gather} \eta _t + 6\eta \eta_x + \eta_{xxx} = F(\eta) -s(x)\eta , \end{gather}

(4.2)\begin{gather}F(\eta ) ={-}\mu_1 \eta_{xx} - \mu_2 (\eta^3)_{xx} -\mu_3 \eta , \end{gather}

(4.3)\begin{gather}\mu_{1,2,3} = \left\{K \left[3\frac{\rho_{a}}{\rho_w} \beta \frac{W_{r}^2}{C^2} - \frac{12 \kappa}{hC} \right] , \frac{2 }{K}\frac{\rho_a}{\rho_w} c_d \frac{hU_{a}^2 }{A_0 C^2 }, 2K f_r C_D \frac{A_0 }{h} \right\}. \end{gather}

The numerical scheme is a Fourier spectral method, similar to that used in Grimshaw & Maleewong (Reference Grimshaw and Maleewong2013, Reference Grimshaw and Maleewong2016, Reference Grimshaw and Maleewong2019) for the topographically or pressure forced KdV equation. The numerical solution is moved forward time using the fourth-order Runge–Kutta scheme with a time step $\delta t = 0.0002$. In most cases, we set $L=451 {\rm \pi}$ with $8192$ Fourier modes. The term $s(x) \eta$ is a sponge layer to absorb waves approaching the boundaries $x = \pm L$, defined by

(4.4)\begin{equation} s(x) = \frac{s_1}{2} \{ (1+\tanh(s_2(x-x_m))) + (1-\tanh(s_2(x+x_m)))\} . \end{equation}

We set $s_1 = 10$, $s_2=0.1$ and $x_m = L - 2$, so that $\mp x_m$ are near the left and right boundaries. The sponge layer changes the mass and energy expressions. For instance, the mass conservation law (3.16) for a modulated periodic wave becomes

(4.5)\begin{equation} D_t + 6\mathcal{M}_x = \langle F(\eta ) \rangle ={-}\mu_3 D - s(x)D. \end{equation}

Similarly, the wave energy equation (3.19) becomes

(4.6)\begin{equation} E_t + (VE + \mathcal{U})_x +6ED_x = \langle \hat{\eta} \hat{F} \rangle - 2s(x)E . \end{equation}

In both cases, the extra dissipation induced by the sponge layer should only be effective near the boundaries. However, the global mass and energy laws ((2.7a,b) and (2.10)) both show overall dissipation due to the sponge layer.

The initial condition is $\eta (x, t=0) = \eta _{0}(x)$, a wave packet similar to that used in MG24

(4.7)\begin{equation} \eta_{0}(x) = ENV(x)(D_0 + M\cos{( k_0 x )}) , \quad ENV(x) = \hbox{sech}(x / L_E) , \quad -L < x < L , \end{equation}

where $M$ is the constant carrier wave amplitude, $k_0 > 0$ is the constant initial wavenumber and $D_{0}$ is the constant initial mean level; $ENV(x)$ is a slowly varying wave envelope, which tends to zero as $x \to \pm L$. The initial condition (4.7) is expressed in non-dimensional variables. The scale for the wave amplitude is $E_S =2h/3K^2$, so with $h=10 \ \textrm {m}$, we put $K = 2.58$ and then $E_S = 1.0015 \approx 1$, and $\eta _{n} m \approx \eta _{d}$, where subscripts $n$ and $d$ refer to non-dimensional and dimensional variables, respectively. The envelope span $L_E =60$ is chosen so that $2{\rm \pi} /k_0 \ll L_E \ll L = 451 {\rm \pi}/k_0$, where the parameter $k_0 =1$ for a $5$-second wave. The envelope implies that the initial wave amplitude is really $A_{0}(x) = ENV(x) M$ and the initial mean level is really $D_0 (x) = ENV(x) D_0$. The carrier wave amplitude varies over the values $M =0.25, 0.5, 0.75, 1.0$ and $D _0$ is a constant with values $0, \pm M$. Although the initial condition (4.7) is a sinusoidal wave, the modulus $m$ of the emerging cnoidal waves can be estimated from the discussion in § 3.1. We find that $M \approx m K(m)^2/{\rm \pi} ^2$ and the modulus $m \to 1$ as $M$ increases. Even for $M =0.25$, $m \approx 1$, as our simulations will show. For the KdV equation, a balance between nonlinearity and dispersion requires that $3M \sim k^2 \sim 1$, which is achieved here. A summary of the various parameters is in table 1, which yields the values $\mu _1 = 0.00014, \ \mu _2 = 0.000011, \ \mu _3 = 0.00035$ for these dimensionless parameters. We carried out many simulations and only a representative sample are shown here. We set $\mu _1 = \pm 0.001, \pm 0.0001$, $\mu _{2,3} = 0.001, 0.0001$, chosen around the values found above using the specified values in table 1 in order to explore the dependence on these key parameters. For the wave parameters we set $M=0.25, 0.5, 0.75, 1.0$, $D_0 = 0, \pm 0.5$. Some simulations were run for a case of no initial wave envelope, $ENV(x) \equiv 1$ and some with no sponge layer, $s(x) \equiv 0$. All simulations used a periodic boundary condition.

Table 1. Parameters.

First, we consider the unforced KdV equation (4.1) when $F(\eta ) = 0$, that is, when $\mu _{1,2,3}=0$. With $D_0 =0$ we simulated four cases $M = 0.25, 0.5, 0.75, 1.0$. Next, we considered some forced cases by setting the values for each $\mu _{1,2,3}$ separately, that is, $\mu _1 \ne 0, \mu _{2,3}=0$, and so on, followed by a full combination of these. After that we simulated cases with $D_0 = \pm M$ a non-zero constant. This range of values for $M, D_0$ and those for each $\mu _i = 0,\pm 0.0001, \pm 0.001$ are based on data for a 5-second wave in MG24. Although for the parameters in table 1 $\mu _1 >0$ we performed simulations for both $\mu _1$ positive and negative. Many other simulations were performed with various parameter settings and with and without the sponge layer. Although only a sample are shown here, the interpretation is based on many simulations.

An unforced case $\mu _{1,2,3}=0$ is shown in figure 2 for $D_0 =0$ and $M=0.5$ (a) and $M=1.0$ (b); contour plots of $\eta$ with $M=0.25, 0.5, 0.75, 1.0$ are shown in figure 3 on the upper left, upper right, bottom left and bottom right, respectively. The initial wave packet generates transient radiation moving upstream ($x<0$) and a solitary wave train with $N$ solitons moving downstream ($x>0$) following the arrow directions as indicated in figure 2. This outcome is well known and to be expected in the KdV equation; see El (Reference El2007) and Osborne (Reference Osborne2010) for instance in the water wave context. There are many solitons and the solitary wave train resembles a soliton gas (El Reference El2021). The number $N$ of solitons can be estimated from the inverse scattering transform theory, see El (Reference El2007) for instance; there are no solitons ($N =0$) if $\eta _{0}(x) < 0$ for all $x$, and otherwise $N \ge 1$

(4.8)\begin{equation} N \leqq 1 + I_N, \quad I_N = \int_{-\infty}^{\infty} |x| | \eta_{0}(x)| \,{{\rm d}\kern0.07em x} . \end{equation}

For the initial condition (4.7) with $D_0 =0$ we estimate that $I_N < 4M L_E^2$, which, although it is certainly an overestimate, indicates correctly that $N$ is very large due to the envelope span $L_E = 60 \gg 1$. Larger values of $M$ produce higher wave amplitudes both upstream and downstream. The radiation field is mainly determined by the envelope $ENV(x)$ demonstrated by simulations when the envelope is replaced with a constant, $ENV(x)=1$. This case is shown in figure 4(a) for $M=1$, $D_0=0$, $\mu _{1,2,3} = 0$. This case replicates the famous pioneering work of Zabusky & Kruskal (Reference Zabusky and Kruskal1965) when only one sinusoidal wave was in the initial condition, whereas our simulations had a sinusoidal wave train. Cases $\mu _1=0.001$ and $\mu _2=0.001$ are shown in figures 4(b) and 4(c), respectively, showing the effect of each forcing term leading to instability, indicated around $t=50$ and $t=25$ with high oscillations, increasing wave speeds and amplitudes. The case $\mu _3=0.001$ is shown in figure 4(d), where the wave speed and amplitude are decreasing.

Figure 2. Evolution of $\eta$ with initial condition (4.7) for $M=0.5$ (a) and $M=1.0$ (b) when $D_0 =0$ and $\mu _1=\mu _2=\mu _3=0$. Arrows show the direction of the moving waves.

Figure 3. Contour plots of $\eta$ with initial condition (4.7) for $M=0.25$ (a), $M=0.5$ (b), $M=0.75$ (c) and $M=1.0$ (d), with $D_0 =0$, when $\mu _1=\mu _2=\mu _3=0$.

Figure 4. Contour plots of $\eta$ with initial condition (4.7) when $ENV(x)=1$, $M=1$ and $D_0=0$. Case ${\mu _{1,2,3}=0}$ (a), case $\mu _1 = 0.001$ $\mu _2=\mu _3=0$ (b), case $\mu _2 = 0.001$ $\mu _1=\mu _3=0$ (c), case $\mu _3 = 0.001$ $\mu _1=\mu _2=0$ (d).

To investigate the effects of the forcing terms on the right-hand side of (4.1), we fix $M=1, L_E =60$ and vary the magnitude of each $\mu _{1,2,3}$ from values based on our wind shear parameters to slightly larger values to explore the effects on a reasonable time scale. Then, we compare the results with the unforced case and with the asymptotic theory of § 3. The evolution of $\eta$ when $\mu _1=0.001 > 0, \mu _2 = \mu _3 = 0$ is shown in figure 5(a) and the corresponding contour plot is shown in figure 5(b). The solitary wave train grows in amplitude, especially downstream. In contrast the evolution of $\eta$ when ${\mu _1=-0.001 < 0}, \mu _2 = \mu _3 = 0$ is shown in figure 6(a) and the corresponding contour plot is shown in figure 6(b). The solitary wave trains are now decaying, with decreasing wave speeds. Next, the evolution of $\eta$ when $M=1$, $\mu _2=0.001 > 0, \mu _1 = \mu _3 = 0$ is shown in figure 7(a) and the corresponding contour plot is shown in figure 7(b). The solitary wave train now grows in amplitude. Finally, the evolution of $\eta$ when $\mu _3=0.001 >0, \mu _1 = \mu _2 = 0$ is shown in figure 8(a) and the corresponding contour plot is shown in figure 8(b). The solitary wave train and the radiation field now decay, especially downstream.

Figure 5. Evolution of $\eta$ with initial condition (4.7) for $M=1.0$, $D_0 =0$ when $\mu _1=0.001, \mu _2=0, {\mu _3=0}$ (a), and contour plot of $\eta$ (b).

Figure 6. Evolution of $\eta$ with initial condition (4.7) for $M=1.0$, $D_0 =0$, when $\mu _1= -0.001, \mu _2=0, \mu _3=0$ (a), and contour plot of $\eta$ (b).

Figure 7. Evolution of $\eta$ with initial condition (4.7) for $M=1.0$, $D_0=0$, when $\mu _1=0, \mu _2=0.001, \mu _3=0$ (a), and contour plot of $\eta$ (b).

Figure 8. Evolution of $\eta$ with initial condition (4.7) for $M=1.0$, $D_0 =0$, when $\mu _1=0, \mu _2=0, \mu _3=0.001$ (a), and contour plot of $\eta$ (b).

Next, in order to compare these simulations with the unforced cases and with the modulation theory of § 3 we plot the tracking of the maximum value of $\eta$ over time in figure 9. Although this only reveals the total maximum and not the individual wave amplitudes, in most cases the maximum was the amplitude of the leading solitary wave. The agreement is surprisingly good, given that the theory is only for time-dependent modulations. The cases when we vary $\mu _1$ with $\mu _2=\mu _3=0$ are shown on the upper left of figure 9. When $\mu _1>0$ there is wave growth due to critical level instability, while $\mu _1<0$ leads to decaying waves due to dissipation at the water surface. The value of $\mu _1 > 0$ is of order $10^{-4}$ based on the value of $\beta = 1.745$ found in MG24, but $\mu _1 < 0$ is much smaller in magnitude as the kinematic viscosity $\kappa = 10^{-6}\ \textrm {m}^{2}\ \textrm {s}^{-1}$. The predicted growth or decay rates from the modulation theory, (3.31a,b) for $\mu _1 =0.001$ and $-0.001$, are plotted to compare with the numerical results. Note that the simulations cannot be run for too long due to the dissipative effect of the numerical sponge layers. The cases with $\mu _2>0$ and $\mu _1=\mu _3=0$ are shown on the upper right of figure 9. The wind stress effect in $\mu _2$ leads to wave growth as expected, although $\mu _2$ has an order of magnitude $10^{-4}$, which is small compared with $\mu _1$. If $\mu _2$ is of the order of $10^{-3}$ due to a much smaller depth, then there will be significant wave growth. The predicted growth rates from (3.32a,b) for $\mu _2 =0.0001$ and $0.001$ are plotted to compare with the numerical results and again they agree well. The cases describing the effect of bottom stress when $\mu _3>0$ and $\mu _1=\mu _2=0$ are shown on the bottom left of figure 9. For the present parameters with depth $h=10\ \textrm {m}$, $\mu _3$ has order $10^{-4}$ and causes significant wave decay, which will increase as the depth $h$ decreases. The predicted decay rates from (3.33a,b) for $\mu _3 =0.0001$ and $0.001$ are plotted and again agree well. A case when all the parameters $\mu _{1,2,3}$ have the same order of magnitude $0.001$ is shown at the bottom right of figure 9. For these parameters the bottom friction term dominates and leads to wave decay.

Figure 9. With the initial condition (4.7), the tracking maximum value of $\eta$ over time when ${M=1}, {D_0 =0}$, compared with the unforced case $\mu _{1,2,3}=0$. Cases $\mu _1 \ne 0$, $\mu _2=\mu _3=0$ (a), cases $\mu _2 \ne 0$ $\mu _1= \mu _3=0$ (b), cases $\mu _3 \ne 0$, $\mu _1=\mu _2=0$ (c). The respective predictions ((3.31a,b)–(3.33a,b)) are from the modulation theory. A combined case $\mu _1=\mu _2=\mu _3=0.001$ is shown in (d).

Next, we examine cases with a non-zero initial mean level. We set $M=0.5$ and $D_0 = \pm 0.5$, a non-zero constant in the initial condition (4.7). Figure 10 with $M= 0.5$, $D_0 =0.5$ and $\mu _1 = \mu _2 =\mu _3 =0$ can be compared with figures 2(a) and 3(b), where $D_0 = 0$. Again, there is formation of a soliton gas downstream, now with increased $N$, larger amplitudes and faster speeds because the magnitude of $\eta _{0}(x)$ has been increased. Figure 11 is similar to figure 10 but with $\mu _1 = 0.001, \mu _2 =\mu _3 =0$ and can be compared with figure 5. Since $\mu _1 =0.001 > 0$ the soliton amplitudes now grow. The corresponding plots of the tracking of the maximum value of $\eta$ over time are shown in figure 12 and compared with the predictions from ((3.31a,b)–(3.31a,b)), which continue to hold for a constant $D_0$, noting that, in the absence of the wave envelope, a constant $D_0$ is just a translation of $x$. The case $\mu _1=\mu _2=\mu _3=0.001$ with a positive initial mean level $D_0=0.5$ and $M=0.5$ is shown in figure 12(d). The downstream wave amplitude grows exponentially, showing the dominant effect of the wind forcing. In the simulations, the initial value for $D$ is $D_0 (x) = D_0 ENV(x)$. When $D_0 >0$, an expansion fan forms upstream in $-\infty < x < 6D_0 t$ and wave steepening occurs downstream in $x < 6D_0 t < \infty$ with formation of a front. The front evolves into a solitary wave train as in figures 11 and 12. This is in contrast to the case of zero initial mean $D_0 =0$ when $M=1.0$ in figure 9(d), which shows wave decay due to the dominance of the term in $\mu _3$. We also varied $D_0$ over the range $[0.1, 0.4]$ and the tracking maximum amplitude shown in figure 12(d) decays very slightly. However, when $D_0=0.5$, the maximum amplitude increases. Seeking an explanation of this phenomena, we refer to the modulation theory (3.17) for $D$. First, (3.17) holds only when $D^2 \gg 2E$, otherwise (3.16) is the modulation equation for $D$ and this is affected by terms in $E_x$ generating decay if $E_x > 0$; examination of figures 2–9 when $D_0 =0$ suggests this is the case. If the linear sinusoidal limit $E =A^2/16= M^2/4$ (3.40a,b) is used to estimate $E$ initially, then (3.17) holds when $D_0 \gg 0.7 M$; for $M =0.5$, this is $D_0 \gg 0.35$, consistent with the numerical simulations of figure 12(d). Second, the modulation solution for $D$ from (3.17) shows that $D$ will decay if $\mu _3 +6D_{0x} > 0$ for all $x$, but otherwise will increase. Taking account of the envelope $\hbox {max}(-6D_{0x}) =3D_0 / L_E < \mu _3$ when $D_0 <\mu _ 3 L_E / 3$. For the parameters we use, that is $D_0 < 0.02$. Although much less than the numerical predictions, the trend is correct.

Figure 10. Evolution of $\eta$ with initial condition (4.7) for $M=0.5$, $D_0 = 0.5$, when $\mu _1= \mu _2=\mu _3=0$ (a), and contour plot (b).

Figure 11. Evolution of $\eta$ with initial condition (4.7) for $M=0.5$, $D_0 = 0.5$, when $\mu _1= 0.001, \mu _2=\mu _3=0$ (a), and contour plot (b).

Figure 12. With the initial condition (4.7), the tracking maximum value of $\eta$ over time when $M =0.5$, $D_0=0.5$, compared with the unforced case $\mu _{1,2,3}=0$. Case $\mu _1 \ne 0$, $\mu _2=\mu _3=0$ (a), case $\mu _2 \ne 0$, ${\mu _1=\mu _3=0}$ (b), case $\mu _3 \ne 0$, $\mu _1=\mu _2=0$ (c). The respective predictions ((3.31a,b)–(3.31a,b)) are from the modulation theory. A combined case $\mu _1=\mu _2=\mu _3=0.001$ with varying $D_0$ is shown in (d).

Figures 13 and 14 show cases with $M=0.5, D_0 =-0.5$ and $\mu _{1,2,3 } =0$ and $\mu _1 =0.001, \mu _{2,3} =0$, respectively. In these cases $\eta _{0}(x) \le 0$ and so no solitons are expected to form. All radiating waves travel upstream with a larger amplitude when $\mu _1 > 0$. The corresponding plots of the tracking of the maximum value of $\eta$ over time are shown in figure 15. As above we use the modulation theory (3.17) to interpret these results. Recalling that in the simulations the initial value for $D$ is $D_0 (x) = D_0 ENV(x)$, when $D_0 < 0$ we get the opposite scenario to $D_0 > 0$. An expansion fan forms downstream in $6D_0 t < x < \infty$ and a front forms upstream in $-\infty < x < 6D_0 t$ with wave steepening and front formation. The fronts form first where $ENV_x$ has a maximum value, that is, at $x = \pm x_M = D_0 t$, where $x_M = 0.88 L_E$, $\sinh {(x_M/L_E)} =1$. The fronts are artefacts of the modulation theory and in practice are smoothed by higher-order dispersion terms. This interpretation is consistent with the numerical simulations. The front evolves into a bore-like feature as in figures 13 and 14, but solitons are now precluded from forming. Instead, a quasi-periodic wave train forms behind the front. Figures 13 and 14 also show intensifying radiation far upstream, from which we infer that the wavenumber $k$ is increasing in magnitude. We suggest this is also due to the mean level upstream steepening. The wavenumber $k$ is described by (3.2), which reduces to (3.36) for linear sinusoidal waves when $V =6D_0 - k^2$. Since $D_0$ is a constant, $k=k_0$ is constant, but when $D_0$ is replaced by $D$ due to the envelope and varies with $x,t$ as described above, $k >0$ must also vary and we observe an increase in $k$ associated with the formation of fronts in $D$. The intensification moves upstream as increasing $k$ leads to an increase in the speed, $V=6D - k^2$, for a linear sinusoidal wave. Although then the solution for $D$ is available from (3.17), an explicit analytic expression for $k$ from (3.2) is not available. However, as $t \to \infty$, $|D|$ will decrease, leading to a steady state far upstream, as seen in figures 13 and 14. If we suppose that the variation in $V = 6D -k^2$ is due mainly to $D$ then (3.2) yields $k_t \approx -k6D_x$; when $D_0 < 0$ and so $D_{x} > 0$ upstream, there is an increase in $k$ as observed, which is not predicted when $D_0 > 0$. The simulations were repeated for $M =0.5$, $D_0 = 0, \pm 0.5$ but with an increased $k_0 =2$. The outcome (not shown here) was unchanged except that the upstream radiation, especially for $D_0 =-0.5 < 0$, was even more intense. As a check on another possible interpretation, we moved the upstream sponge layer to far downstream to test whether it was involved in this feature, but we found no discernible difference.

Figure 13. Evolution of $\eta$ with initial condition (4.7) for $M=0.5$, $D_0 = -0.5$, when $\mu _1=\mu _2=\mu _3= 0$ (a), and contour plot (b).

Figure 14. Evolution of $\eta$ with initial condition (4.7) for $M=0.5$, $D_0 = - 0.5$, when $\mu _1= 0.001, \mu _2=\mu _3=0$ (a), and contour plot (b).

Figure 15. With the initial condition (4.7), the tracking maximum value of $\eta$ over time when $M =0.5$, $D_0=-0.5$, compared with the unforced case $\mu _{1,2,3}=0$. Cases $\mu _1 = \pm 0.001$, $\mu _2=\mu _3=0$ (a), case $\mu _3 = 0.001$, $\mu _1=\mu _2=0$ (b).

5. Summary and conclusions

In this paper we have used the well-known KdV equation, supplemented by the addition of various forcing/friction terms, to describe the evolution of wind-generated water wave packets in shallow water. The outcome is a type of KdV–Burgers equation, given in dimensional form by ((2.1) and (2.2)) and in canonical non-dimensional form by ((2.4) and (2.5)) in § 2. The forcing/friction terms are the shallow-water limit of terms derived by Maleewong & Grimshaw (Reference Maleewong and Grimshaw2024) (MG24) and are conveniently described here in (2.6). There are three terms each with a non-dimensional coefficient $\mu _{1,2,3}$. The term with coefficient $\mu _1$ is a combination of wave growth due to the critical level instability theory of Miles (Reference Miles1957) and wave decay due to laminar friction at the water surface; $\mu _1$ can be positive or negative, describing either an unstable KdV–Burgers equation with wave growth, or a stable KdV–Burgers equation with wave decay. The term with coefficient $\mu _2$ describes wave growth due to a representation of wave stress in the air at the water surface induced by a turbulent wind. The term with coefficient $\mu _3$ describes wave decay due to a representation of a turbulent bottom boundary layer. The modified KdV–Burgers equation (2.4) has two principal conservation laws for mass and energy, expressed here by (2.4) itself and (2.9), respectively. In the absence of forcing/friction, wave energy is conserved; otherwise, wave energy either grows or decays depending on the relative signs and magnitudes of $\mu _{1,2,3}$. The term with coefficient $\mu _1$ leads to wave growth or decay according to whether $\mu _1$ is positive or negative; the term with coefficient $\mu _2 > 0$ leads to wave growth; the term with coefficient $\mu _3 > 0$ leads to exponential decay. For wind-induced wave packets, the rate of change of wave energy for wave growth is proportional to the square of the wave slope. Thus the wave amplitude grows slowly over time for initially small-amplitude waves, but then grows rapidly when the amplitude becomes large. The strength of wind shear and the wind speed near the water surface appear in the parameters $\mu _1, \mu _2$, implying that a stronger wind has a greater effect on water wave growth if $\mu _1>0$. We found that $\mu _1$ decreases from positive to negative as the depth decreases, $\mu _2$ is independent of the depth and $\mu _3$ increases rapidly as the depth decreases.

We pursued two approaches to investigate wave packet solutions of the KdV–Burgers equation (2.4). First, in § 3 we adapted the Whitham modulation theory to study slowly varying periodic wave trains, with an emphasis on the solitary wave train limit. We applied a wave average to the conservation laws for mass and energy, supplemented by (3.2) for conservation of waves. We obtained analytical solutions for modulated solitary wave trains, amongst which we noted especially the expressions ((3.31a,b)–(3.33a,b)) for time-dependent modulations of a solitary wave due to the forcing/friction terms with coefficients $\mu _{1,2,3}$, respectively. Second, in § 4, we simulated (4.1) numerically with the wave packet initial condition (4.7) modelled by a sinusoidal wave on a mean level, with a slowly varying envelope. The simulations used a sponge layer (4.4) at each end of the domain to absorb outgoing waves. In general, the outcome was a solitary wave train downstream with many waves resembling a soliton gas, with some radiation upstream of sinusoidal waves. Where appropriate, the numerical solutions agreed well with the modulation theory for the wave growth and decay rates.

We conclude that, in this shallow-water limit, using the KdV–Burgers model with forcing/friction modifications, (2.4) provides a useful model for wind-generated waves. This has been noted in the literature, notably by Costa et al. (Reference Costa, Osborne, Resio, Alessio, Chriv, Saggese, Bellomo and Long2014) in an examination of ocean data in a shallow sea using just the KdV equation, and more recently by Zdyrski & Feddersen (Reference Zdyrski and Feddersen2021) and Manna & Latifi (Reference Manna and Latifi2023) using an unstable KdV–Burgers equation, numerically and analytically, respectively. The main disadvantage of this model is that it is uni-directional. For weak two-dimensional dispersion, the KdV part can be replaced by a Kadomtsev–Petviashvili (KP)-based model. More generally, to take account of both bi-directional behaviour in the wind direction and two-dimensional dispersion, a model based on a suitable Boussinesq equation could be developed.