Evolution of water wave packets by wind in shallow water

Abstract

We use the Korteweg–de Vries (KdV) equation, supplemented with several forcing/friction terms, to describe the evolution of wind-driven water wave packets in shallow water. The forcing/friction terms describe wind-wave growth due to critical level instability in the air, wave decay due to laminar friction in the water at the air–water interface, wave stress in the air near the interface induced by a turbulent wind and wave decay due to a turbulent bottom boundary layer. The outcome is a modified KdV–Burgers equation that can be a stable or unstable model depending on the forcing/friction parameters. To analyse the evolution of water wave packets, we adapt the Whitham modulation theory for a slowly varying periodic wave train with an emphasis on the solitary wave train limit. The main outcome is the predicted growth and decay rates due to the forcing/friction terms. Numerical simulations using a Fourier spectral method are performed to validate the theory for various cases of initial wave amplitudes and growth and/or decay parameter ranges. The results from the modulation theory agree well with these simulations. In most cases we examined, many solitary waves are generated, suggesting the formation of a soliton gas.

1. Introduction

The evolution of a wind-generated water wave is a fundamental and much-studied problem of both scientific and operational concern. Oceanic wind waves affect the weather and climate through transfer processes across the air–water interface, generate large forces on marine structures, ships and submersibles and lead to extreme events such as rogue waves and storm surges. But despite much theoretical research over many years, combined with more recent in situ observations and numerical simulations using powerful computers, the theoretical mechanism for wind-wave formation and evolution remains a challenging problem. This was very evident at the IUTAM Symposium on Wind Waves held in London in September 2017 (Grimshaw, Hunt & Johnson 2018), where a wide range of theories, simulations, observations and experiments were presented with a very lively discussion.

In particular there are only tentative theories about how wind affects the dynamics of wave groups, where the issue is how, in the presence of wind, do water waves form into characteristic wave groups, and what are their essential properties, depending on the local atmospheric and oceanic conditions; see Zakharov et al. (2015), Zakharov, Resio & Pushkarev (2017) and Zakharov (2018) amongst many similar critical comments. This is the issue we have been looking at; see Grimshaw (2018, 2019a,b) and Maleewong & Grimshaw (2022a,b, 2024), the latter hereafter referred to as MG24. The key feature in our analysis is that, in the linearised limit, a wave group moves with a real-valued group velocity even for unstable waves when the wave frequency and the wavenumber may both be complex valued. In the absence of wind forcing, it is well known that the nonlinear Schrödinger equation (NLS) describes wave groups in the weakly nonlinear asymptotic limit where wave groups are initiated by modulation instability and then represented by the soliton and breather solutions of the NLS model, see Grimshaw (2007) and Osborne (2010) for instance. The effect of wind forcing is captured by the addition of a linear growth term, based on the critical level instability theory of Miles (1957), which leads to the forced NLS, see Maleewong & Grimshaw (2022a,b, 2024). Recently, in MG24, we examined how the growth rate depended on the water depth and the parameters of the wind shear profile, utilising a weak frictional modification of potential flow introduced by Dutykh & Dias (2007b), subsequently used by Kharif et al. (2010) amongst others, and extending it to take account of critical level instability and wave stress turbulent forcing in the air near the water surface.

As part of the analysis in MG24 we reduced the fully nonlinear modified Euler system in the limit as the depth goes to zero to a modified version of the shallow-water equations. These lead to the possibility of wave steepening, which indicates the necessity for higher-order dispersive terms to control this steepening. In MG24 we analysed the modified shallow-water system using right-going and left-going Riemann invariants. We noted there that, in a one-way setting for the right-going waves, inclusion of higher-order dispersion would lead to a Korteweg–de Vries (KdV) equation, supplemented by the addition of forcing terms describing critical level instability and wave stress driving at the air–water interface, and frictional terms due to laminar friction at the water surface and turbulent friction at the water bottom. That step is the subject of this paper and the outcome is a KdV–Burgers type of equation. We note that retaining two-way propagation would lead to a Boussinesq system, as shown by Dutykh & Dias (2007a,b) and Dutykh (2009), who extended the frictionally modified shallow-water equations to such a Boussinesq system. We have not pursued that here, as the modulation theory for nonlinear waves that we use in § 3 is not as well developed for a Boussinesq equation as for the KdV equation. Other one-way alternatives such as the Benjamin–Bona–Mahony equation were not considered for similar reasons.

The KdV equation arises in a huge variety of physical contexts as it is the canonical equation expressing balance between weak nonlinearity and weak dispersion. Indeed, it is justifiably famous for that reason just as much as for its integrability, multi-soliton structures and astonishingly rich dynamics. It has been and continues to be widely used as a model for water waves in shallow water. It is therefore somewhat surprising that applications to wind-driven waves are relatively sparse. An important recent exception is the work by Costa et al. (2014), who used the multi-soliton dynamics to analyse wind wave data in a shallow bay. That analysis is based on the theory of a soliton gas; see El (2021) for a recent review which touches on the water wave application. As we will show, the extension to a KdV–Burgers type of equation builds on the well-known KdV dynamics.

In § 2 we present the KdV–Burgers equation with the additional forcing/friction derived in MG24. These forcing/friction terms are small, but important and essential perturbations to the KdV equation, and so in § 3 we describe modulation asymptotic theory for this perturbed KdV equation, focussing on the solitary wave train limit. In § 4 we present some numerical simulations of the full KdV–Burgers equation for an initial condition describing a wave packet. We conclude with a summary and discussion in § 5.

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 (1957), $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 (1957) 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 (1957), the term in $\alpha$ is omitted to derive the first term in (2.2); see Kharif et al. (2010). Finally, we put ${\rm i} k \to \partial /\partial x$. It is pertinent to note here that Miles (1957) 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 (1957) 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 (2007a,b), Dutykh (2009) and Kharif et al. (2010)); 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 (1995) and Tang et al. (1996) 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 (1995) and Tang et al. (1996) 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 (1995) and Tang et al. (1996) 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 (2007a,b)) 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 (1957) and Kharif et al. (2010) 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 (1969) 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 (1970) 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 (2021) and Manna & Latifi (2023) 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 (2021) used numerical simulation to examine wave shape, while Manna & Latifi (2023) 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 (1974) and Whitham (1974). It has been heavily used since, see for instance Grimshaw & Smyth (1986), Smyth (1987, 1988), Myint & Grimshaw (1995), Grimshaw (2010), El, Grimshaw & Tiong (2012) and Grimshaw & Yuan (2016, 2018) 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 (1986), El (2007), Grimshaw (2010) and Grimshaw & Yuan (2016). 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 (1979) and Grimshaw & Mitsudera (1993) 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 (1986), Smyth (1987), El (2007) and Grimshaw (2010); perturbative effects such as friction and variable depth have been examined by Smyth (1988), Myint & Grimshaw (1995), El et al. (2012) and Grimshaw & Yuan (2016, 2018), 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 (2016) 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 (1974) and Gurevich & PitaevskiI (1974) and the comprehensive review by El (2007). 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 (1986), Smyth (1987) and Grimshaw & Yuan (2016, 2018) 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 (2023) 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 (1970) 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 (2013, 2016, 2019) 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 (2007) and Osborne (2010) for instance in the water wave context. There are many solitons and the solitary wave train resembles a soliton gas (El 2021). The number $N$ of solitons can be estimated from the inverse scattering transform theory, see El (2007) 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 (1965) 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 (2024) (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 (1957) 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. (2014) in an examination of ocean data in a shallow sea using just the KdV equation, and more recently by Zdyrski & Feddersen (2021) and Manna & Latifi (2023) 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.