2.1. Nonlinear equations for water waves

The water waves are modelled using incompressible two-dimensional flow in the domain $-h < y< \eta$, where $y = -h$ is a rigid fixed bottom and $y =\eta (x,t)$ is the free surface. We adapt the weak modification of potential flow used by Longuet-Higgins (Reference Longuet-Higgins1992), Dutykh & Dias (Reference Dutykh and Dias2007) and Dias et al. (Reference Dias, Dyachenko and Zakharov2008) to account for thin laminar frictional boundary layers at the free surface and at the water bottom, and modified by Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010) to account for wave growth due to the Miles mechanism of critical-level instability. To these effects we add a parametrisation of wave stress due to turbulent air flow near the air–water interface, and a parametrisation of turbulent bottom stress, an alternative to laminar friction at the bottom. The equations are expressed in terms of a velocity potential $\phi (x, y, t)$, where the fluid velocity $(u,v) = ( \phi _x, \phi _y ) + (\psi _y , -\psi _x)$ and the terms in $\psi$ are there to take account of weak frictional effects:

(2.1)\begin{equation} \phi_{xx} + \phi_{yy}=0 , \quad -h < y < \eta . \end{equation}

In the absence of surface tension the kinematic and dynamic equations at the free surface are

(2.2)\begin{gather} \eta_t +\phi_x \eta_x -\phi_y= 2\kappa \eta_{xx} - \int_{0}^{t}\left(\frac{\tau_a}{\rho_w }\right)_x \,{\rm d}t, \quad \text{on} \ y=\eta , \end{gather}

(2.3)\begin{gather} \phi_t +\frac{1}{2}[(\phi_x)^2 + (\phi_y)^2] + g\eta ={-} \frac{P_a}{\rho_w}+ 2\kappa \phi_{xx}, \quad \text{on} \ y=\eta, \end{gather}

where $\rho _w$ is the constant water density, $\kappa$ is the constant kinematic water viscosity, $P_a(x,t)$ is the surface air pressure and $\tau _a$ is the air tangential stress at the free surface. From Dutykh & Dias (Reference Dutykh and Dias2007) the bottom boundary condition is

(2.4)\begin{equation} \phi_y = \left(\frac{\kappa }{{\rm \pi} }\right)^{1/2} \int^{t}_{0} \frac{\phi_{xx}(x,y={-}h, t-\sigma)}{\sigma^{1/2}}\,{\rm d}\sigma, \quad \text{on} \ y ={-}h . \end{equation}

Since friction is assumed to be small it is sufficient to assume $\psi$ satisfies a linear Navier–Stokes equation:

(2.5)\begin{equation} \psi_t = \kappa (\psi_{xx} + \psi_{yy}) . \end{equation}

The bottom boundary condition (2.4) arises from the requirement that $u=0, v=0$ at $y=-h$ and that using a thin boundary layer analysis at $y=-h$,

(2.6)\begin{equation} \psi = \left(\frac{\kappa }{{\rm \pi} }\right)^{1/2} \int^{t}_{0} \frac{\phi_{x}(x,y={-}h, t-\sigma)}{\sigma^{1/2}}\,{\rm d}\sigma ,\quad \text{on} \ y ={-}h . \end{equation}

In the following we consider a small-amplitude periodic wave given by

(2.7)\begin{equation} \eta \approx A\exp{({\rm i}kx -{\rm i}\omega t)} + \text{c.c.} +\cdots , \end{equation}

where $\omega$ and $k$ are the frequency and wavenumber of the carrier wave, which at leading order, in the absence of friction and forcing, are both real-valued and satisfy the linear dispersion relation

(2.8)\begin{equation} \omega^2 = gk\tanh{(kh)} . \end{equation}

The envelope amplitude $A(x,t)$ is a slowly varying function of $x, t$ and the omitted terms $[\cdots ]$ in (2.7) are nonlinear and weak dispersive corrections.

The equation set (2.1)–(2.6) and the small-amplitude periodic wave approximation ((2.7) and (2.8)) are expressed in non-dimensional form using time and length scales $S_T, S_L$ to remove the physical parameters $g, h$. We choose $S_T = \varOmega ^{-1}, S_L = K^{-1}$ where $\varOmega$ and $K$ are a characteristic frequency and wavenumber of the carrier wave of a wave group and are subject to the constraint

(2.9)\begin{equation} \varOmega^2 = gK . \end{equation}

This is the free-surface dispersion relation (2.8) when the depth is infinite, chosen so that the scales $\varOmega, K$ do not depend on the depth $h$. The dimensional depth $h$ is replaced by the non-dimensional depth $H= Kh$. In the following we choose $\varOmega$ as the frequency of a typical carrier wave, with $K$ then determined by (2.9). In practice we usually select a frequency for a $5$ s wave so that $\varOmega = 1.257\,{\rm s}^{-1}, K = 0.161\,{\rm m}^{-1}, \varOmega K^{-1} = g \varOmega ^{-1} = 7.81 \, {\rm m}\,{\rm s}^{-1}$. In detail, if the subscript $d$ denotes the dimensional variable, and $n$ denotes the non-dimensional variable, then $x_n = Kx_d$, $y_n = Ky_d$, $t_n = \varOmega t_d$. Similarly the dimensional frequency and wavenumber $\omega _d, k_d$ have non-dimensional frequency and wavenumber counterparts, $\omega _n = \omega _d /\varOmega,\, k_n = k_d /K$, while the parameter $Q=k_d h = k_nH$ is dimensionless. Our motivation for the choice $S_T = \varOmega ^{-1}, S_L = K^{-1}$ is to ensure that $\omega _n, k_n$ are of the order of unity. There is no obligation or requirement to choose $\omega _n =1$ although that is natural and mostly what we do. The subscripts $d, n$ are subsequently omitted; where needed, dimensional variables are distinguished by having units associated with them.

The non-dimensional equations replacing (2.1)–(2.6) are

(2.10)\begin{gather} \phi_{xx} + \phi_{yy}=0 , \quad -H < y < \eta ,\quad H = Kh , \end{gather}

(2.11)\begin{gather} \eta_t +\phi_x \eta_x -\phi_y= 2\nu \eta_{xx} - \int_{0}^{t} \left(\frac{\tau_a}{\rho_w } \right)_x \,{\rm d}t ,\quad \text{on} \ y=\eta , \end{gather}

(2.12)\begin{gather} \phi_t +\frac{1}{2}[(\phi_x)^2 + (\phi_y)^2] + \eta ={-} \frac{P_a}{\rho_w}+ 2\nu \phi_{xx} ,\quad \text{on} \ y=\eta , \end{gather}

(2.13)\begin{gather} \phi_y=\left(\frac{\nu}{\rm \pi} \right)^{1/2}\int^{t}_{0} \frac{\phi_{xx}(x,y={-}H, t-\sigma)}{\sigma^{1/2}}\,{\rm d}\sigma ,\quad \text{on} \ y ={-}H , \end{gather}

(2.14)\begin{gather} \psi_t = \nu(\psi_{xx} + \psi_{yy}) , \quad \nu= \frac{K^2}{\varOmega}\kappa , \end{gather}

(2.15)\begin{gather} \psi = \left( \frac{\nu}{{\rm \pi} } \right)^{1/2} \int^{t}_{0} \frac{\phi_{x}(x,y={-}H, t-\sigma)}{\sigma^{1/2}}\,{\rm d}\sigma ,\quad \text{on} \ y ={-}H . \end{gather}

The linear dispersion relation (2.8) becomes, in the absence of friction and forcing,

(2.16)\begin{equation} \omega^2 = k\tanh{Q} , \quad Q =kH . \end{equation}

The deep-water limit $H \to \infty$ is where $\tanh {Q} \approx 1$, and requires $Q > 2.5$ when the error is less then $1.5\,\%$. In contrast, the shallow-water limit $H \to 0$ is where $\tanh {Q} \approx Q$, and requires $Q < 0.4$ when the error is less then $2\,\%$. For a given fixed $\omega$, the wavenumber $k$ depends on the non-dimensional depth $H$ through (2.16) and

(2.17a,b)\begin{equation} k_H ={-}\frac{ k^2 \operatorname{sech}^2 (Q) }{ [ \tanh{(Q)} + Q\operatorname{sech}^2 (Q)] } , \quad Q_H = \frac{ k\tanh{(Q)}}{[\tanh{(Q)} + Q\operatorname{sech}^2 (Q)] }. \end{equation}

Thus for a fixed frequency $\omega$, $k$ increases as $H$ decreases, while the phase velocity $c = \omega /k$ decreases. Since $Q$ decreases monotonically as $H$ decreases, the domain $0 < H < \infty$ can be replaced by the domain $0 < Q < \infty$, or even by $0 < \mathcal {T} = \tanh {Q} < 1$, useful as $Q, \mathcal {T}$ are dimensionless.

The envelope amplitude $A(x,t)$ at leading order, in the absence of forcing and friction, propagates with the group velocity $c_g$:

(2.18)\begin{equation} c_g = \omega_k = \frac{c}{2}\bigg(1 + \frac{Q(1- \tanh^2{(Q))}}{\tanh{(Q)}}\bigg) , \quad c = \frac{\omega}{k}. \end{equation}

Wave energy is defined by

(2.19)\begin{equation} \mathcal{E} = \left\langle\frac{1}{2} \left(\int^{\eta}_{{-}H}[\phi_x^2 +\phi_y^2 ] \,{{\rm d} y} + \eta^2 \right) \right\rangle, \end{equation}

where $\langle \cdots \rangle$ denotes a phase average for waves periodic in a phase $\theta =kx-\omega t$ given by (2.7). Note that $\mathcal {E}$ is the inviscid part of the wave energy and to leading order $\mathcal {E} =2|A|^2$, and depends on $(x,y,t)$ through the dependence on $|A|^2$. Then, using (2.10)–(2.13),

(2.20)\begin{align} \mathcal{E}_t & = \left\langle\int^{\eta}_{{-}H} [\phi_x \phi_{xt} +\phi_y \phi_{yt} ] \,{{\rm d} y} + \left[\eta + \frac{1}{2}[\phi_x^2 +\phi_y^2]\right]\eta_t (\kern0.7pt y=\eta) \right\rangle \nonumber\\ & = \left\langle\int^{\eta}_{{-}H}\left[\phi_x \phi_{xt} +\phi_y \phi_{yt} \right] \,{{\rm d} y} + \left[-\phi_t - \frac{P_a}{\rho_w}+ 2\nu \phi_{xx}\right]\eta_t (\kern0.7pt y=\eta) \right\rangle \nonumber\\ & = \left\langle\vphantom{\left[-\phi_t - \frac{P_a}{\rho_w }+ 2\nu \phi_{xx}\right]} \left[\phi_y \phi_{t}\right](\kern0.7pt y=\eta)- [\phi_y \phi_{t}](\kern0.7pt y={-}H) - \left[\phi_t \phi_x \eta_x\right](\kern0.7pt y=\eta)\right. \nonumber\\ &\left.\quad + \left[-\phi_t - \frac{P_a}{\rho_w }+ 2\nu \phi_{xx}\right] \eta_t (\kern0.7pt y=\eta) \right\rangle \nonumber\\ & = \left\langle - [\phi_y \phi_{t}](\kern0.7pt y={-}H) + \left[-\eta_t \frac{P_a}{\rho_w } + 4\nu \eta_t \phi_{xx} + \phi_t \int_{0}^{t} \left(\frac{\tau_a}{\rho_w } \right)_x \,{\rm d}t \right](\kern0.7pt y=\eta) \right\rangle . \end{align}

This can be used to define a growth rate $\varDelta$:

(2.21)\begin{equation} \varDelta = \frac{\mathcal{E}_t}{2 \mathcal{E}} . \end{equation}

Parameter $\varDelta$ is non-zero due to (1) the wind forcing term from $P_a$, (2) the water friction and air wave stress terms at the free surface $y=\eta$ and (3) the bottom stress term at $y=-H$. Correspondingly we write

(2.22)\begin{equation} \varDelta = \varDelta_1 + \varDelta_2 +\varDelta_3 . \end{equation}

More details on the definitions and properties of each $\varDelta _{i}, i= 1,2,3$, are given in the following sections. The derivations using (2.20) follow in §§ 2.2 and 2.3; see (2.39) for $\varDelta _1$, (2.40) and (2.45) for $\varDelta _2$ and (2.32) and (2.35) for $\varDelta _3$, expressed in non-dimensional variables. The outcomes are summarised in (4.5) in § 4 expressed there in dimensional variables.

2.2. Wave stress

In this subsection we examine the stress boundary conditions at the free surface and at the bottom and their impact on the wave field. The continuity of normal stress at the free surface yields at leading order

(2.23)\begin{equation} p = P_a - 2\rho_w \nu \phi_{xx} ,\quad \text{on} \ y = \eta , \end{equation}

where $p$ is the water pressure. In this modified potential flow development, the expression (2.23) yields the boundary condition (2.12). The tangential stress is defined by the tensor component $\tau$ where

(2.24)\begin{equation} \left. \begin{aligned} & \frac{\tau }{\rho_w } = \nu (u_y + v_x ) = \nu(2\phi_{xy} + \psi_{yy} - \psi_{xx} )\\ \text{or also} \quad & \frac{\tau }{\rho_w } = \nu (2v_{x} + \psi_{yy} + \psi_{xx} ) = 2\nu v_{x} + \psi_t . \end{aligned} \right\} \end{equation}

The shear stress at the free surface is denoted by $\tau _s$ and satisfies the boundary condition

(2.25)\begin{equation} \tau_s = \tau_a ,\quad \text{on} \ y = \eta , \end{equation}

where $\tau _a$ is the shear stress in the air. The boundary condition (2.25) implies that

(2.26)\begin{equation} \psi ={-} 2\nu \eta_x + \int_{0}^{t}\frac{\tau_a}{\rho_w } \,{\rm d}t ,\quad \text{on} \ y = \eta . \end{equation}

The free-surface boundary layer in the air is usually assumed to be turbulent due to the wind action, with the parametrisations

(2.27)\begin{equation} \tau_a = \rho_a c_d |W_a|W_a , \end{equation}

where $\rho _a$ is the air density, $W_a$ is a measure of the wind speed near the surface and $c_d$ is a dimensionless drag coefficient which increases as $|W_a |$ increases (see e.g. Tang et al. (Reference Tang, Grimshaw, Sanderson and Holland1996) amongst many papers from the storm surge literature).

It is useful to note here that conservation of mass is expressed by

(2.28)\begin{equation} \eta _t + \left(\int_{{-}H}^{\eta} \phi_x \,{{\rm d} y} \right)_x = 2\nu \eta_{xx} - \int_{0}^{t} \left(\frac{\tau_a}{\rho_w } \right)_x \,{\rm d}t + \int_{0}^{t} \left(\frac{\tau_b}{\rho_{w}} \right)_x \,{\rm d}t, \end{equation}

where $\tau _b$ is the bottom stress (see (2.29) and (2.30) in the next paragraph).

The bottom shear stress is denoted by $\tau _b$. From (2.24), since at the bottom $y=-H, u=0, v=0$, we get that

(2.29)\begin{equation} \frac{\tau_b }{\rho_w } = \psi_t (\kern0.7pt y ={-} H) = \left(\frac{\nu }{{\rm \pi} }\right)^{1/2}\int^{t}_{0} \frac{\phi_{xt}(x,y ={-}H,t-\sigma )}{\sigma^{1/2}}\,{\rm d}\sigma, \end{equation}

on using the expression (2.15) for $\psi$ at $y=-H$. Again, since at the bottom $y=-H$, $u=0, v=0$, where the latter implies that $\phi _y = \psi _x$. Thus, using (2.29), the bottom boundary condition (2.13) can be written as

(2.30)\begin{equation} \phi_y=\int_{0}^{t}\left(\frac{\tau_b}{\rho_{w}}\right)_x \, {\rm d}t ,\quad \text{on} \ y ={-}H . \end{equation}

For a plane small-amplitude periodic wave given by (2.7) the bottom boundary condition (2.13) becomes, assuming $\omega > 0$ without loss of generality,

(2.31) \begin{equation} \left. \begin{gathered} \phi_y =\left(\frac{\nu}{{\rm \pi} \omega}\right)^{1/2} \left[\frac{{\rm i}k^2 A }{\omega\cosh{(kH)}}\exp{({\rm i}kx-{\rm i}\omega t) I(t)} +\text{c.c.} \right],\quad \text{on} \ y ={-}H ,\\ I(t) = \int_{0}^{\omega t} \frac{\exp{({\rm i} \varTheta)}}{\varTheta^{1/2}}\,{\rm d}\varTheta \quad\to \left( \frac{{\rm \pi} }{2} \right)^{1/2}(1+{\rm i}) \quad \text{as} \ \omega t \to \infty. \end{gathered} \right\} \end{equation}

Bottom friction leads to loss of wave energy; for small-amplitude slowly varying periodic waves expressed by (2.7) this is given by (2.19) and (2.20) with $\mathcal {E}=2|A|^2$. Using the expression (2.31) in (2.20) and keeping only the term for the bottom stress, we define the corresponding $\varDelta _{3[s]}$ by

(2.32) \begin{align} \left. \begin{gathered} \mathcal{E}_t = \left\langle - \phi_y \phi_t (\kern0.7pt y={-}H) \right\rangle ,\text{ where} \\ \phi_y (\kern0.7pt y={-}H) = \left(\frac{\nu}{ {\rm \pi} \omega}\right)^{1/2} \left[\frac{{\rm i}k^2 A}{ \omega \cosh{(kH)}}\exp{({\rm i}kx-{\rm i}\omega t) I(t)} +\text{c.c.} \right],\\ \phi_t (\kern0.7pt y={-}H) ={-} \left[\frac{ A }{\cosh{(kH)}}\exp{({\rm i}kx-{\rm i}\omega t)} +\text{c.c.} \right],\\ \varDelta_{3[s]} =\frac{\mathcal{E}_t}{2\mathcal{E}} ={-} \left(\frac{\nu}{{\rm \pi} \omega}\right)^{1/2} \frac{k^2}{2\omega}\frac{1}{\cosh^2{(kH)}}J(t),\\ J(t) = \text{Im} \,I(t) = \int_{0}^{\omega t} \varTheta^{{-}1/2}\sin{(\varTheta)}\,{\rm d}\varTheta \to \left( \frac{{\rm \pi} }{2}\right)^{1/2}, \quad \text{as} \ \omega t \to \infty. \end{gathered} \right\} \end{align}

In order to take account of bottom turbulence (2.29) is replaced by (see e.g. Hasselmann & Collins Reference Hasselmann and Collins1968)

(2.33)\begin{equation} \tau_{b[t]} = \rho_w C_D |\phi_x | \phi_x ,\quad \text{on} \ y ={-}H , \end{equation}

where $C_D$ is a dimensionless drag coefficient (often $C_D = 0.0015$ is constant). In this case, the bottom boundary condition is (2.30) with $\tau _b$ now given by $\tau _{b[t]}$. Here, we note again that at the bottom $y= - H , v =0,\phi _y = \psi _x$ and so

(2.34)\begin{equation} \mathcal{E}_{t} (\kern0.7pt y={-}H) = \langle - [\phi_y \phi_{t}] \rangle = \left\langle - [\psi_x \phi_{t}] \right\rangle = \left\langle - [\psi_t \phi_{x}] \right\rangle, \quad \text{where}\ \psi_t = \frac{\tau_b }{\rho_w } . \end{equation}

To proceed, the expression in (2.32) is replaced by an expression for $\varDelta _{3[t]}$:

(2.35) \begin{equation} \left. \begin{aligned} & 2\mathcal{E} \varDelta_{3[t]} = \mathcal{E}_t ={-} \Big\langle \Big[\phi_x \frac{\tau_{b[t]}}{\rho_w}\Big](\kern0.7pt y={-}H) \Big\rangle\\ & \qquad\quad\, ={-} \left\langle \left[\phi_x C_D |\phi_x | \phi_x \right](\kern0.7pt y={-}H) \right\rangle , \\ & \phi_x (\kern0.7pt y={-}H) = \Big[\frac{\omega A }{\sinh{(H)}}\exp{({\rm i}kx-{\rm i}\omega t)} +\text{c.c.} \Big],\\ &\varDelta_{3[t]} ={-} \frac{16 C_D \omega k|A| }{ 3 \cosh{kH} \sinh^2{(kH)}}, \end{aligned} \right\} \end{equation}

which gives a connection between $\nu$ and $C_D$.

2.3. Wind forcing

In general an analogous set of equations is required for the air flow, with coupling to the water flow at the common free surface. Rather than deal with a complicated two-fluid system, we follow the pioneering approach of Miles (Reference Miles1957) and seek an expression for $P_a (x,t)$ directly in terms of $\eta (x,t)$, as this would then close the water wave modified potential flow equations (2.10)–(2.13). For small-amplitude waves, it is sufficient to use the linearised air-flow equations for the flow about a basic horizontal wind shear profile $W(\kern0.7pt y) > 0$ which vanishes at the air–water interface, $W(0)=0$. Miles (Reference Miles1957) showed that if $\eta$ is given by the small-amplitude periodic wave (2.7) then

(2.36)\begin{equation} \frac{P_{a}}{\rho_a} = (\alpha k\eta+ \beta \eta_x) W_{r}^2 = (\alpha + {\rm i} \beta ) kW_{r}^2 A\exp{({\rm i}kx -{\rm i}\omega t)} + \text{c.c.}, \end{equation}

where $\alpha, \beta$ are dimensionless parameters determined from the air-flow solution (see § 3), $\rho _a$ is the air density and $W_r$ is an appropriate reference velocity measuring the wind speed near the surface. Following Miles (Reference Miles1957) and many subsequent works, one choice is that $W_r = u_{*} \varkappa ^{-1} \ \textrm {m}\ \textrm {s}^{-1}$, where $u_{*} \, \textrm {m}\ \textrm {s} ^{-1}$ is the friction velocity for wind over water and $\varkappa$ is the von Kármán constant. The friction velocity $u_{*}$ is determined empirically, and has typical values ranging from about $0.05 \ \textrm {m}\ \textrm {s}^{-1}$ in light winds to about $1.0 \ \textrm {m} \ \textrm {s}^{-1}$ in strong winds. We will choose $u_{*} = 0.36 \ \textrm {m}\ \textrm {s}^{-1}$ so that the reference velocity $W_{r} = 0.9 \ \textrm {m} \ \textrm {s}^{-1}$. Estimates for the parameters $\alpha, \beta$ were given by Miles (Reference Miles1957) and in many subsequent papers by him and others, and depend on the wind shear profile. Using a logarithmic wind shear profile (see § 3.2), Miles (Reference Miles1957) initially estimated that $\beta \approx 10$ at its maximum value when $c \sim 4 W_{r}$, but based on more accurate numerical calculations this was later revised to $\beta \approx 3$ (Miles Reference Miles1959). We re-examine this again in § 3.

As in Miles (Reference Miles1957), Janssen (Reference Janssen2004), Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010) and many other works, it is useful to expand the dispersion relation (2.16) to take account of the weak frictional and forcing effects. The outcome is

(2.37)\begin{align} \mathcal{D}(\omega, k)\eta \equiv \{\omega^2 - k\tanh{(kH)}\}\eta = k\tanh{(kH)} \left(\frac{P_a}{\rho_w}\right) + 4\nu k^2 \eta_{t} + \left(\frac{\tau_a}{\rho_w } \right)_x - \operatorname{sech}(kH) \left(\frac{\tau_b}{\rho_{w}} \right)_x . \end{align}

The wave stress $\tau _a$ at the free surface is given by the turbulent expression (2.27), and at the bottom the wave stress $\tau _b$ is given by the laminar expression (2.29) or by the turbulent expression (2.33). Here we note that the term proportional to $\beta$ in (2.36) for $P_a$ implies by itself that the frequency must be complex-valued, $\omega = \omega _r +\textrm {i}\omega _i$ where $\omega _i > 0$ implies instability. Since $\rho _a/\rho _w \ll 1$, it follows that $\omega _i / \omega _r \ll 1$ and at leading order $\omega _r , k$ satisfy the dispersion relation (2.16). Considering just $P_a$, substituting from (2.36) and expanding $\mathcal {D}$ for small $\omega _i$ leads to

(2.38)\begin{equation} \frac{\omega_i }{\omega_r } = \frac{\rho_a}{\rho_w}\frac{1}{2}\beta k W_{r}^2 . \end{equation}

Similar expressions generating $\omega _i \ne 0$ can be found from $\tau _a, \tau _b$ and are described later. In (2.38) only the contribution from $P_a$ which is in phase with the free-surface slope is needed (see Kharif et al. Reference Kharif, Kraenkel, Manna and Thomas2010), showing that $\beta > 0$ for instability. Here we have assumed for simplicity that $k$ is real-valued, and that $k, \omega _r, c_r = \omega _r / k > 0$. But we note that for wave groups the group velocity $c_g =\omega _k$ must be real-valued, implying that the wavenumber $k$ should also be complex-valued with $\omega _i \approx c_g k_i$ (see Grimshaw Reference Grimshaw2018, Reference Grimshaw2019b). However, for simplicity and in accordance with the usual convention we assume here that $k$ is real-valued, with both $\omega _r, k > 0$. If we omit the term in $\alpha$ in (2.36) then the system (2.10)–(2.13) is closed with the pressure condition $P_{a} = \rho _a W_{r}^2 \beta \eta _x$ independently of the periodic wave assumption (2.7). The growth rate $\varDelta _1$ arises from this modified pressure term (2.36). Using (2.20)–(2.22), or more directly from (2.38) it is expressed here in non-dimensional variables:

(2.39)\begin{equation} \varDelta_1 = \frac{\rho_a}{2\rho_w }\beta \omega_r k W_{r}^2 , \end{equation}

which agrees with that in Maleewong & Grimshaw (Reference Maleewong and Grimshaw2022a,Reference Maleewong and Grimshawb). Note that the critical layer instability mechanism of Miles (Reference Miles1957) requires that $c_{r} =\omega _{r} /k > 0$ and so $\varDelta _1 > 0$ in (2.39) as required.

Here we add to that the friction terms with coefficient $\nu$ and wind stress driving term due to $\tau _a$; see (2.25) and (2.27). For the first part, it is sufficient to use the linearised equations and replace $2\nu \eta _{xx}$, $2\nu \phi _{xx}$ in (2.11) and (2.12) with $-2k^2 \nu \eta$, $-2k^2 \nu \phi$, respectively. In effect $-\textrm {i}\omega$ is replaced with $-\textrm {i}\omega + 2k^2\nu$ and so $A_t$ in (4.2) (see § 4) is replaced with $A_t + 2k^2\nu A$. This also follows from the expanded dispersion relation (2.37) where the last term on the right-hand side is $-4\nu \textrm {i}\omega k^2 \eta$. Then expanding the left-hand side as in (2.38) shows that $\omega _i \approx - 2\nu \omega k^2$ as before. The contribution $\varDelta _{2[s]}$ to the growth rate $\varDelta _2$ is found from this expression, or from (2.20), and expressed here in non-dimensional coordinates:

(2.40)\begin{equation} \varDelta_{2[s]} ={-} 2k^2\nu . \end{equation}

This agrees with the expression in Leblanc (Reference Leblanc2007) and Kharif et al. (Reference Kharif, Kraenkel, Manna and Thomas2010). The second part generates a contribution $\varDelta _{2[t]}$, where again using (2.20) and that $\tau _a$ is given by (2.27):

(2.41)\begin{align} 2\mathcal{E} \varDelta_{2[t]} & = \left\langle\left[\phi_t \int_{0}^{t}\left(\frac{\tau_a}{\rho_w }\right)_x \,{\rm d}t\right](\kern0.7pt y=\eta) \right\rangle\nonumber\\ & = \left\langle\left[\phi_t \int_{0}^{t}\frac{\rho_a}{\rho_w}(c_d |W_a| W_a)_x \,{\rm d}t \right](\kern0.7pt y=\eta) \right\rangle\nonumber\\ & = \left\langle\left[\phi_x \frac{\rho_a}{\rho_w}c_d |W_a| W_a \right] (\kern0.7pt y=\eta) \right\rangle . \end{align}

To estimate the near-surface wind $W_a$ we put $W_a \approx u_a$, the horizontal air-flow velocity at the free surface found from the linearised air flow equations (see (3.1) in § 3.1). As in Miles (Reference Miles1957), we use these to show that $\rho _a c u_a = P_{a}$, where $P_a$ is given by (2.36), so that

(2.42)\begin{equation} u_a = \frac{\alpha W_{r}^2}{c} (\alpha k \eta +\beta \eta_x ). \end{equation}

At the free surface $c\phi _x = \eta$ to leading order, see (2.12), and so from the expression (2.41) we see that the maximum contribution to $\varDelta _{2[t]}$ will come from the terms in $u_{a}$ proportional to $\eta$. Hence as an estimate we set

(2.43)\begin{equation} u_a = U_{a} k \eta, \quad \text{where} \ U_{a} = \frac{\alpha W_{r}^2}{c}. \end{equation}

As $\alpha$, like $\beta$, is not known at present, but both are $O(1)$ dimensionless parameters, we maximise this contribution by choosing $\alpha = \alpha _M$ the maximum positive value found in our long-wave approximation described in § 3.1. Thus finally (2.27) becomes

(2.44)\begin{equation} \frac{\tau_a }{ \rho_a} = c_d U_{a}^2 k^2 |\eta| \eta , \end{equation}

and then (2.41) is given by

(2.45)\begin{equation} \varDelta_{2[t]} = \frac{\rho_a}{\rho_w} c_d U_{a}^2 \frac{16k^2|A|}{3c } , \end{equation}

where $c_d$ is evaluated at a representative value $c_d = 0.0028$ (see Tang et al. Reference Tang, Grimshaw, Sanderson and Holland1996). We note here that Branger et al. (Reference Branger, Manna, Luneau, Abid and Kharif2022) showed from several laboratory experiments that $c_d$ had a dependence on the depth for very shallow water, $kh <1.8$. We attribute this to the higher and steeper waves that form in shallow water and take account of this through the dependence of $W_a =u_a$ on the wave amplitude that we have developed here in (2.43). Note that $\varDelta _2 = \varDelta _{2[s]} + \varDelta _{2[t]}$ is composed of two contributions. The first $\varDelta _{2[s]}$ (2.40) is due to frictional boundary layer in the water and is negative. The second $\varDelta _{2[t]}$ (2.45) is due to wave stress in the air boundary layer at the free surface and is positive, increasing with wave amplitude $A$. The sum $\varDelta _2$ can be positive or negative.

3.1. Formulation and long-wave approximation

As in Miles (Reference Miles1957) and in many works by him and others, the air flow is described by equations linearised about a wind shear profile $W(\kern0.7pt y)$ which in the inviscid limit leads to the well-known Rayleigh equation. Here we follow the approach used by Grimshaw (Reference Grimshaw2018, Reference Grimshaw2019a). In standard notation for air of constant density $\rho _a$ and for inviscid flow, using our non-dimensional variables the equations are

(3.1)$$\begin{gather} \rho_a(D_W u + v W_{y}) + p_x = 0 , \end{gather}$$

(3.2)$$\begin{gather}\rho_a D_W v + \rho_a + p_y = 0 , \end{gather}$$

(3.3)$$\begin{gather}u_x + v_y = 0 , \end{gather}$$

(3.4)$$\begin{gather}\text{where} \ D_W = \frac{\partial }{\partial t} + W \frac{\partial }{\partial x } . \end{gather}$$

To ensure that the wind generates wind waves in the positive $x$ direction, we assume that $W(\kern0.7pt y) >0, y>0$, with $W(0)=0$, and further that the wind is monotonically increasing with height, $W_y > 0$. Further as $y \to \infty$ we require that $W \to W_0$, a constant. Next we introduce the vertical particle displacement $\zeta$ defined in this linearised formulation by

(3.5)\begin{equation} D_W \zeta = v . \end{equation}

The substitution of $v$ by $\zeta$ and eliminating $u, p$ yields a single equation for $\zeta$:

(3.6)\begin{equation} (D_{W}^2 \zeta_y )_y + (D_{W}^2 \zeta)_{xx} = 0 . \end{equation}

This equation is supplemented with the boundary conditions that as $y \to \infty$, $\zeta \to 0$ and as $y \to 0$, $\zeta \to \eta$ and $p \to P_a$.

Since our concern is with wave groups, we seek an asymptotic solution of (3.6) in the form, analogous to (2.7) for $\eta$,

(3.7)\begin{equation} \zeta = A(x,t)\varphi(\kern0.7pt y)\exp{({\rm i}kx -{\rm i}\omega t)} + \text{c.c.} \end{equation}

Here $A(x,t)$ is a slowly varying wave packet amplitude of a carrier wave with frequency $\omega = kc$, wavenumber $k$ and phase velocity $c$. Substitution of (3.7) into (3.6) yields at leading order the modal equation for the modal function $\varphi (\kern0.7pt y)$:

(3.8)\begin{equation} ((W-c)^2 \varphi_y)_y - k^2 (W-c)^2 \varphi = 0 . \end{equation}

Expressed in terms of $v$ from (3.5) this is the well-known Rayleigh equation for a shear flow. The modal function $\varphi$ is normalised by $\varphi (\kern0.7pt y =0) = 1$ so that $\zeta (\kern0.7pt y=0) = \eta$. Then the required relation between the surface pressure $P_a$ and the surface displacement $\eta$ is given by

(3.9) \begin{equation} \left. \begin{gathered} P_a = \rho_a c^2 \varphi_{y}(0) A(x, t)\exp{({\rm i}kx -{\rm i}\omega t)} + \text{c.c.},\\ \eta = A(x,t) \exp{({\rm i}kx -{\rm i}\omega t)} + \text{c.c.} \end{gathered} \right\} \end{equation}

Within this linear approximation elimination of $A(x,t) \exp {(\textrm {i}kx -\textrm {i}\omega t)}$ yields the desired relation (2.36) in terms of $\varphi _{y}(0)$:

(3.10)\begin{equation} (\alpha + {\rm i} \beta ) k W_{r}^2 = c^2 \varphi_{y}(0) ={-} \int_{0}^{\infty} k^2 (W-c)^2 \varphi \, {{\rm d} y} . \end{equation}

Here the alternative integral expression, obtained by integration of (3.8), is related to that used by Miles (Reference Miles1957) but will not be used here. The parameters $\alpha, \beta$ can now be obtained in terms of the available physical parameters. Importantly the frequency $\omega$ and wavenumber $k$ in the modal (3.8) must be complex-valued with albeit small imaginary parts $\omega _i, k_i$ while keeping the group velocity $\omega _k$ real-valued (see Grimshaw Reference Grimshaw2018, Reference Grimshaw2019b). But here for simplicity we take the conventional view that $k$ is real-valued.

The modal equation only rarely has analytic solutions for the considered realistic wind shear profiles $W(\kern0.7pt y)$. We assume that $W(\kern0.7pt y)$ is smooth, vanishes at $y=0, W(0)=0$ and is monotonically increasing with height $y$ in a domain $0 \le y \le y_0$ where $W_{y} > 0$, and in addition we assume that for $y \ge y_0 , W = W(\kern0.7pt y_0) = W_0$ is a constant. In particular there are no exact analytic solutions for the commonly used logarithmic wind shear profile of § 3.2 (see Miles Reference Miles1957; Janssen Reference Janssen2004). Thus the modal equation requires a detailed asymptotic analysis to yield explicit expressions for $\alpha, \beta$. To achieve this, it is customary to take the limit $\rho _a / \rho _w \to 0$, $\rho _a$ and $\rho _w$ being the air and water densities, when $c_i = \textrm {Im}\, c\to 0$. Then various approximations have been used which generally require evaluation of the modal function near a critical level $y_c$ where $W(\kern0.7pt y_c ) = c_r = \textrm {Re}\, c$ and there is a singularity in the modal equation (3.8). Indeed a critical-level singularity causes wind wave generation, first shown by Miles (Reference Miles1957) in his pioneering seminal paper.

Before proceeding with our asymptotic analysis, it is useful to note the invariant

(3.11)\begin{equation} \mathcal{N} \equiv \int_{0}^{\infty} (W - c)^2 [\varphi^2_y +k^2 \varphi^2] \,{{\rm d} y} + c^2 \varphi_{y}(0)= 0 . \end{equation}

This is an expanded dispersion relation, equivalent to (2.37) when only the forcing term due to $P_a$ is included. Next, in $y \ge y_0$ where $W =W_0$ is a constant,

(3.12)\begin{equation} \varphi (\kern0.7pt y) = D_0 \exp{({-}k(\kern0.7pt y-y_0))} , \quad y \ge y_0 , \end{equation}

where $D_0$ is the constant amplitude at $y=y_0$, and we assume that $k$ is real-valued and positive. It then follows that

(3.13)\begin{equation} \mathcal{N} \equiv \int_{0}^{y_0}\, (W - c)^2 [\varphi^2_y +k^2 \varphi^2] \,{{\rm d} y} + D_{0}^{2} k (W_0 - c)^2 + c^2 \varphi_{y}(0) = 0 . \end{equation}

As mentioned above, the modal equation (3.8) usually does not have analytic solutions and hence approximations are needed to find explicit expressions for $\alpha, \beta$. For the logarithmic profile of § 3.2, Miles (Reference Miles1957) used the simple approximation that $\varphi (\kern0.7pt y) \propto \exp {(-ky)}$ expressed here in our formulation and notation. This is exact for $y \to \infty$ but does not capture the critical-level structure. Miles (Reference Miles1957) addressed this failing by using the integral expression in (3.10) to find an expression for $\beta$, which suggested a maximum value of $\beta \approx 10$. Later this was reduced by Miles (Reference Miles1959) to a maximum value of $\beta \approx 3$ based mainly on a numerical solution. In the literature there are several other approximations that have been used, many based on moving the lower air-flow boundary from $y=0$ to $y= y_{*}$, the roughness length scale, and then analysing the near-surface turbulent air flow. A series of papers (Montalvo et al. Reference Montalvo, Dorignac, Manna, Kharif and Branger2013a,Reference Montalvo, Kraenkel, Manna and Kharifb; Latifi et al. Reference Latifi, Manna, Montalvo and Ruivo2017; Branger et al. Reference Branger, Manna, Luneau, Abid and Kharif2022) used a two-term critical-level Bessel function approximation to represent the modal function, sometimes combined with the lower boundary modification. These and other approximations give similar $O(1)$ estimates for $\beta$.

Here we use a long-wave approximation similar to those used by Janssen (Reference Janssen2004) and Grimshaw (Reference Grimshaw2018, Reference Grimshaw2019a). We assume that in the domain $0 < y< y_0$ the term $(W-c)^2 k^2 \varphi$ in (3.8) is neglected, formally valid when $|ky| \ll 1$, and in particular, $|ky_c | \ll 1$. Then an approximate asymptotic solution is

(3.14)\begin{equation} \varphi \approx 1 + D\int^{y}_{0}\frac{{{\rm d} y} }{(W(\kern0.7pt y)-c)^2} , \quad 0 < y < y_0 , \end{equation}

where the constant $D$ is determined by matching with the exact solution (3.12) in $y \ge y_0$. Across $y= y_0$ both $\varphi, \varphi _y$ are continuous, since $\zeta, p_a$ must be continuous. This also follows from integration of the modal equation (3.8) over a small interval about $y = y_0$. The approximation (3.14) is expressed in non-dimensional variables; when put back into dimensional variables the constant $D_n$ in (3.14) is replaced by $D_d = g D_n$. We note here that the constant $D_0$ in (3.12) is dimensionless. Substitution of (3.14) into (3.13) and neglecting the corresponding term proportional to $k^2$ in (3.13) yields

(3.15)\begin{equation} \mathcal{N} \equiv D^2 \int_{0}^{y_0} \frac{{\rm d} y}{(W - c)^2 } + \frac{D^2}{k (W_0 - c)^2} + D = 0 . \end{equation}

Here we have used the continuity of $\varphi _y$ across $y=y_0$ so that $kD_0 (W_{0}-c)^2 =-D$ and noted that $c^2 \varphi _y(0) = D$. The invariant (3.15) can be used to give an expression for $D$. This approximation (3.14) captures the two leading terms at the critical level, is formally valid for all $0 \le y \le y_0$ and matches to the exact solution in $y > y_0$. As such, we believe it is an improvement over similar approximations in the literature.

Since in the domain of interest $W_y > 0$ we can change the integration variable in (3.14) from $y$ to $W$ to yield

(3.16)\begin{equation} \int^{y}_{0}\frac{{{\rm d} y} }{(W(\kern0.7pt y) - c)^2} = \mathcal{M}(W) = \int^{W}_{0}\frac{\mathcal{S}(W)\,{\rm d}W}{(W-c)^2} , \quad \mathcal{S}(W) = \frac{1}{W_{y}} , \quad 0 < W < W_0 , \end{equation}

where $W_y$ is expressed as a function of $W$. Then an integration by parts yields

(3.17) \begin{equation} \left. \begin{gathered} \mathcal{M}(W) ={-}\left[\frac{\mathcal{S}}{(W-c)}\right] _{0}^{W} - I(W) ,\\ I(W) = \int^{W}_{0}\frac{\mathcal{K}\,{\rm d}W}{(W-c)} , \quad \mathcal{K} ={-} \mathcal{S}_W =\frac{W_{yy}}{W_{y}^3} . \end{gathered} \right\} \end{equation}

Note that the term $I(W)$ in (3.17) depends on the curvature expression $\mathcal {K}$. The invariant approximation (3.15) becomes

(3.18)\begin{equation} \mathcal{N} \equiv D^2 \left(\mathcal{M}(W_0) + \frac{1}{k (W_0 - c)^2} \right) + D = 0 . \end{equation}

As noted above, (3.18) yields an expression for $D$:

(3.19)\begin{equation} D \left( \mathcal{M}(W_0) + \frac{1}{k (W_0 - c)^2} \right)={-}1 . \end{equation}

In the limit $c_{i} \to 0$ the integral term $\mathcal {I}$ has a singularity at $W = c_r, y=y_c$ and is evaluated by assuming that $c_i > 0$, and then taking the limit $c_{i} \to 0^{+}$. This yields the Frobenius expansion when $W_0 >W > W_c = c_r$:

(3.20)\begin{equation} I(W) = \mathcal{P}\int^{W}_{0}\frac{\mathcal{K} \,{\rm d}W}{(W-c_r)} + {\rm i} {\rm \pi}\mathcal{K}_c , \quad \mathcal{K}_c = \mathcal{K}(W = W_c =c_r) . \end{equation}

The imaginary term is proportional to the curvature expression $\mathcal {K}_c$ and in the real term $\mathcal {P}\int$ denotes the principal value integral. It is evaluated by putting the curvature expression $\mathcal {K} = \mathcal {K}_c + (\mathcal {K} - \mathcal {K}_c)$ so that

(3.21)\begin{equation} I(W) = \int^{W}_{0}\frac{(\mathcal{K} - \mathcal{K}_c)}{(W-c_r)} \,{\rm d}W + \mathcal{K}_c \ln{\left[\frac{(W - c_r)}{c_r }\right]} + {\rm i} {\rm \pi}\mathcal{K}_c . \end{equation}

In $y > y_0$ where $W =W_0$ is a constant, $\varphi (\kern0.7pt y)$ is given by (3.12). Across $y = y_0$ both $\varphi$ and $\varphi _y$ are continuous and then from (3.16)

(3.22) \begin{equation} \left.\begin{gathered} D_0 = 1 + D\Big(\Big[-\frac{\mathcal{S}}{(W-c)}\Big] _{0}^{W_0} - I(W_0)\Big),\\ -kD_ 0 = \frac{D}{(W_0 -c)^2 }. \end{gathered} \right\} \end{equation}

Elimination of $D_0$ yields the expression for $D$:

(3.23) \begin{equation} D\Big( \Big[ \frac{\mathcal{S}}{(W-c)} \Big]_{0}^{W_0} - \frac{1}{k(W_0 - c)^2}+ I(W_0) \Big) = 1, \end{equation}

where $I(W_0)$ is given by (3.21). This agrees with the expression (3.19) from the invariant $\mathcal {N}$. In the limit $c_{i} \to 0^{+}$, $I(W_0)$ is evaluated from (3.21) and the only imaginary part is the second term proportional to the curvature expression $\mathcal {K}_c$. Thus,

(3.24) \begin{align} \left. \begin{gathered} D = \frac{1}{E + i{\rm \pi} \mathcal{K}_c } , \\ E = \left[\frac{\mathcal{S}}{(W-c_r)}\right]_{0}^{W_0} - \frac{1}{k(W_0 - c_r)^2} +\int^{W_0}_{0}\,\frac{(\mathcal{K} - \mathcal{K}_c)}{(W-c_r)} \,{\rm d}W + \mathcal{K}_c \ln{\left[\frac{(W_0 - c_r)}{c_r }\right]}, \end{gathered} \right\} \end{align}

where $E$ is real-valued. In order to obtain an explicit expression for the integral term in (3.24) we use a simple approximation based on a one-term Taylor series expansion:

(3.25)\begin{align} \int^{W_0}_{0}\frac{(\mathcal{K} - \mathcal{K}_c)}{(W-W_c)} \,{\rm d}W \approx \mathcal{K}_{Wc}W_0 , \quad \mathcal{K}_{Wc} = \mathcal{K}_{W} (W=W_c ), \quad \mathcal{K}_{W} =\frac{\partial \mathcal{K}(W)}{\partial W} = \frac{\mathcal{K}_y}{W_y}. \end{align}

For the typical wind shear profiles that we consider, this is an underestimate, although it is exact for certain algebraic profiles for which $\mathcal {K}_{W}$ is a constant (see § 3.3). Thus the final expression for $E$ that we will use is

(3.26)\begin{align} E &= \left[\frac{\mathcal{S}(W_0)}{(W_0-c_r)} + \frac{\mathcal{S}(0)}{c_r}\right] - \frac{1}{k(W_0 - c_r)^2} \nonumber\\ &\quad + \mathcal{K}_{Wc}W_0 + \mathcal{K}_c \ln{\left[\frac{(W_0 - c_r)}{c_r }\right]}. \end{align}

Before proceeding we note that (3.26) is expressed in non-dimensional variables. In dimensional variables $E_d = E_n g^{-1}$, where the same expression as (3.26) holds for $E_d$.

Since $c^2 \varphi _{y}(0) = D$ from (3.14), substitution into (3.9) and using (3.24) yields expressions for $\alpha, \beta$:

(3.27)\begin{equation} (\alpha + {\rm i} \beta ) k W_{r}^2 = D , \quad \alpha k W_{r}^2 = \frac{E}{E^2 + ({\rm \pi} \mathcal{K}_c )^2} , \quad \beta k W_{r}^2={-} \frac{{\rm \pi} \mathcal{K}_c }{E^2 + ({\rm \pi} \mathcal{K}_c )^2} . \end{equation}

Then since $\beta > 0$ for instability we require that $\mathcal {K}_c < 0$, that is, from (3.17), $W_{yy}(W=W_c =c_r) <0$, as is now very well known from the work of Miles (Reference Miles1957). This is the essential condition for a critical-level instability. Note that with $W_{yy} < 0$, $W_{y} > 0$ decreases as $y$ increases. The expression (3.27) shows that $\alpha, \beta$ for a fixed wave frequency $\omega$ depend on $c_r = \omega k^{-1}$ and the parameters of the wind shear profile, these latter occurring only in $E$. One important parameter of the wind shear profile is the choice of the limiting velocity $W_0$. It is constrained so that $W_0 > W_c$ for all depths to ensure that there is always a critical level. As we discuss further below and in §§ 3.2–3.4 for specific wind shear profiles, $W_0$ is chosen to yield the maximum value of $\beta$. For the wave parameters we select a frequency $\omega$ for a $5$ s wave when $\omega = 1.26 \ \textrm {s}^{-1}$; then the wavenumber $k$ comes from the linear dispersion relation (2.16) where $k$ depends on the depth $H$, (2.17a,b). We recall that then the scaling parameters are $\varOmega = 1.257 \ \textrm {s}^{-1}, K = 0.161 \ \textrm {m}^{-1}, \varOmega K^{-1} = g \varOmega ^{-1} = 7.81 \ \textrm {m}\ \textrm {s}^{-1}$. The critical level $y=y_c$ is where $W(\kern0.7pt y_c)=W_c = c_r$. In deep water for a $5$ s wave, $c_r = 7.81 \ \textrm {m}\ \textrm {s}^{-1}$ and decreases as $H$ decreases. Then we require that $W_0 > 7.81 \ \textrm {m}\ \textrm {s}^{-1}$. In practice we impose a lower bound $W_0 \ge W_{0L} =W(2y_{c})$, where here $y_{c}$ is the critical level in the deep-water limit, $H \to \infty$. For an upper bound we impose $W_0 \le W_{0U} = W(200 y_c)$, where $W_{0U} < W_{\infty } = W(\kern0.7pt y \to \infty )$. These can be relaxed if needed provided that always $c_r < W_{0} < W_{\infty }$.

Since the wind shear profile $W(\kern0.7pt y)$ does not depend on the wave parameters $\omega, k$, or on the depth parameters $H, Q$, the dependence of $E$ (3.26) on these parameters and especially on $H, Q$ is entirely through its dependence on the phase velocity $c_r$ including the terms $\mathcal {K}_c, \mathcal {K}_{Wc}$. Since $c_{r} =\omega k^{-1}$ where $\omega$ is fixed, $c_r$ decreases monotonically as $H, Q$ decrease since $k$ increases monotonically; see (2.17a,b). Based on the explicit results for the specific wind shear profiles, logarithmic, algebraic and exponential, that we consider in §§ 3.2–3.4, we expect $E$ to initially increase as $H$ decreases from infinity and to eventually reach a maximum value when $H$ is quite small. For very small $H$ with a fixed $\omega$, $k \propto H^{-1/2}, c_r \propto H^{1/2}$. The behaviour of $E$ as $H \to 0$ is controlled by two singularities: the first is the term $\mathcal {S}(0)/c_r$ which by itself causes $E$ to increase to positive infinity; the second is the term $-\mathcal {K}_c \ln {c_r }$ which by itself causes $E$ to decrease to negative infinity. Although the second singularity is weaker than the first, we find that initially it dominates as $c_r$ decreases, but at a certain very small value of $c_r$ ($W_r$ for the logarithmic profile) there is a turning point in $E$ after which $E$ increases to positive infinity as $H^{-1/2}$. Correspondingly in this final limit $\alpha, \beta \to 0$ as $H, H^{3/2}$, respectively. But we emphasise that these asymptotic limits occur in practice only for unrealistic small depths.

For fixed wave parameters and depth parameters $H, Q$, the dependence of $\alpha, \beta$ on the shear profile parameters, such as $W_0$, is entirely through $E$ (3.26). In particular through the choice of $W_0$, where the extreme values of $\alpha, \beta$ as $W_0$ varies are found from (3.27):

(3.28)\begin{equation} \frac{\partial \beta }{\partial W_0 }k W_{r}^2 = \frac{2{\rm \pi} \mathcal{K}_c EE_{W_0}}{(E^2 + ({\rm \pi} \mathcal{K}_c )^2 )^3} , \quad \frac{\partial \alpha}{\partial W_0 }k W_{r}^2 = \frac{(({\rm \pi} \mathcal{K}_c )^2 - E^2)E_{W_0}}{(E^2 + ({\rm \pi} \mathcal{K}_c )^2)^3} , \end{equation}

where

(3.29)\begin{align} E_{W_0} = \frac{\partial E }{\partial W_0 } ={-} \frac{\mathcal{S}(W_0)}{(W_0-c_r)^2} - \frac{\mathcal{K}(W_0) }{(W_0 - c_r)} + \frac{2}{k(W_0 - c_r)^3} + \mathcal{K}_{Wc} + \mathcal{K}_c \frac{1}{(W_0 - c_r)}. \end{align}

One set of extreme values is determined by those $W_0$ such that $E_{W_0} =0$ and then $\alpha, \beta$ have simultaneous extreme values. But for typical wind shear profiles and in our parameter range for $W_0$, $E_{W_0} \ne 0$ (see §§ 3.2–3.4 for logarithmic, algebraic and exponential profiles). Otherwise the extreme values of $\alpha, \beta$ are determined by particular values of $E$. Parameter $\beta$ has a single maximum at $E=0$ (where $\alpha =0$):

(3.30)\begin{equation} {[\max]} \beta = \beta_M ={-}\frac{1}{({\rm \pi} \mathcal{K}_c )(k W_{r}^2)} , \end{equation}

and $\alpha$ has maximum, minimum values at $E= \mp ({\rm \pi} \mathcal {K}_c )$ where

(3.31)\begin{equation} {[\max, \min]} \alpha = \alpha_{M, m} ={\mp} \frac{1}{2({\rm \pi} \mathcal{K}_c )(k W_{r}^2)} ={\pm} \frac{\beta_{M}}{2} . \end{equation}

Parameter $\beta _{M}$ is independent of $W_{0}$ as it occurs at a specific value of $W_0 = W_{\beta M}$. But it does depend on the depth $H$, and as $H$ decreases with $\omega$ fixed, $k$ increases, see (2.17a,b), and so $c_r =\omega k^{-1}$ decreases:

(3.32)\begin{equation} \frac{\partial \beta_M }{\partial c_r}= \frac{\mathcal{K}_{Wc}}{({\rm \pi} \mathcal{K}_{c}^2 )(k W_{r}^2)} -\frac{1}{({\rm \pi} \mathcal{K}_c )(\omega W_{r}^2)} = \frac{\mathcal{K}_{Wc} c_r - \mathcal{K}_c}{({\rm \pi} \mathcal{K}_{c}^2 )(\omega W_{r}^2)} . \end{equation}

This is negative for $c_r > c_M = \mathcal {K}_c /\mathcal {K}_{Wc}$ and becomes positive for $c_r < c_M$. For the typical wind shear profiles, logarithmic, algebraic and exponential that we consider in §§ 3.2–3.4, $c_M < c_r (H \to \infty )$ and so $\beta _M$ at first increases as $H, c_r$ decrease until $c_r = c_M$ and then decreases to $0$ as $H \to 0$, where $\beta _M \propto H^{1/2}$. Also for these wind shear profiles $W_{\beta M}$ decreases as $H$ decreases.

3.2. Logarithmic wind shear profile

The most commonly used wind shear profile is the logarithmic wind shear profile, used by Miles (Reference Miles1957, Reference Miles1959), as for turbulent flow near the interface it is supported by theory and experiment. It is given here by

(3.33) \begin{align} \left. \begin{gathered} W (\kern0.7pt y) = W_r \ln{\Big[1 + \frac{y}{y_s}\Big]} , \quad 1 +\frac{y}{y_s}= \exp{(W/W_{r})},\\ W_{y} = \frac{W_{r}}{(\kern0.7pt y_s + y)} , \quad \mathcal{S} = \frac{y_s}{W_r }\exp{(W/W_{r})}, \quad \mathcal{K} ={-}\frac{y_{s}}{W_r^2 } \,\exp{(W/W_{r})}, \quad \mathcal{K}_W = \mathcal{K}/W_{r} . \end{gathered} \right\} \end{align}

Miles (Reference Miles1957) and many others put $W(\kern0.7pt y) =W_r \ln {(\kern0.7pt y/y_s )}$ rather than (3.33), requiring the lower boundary to be moved from $y=0$ to $y=y_{*}$ to avoid the singularity at $y=0$, and then needing analysis of the region $0 < y < y_{*}$ with matching to $y > y_s$ (see the recent work by Abid et al. (Reference Abid, Kharif, Hsu and Chen2022) for instance). The form (3.33), also used by Janssen (Reference Janssen2004), avoids this technical difficulty and has $W(0) =0$ as required. The two expressions differ by less than $1\,\%$ for $y/y_s > 40$. The expression (3.26) for $E$ becomes

(3.34) \begin{align} \left. \begin{aligned} & E = \frac{y_s}{W_r}\left[\frac{ \exp{(W_0 /W_r )}}{(W_0 -c_r)} + \frac{1}{c_r}\right] - \frac{1}{k(W_0 - c_r)^2} + \mathcal{K}_{Wc}W_0 + \mathcal{K}_c\ln{\left[\frac{(W_0 - c_r)}{c_r } \right]} ,\\ &\text{where} \quad \mathcal{K}_c ={-}\frac{y_{s}}{W_r^2 } \exp{(c_r /W_{r})}, \quad \mathcal{K}_{Wc} ={-} \frac{y_s \exp{(c_r/W_r )}}{W_{r}^3 } . \end{aligned} \right\} \end{align}

There are two important parameters: the length scale $y_s$ and the limiting velocity $W_0$. A common choice for $y_s$ is $y_{*}$, the roughness length scale, usually determined empirically but can be estimated as $y_{*} = c_{*} u_{*}^2 g^{-1}\ \textrm {m}$, where the Charnok parameter $c_{*} = 0.015$. Here we set $u_{*} = 0.36 \ \textrm {m}\ \textrm {s}^{-1}$ , and then $y_s = 0.0002\ \textrm {m}$. As this may be too small for our purposes, we also considered another case with an increased value $y_s = 0.001\ \textrm {m}$, based on the ad hoc criterion that $y_s$ is approximately one-tenth of the wave height, where here we consider small-amplitude waves. The results overall were similar, but with a smaller value for $\beta$ by a factor of around $5$, to be expected as then $\mathcal {K}_c , E$ increase by around the ratio of the respective $y_s$ values. In the following we put $y_s = 0.0002\ \textrm {m}$.

We use the parameters of a $5$ s wave, that is, $\varOmega = 1.257 \ \textrm {s}^{-1},$ $K = 0.161 \ \textrm {m}^{-1},$ $\varOmega K^{-1} = 7.811 \ \textrm {m}\ \textrm {s}^{-1}$, so that in deep water we set $\omega = 1.257 \ \textrm {s}^{-1}$ (1); $k = 0.161 \ \textrm {m}^{-1}$ (1); $c_r = 7.811 \ \textrm {m}\ \textrm {s}^{-1}$ (1), where the non-dimensional values are $({\cdot })$. With these parameters, the lower bound $W_{0L} = 8.43\ \textrm {m}\ \textrm {s}^{-1}$ (1.08). The plots of $\alpha, \beta$ as a function of $W_0$ are shown in figure 1 for the deep-water limit $H \to \infty$. In agreement with the preceding discussion around equations (3.30) and (3.31) there is a pronounced maximum value of $\beta = \beta _M = 1.690$ at $W_0 = W_{\beta M} =11.25 \ \textrm {m}\ \textrm {s}^{-1}$ (1.44) where $E = 0$, and for $\alpha$ there is a maximum of $0.845$ at $W_0=W_{\alpha M} =11.56 \ \textrm {m}\ \textrm {s}^{-1}$ (1.48) and a minimum of $-0.845$ at $W_0=W_{\alpha m} =10.81 \ \textrm {m}\ \textrm {s}^{-1}$ (1.40) where $E = \mp ({\rm \pi} \mathcal {K}_c )= \pm 4.554 \ \textrm {m}^{-1} \ \textrm {s}^{2}$, respectively. Hence we set $W_0 = W_{0B} =11.25 \ \textrm {m}\ \textrm {s} ^{-1}$ (1.44) for all depths $H$ which satisfies the constraint $c_r < W_{0L} < W_0 < W_{0U} < w_{\infty }$. In summary we consider a $5$ s wave, with $\omega$ fixed and $k(H)$ found from the linear dispersion relation (2.16) with $y_s = 0.0002 \ \textrm {m}$ (0.000032) and $W_0 = W_{0B} =11.25 \ \textrm {m}\ \textrm {s}^{-1}$ (1.44).

Figure 1. Plots of $\alpha, \beta$ for the logarithmic wind shear profile (3.33), with $y_s=0.0002 \ \textrm {m}$, for a $5$ s wave in deep water. (a) Plot of $\beta$ versus $W_0$; maximum of $\beta = 1.690$ at $W_0=11.252 \ \textrm {m}\ \textrm {s}^{-1}$. (b) Plot of $\alpha$ versus $W_0$; maximum of $\alpha = 0.845$ for $W_0=11.563 \ \textrm {m}\ \textrm {s}^{-1}$ and minimum of $\alpha = -0.845$ for $W_0 = 10.805 \ \textrm {m}\ \textrm {s}^{-1}$.

The outcomes for $E, \alpha, \beta$ as functions of $H$ are shown in figures 2 and 3 for several $W_0$ around the benchmark value $W_0 =W_{0B} = 11.25 \ \textrm {m}\ \textrm {s} ^{-1}$. As expected from the preceding discussion, in all cases $E$ is constant in the deep-water limit $H \to \infty$, zero for the benchmark $W_{0B}$, positive for $W_{0}> W_{0B}$ and negative for $W_{0} < W_{0B}$. Then as $Q, H, c_r$ decrease $E$ increases to a positive maximum and then decreases. The final stage when $E$ again increases occurs for a very small $H$, too small to show on these plots. For this logarithmic wind shear profile this last turning point is less than $c_M = \mathcal {K}_c /\mathcal {K}_{Wc} = W_{r} =0.9 \ \textrm {m}\ \textrm {s}^{-1}$ which in turn is less than the smallest $c_r =1.39 \ \textrm {m}\ \textrm {s}^{-1}$ shown in figure 2. This smallest value shown is for $Q= 0.18$ when the depth $h = 0.2 \ \textrm {m}$, $H = 0.032$; this is physically unrealistic for the nearly frictionless model we are using. Both $\alpha$ and $\beta$ approach constant values in the deep-water limit $H \to \infty$. For $W_0 =W_{0B} = 11.25\ \textrm {m}\ \textrm {s}^{-1}$ when $E=0$ in the deep-water limit, $\beta$ reaches a maximum $\beta _M = 1.745$ at $H=2.769$, while $\alpha$ has a maximum $\alpha _M = 1.176$ at $H=1.738$ and a minimum $\alpha _m = -0.027$ at $H=5.650$ as predicted; see (3.30) and (3.31). For $W_0 > W_{0B}$, $E > 0$, $\beta$ decreases as the depth decreases, while $\alpha$ at first increases with a barely discernible maximum before decreasing. Both $\alpha, \beta \to 0$ as $H \to 0$. When $W_0 < W_{0B}$, $E < 0$, in the deep-water limit $\beta$ tends to a positive value less than for $W_0 = W_{0B}$, while $\alpha$ tends to a negative value. Interestingly, however, there are now marked increases in magnitude to a maximum for both $\alpha, \beta$, more marked the further $W_0$ is decreased.

Figure 2. Plot of $E$ as a function of $H$ as $W_0$ is varied around $W_{0B} =11.25 \ \textrm {m}\ \textrm {s}^{-1}$ for the logarithmic wind shear profile (3.33), with $y_s=0.0002 \ \textrm {m}$, for a $5$ s wave.

Figure 3. Plots of $\alpha, \beta$ as functions of $H$ for the logarithmic wind shear profile (3.33) with $y_s=0.0002 \ \textrm {m}$, for a $5$ s wave, as $W_0$ is varied around $W_{0B} =11.25 \ \textrm {m}\ \textrm {s}^{-1}$: (a) $\beta$ versus $H$; (b) $\alpha$ versus $H$.

We will use $\beta = 1.745$ hereafter as it is $O(1)$ and consistent with the many estimates in the literature. Miles (Reference Miles1959) found a maximum for $\beta$ of around $3$ from a numerical solution, but that was for a smaller value of $c_r$. Using the same $c_r$ as here, the value of $\beta$ in Miles (Reference Miles1959) is reduced by about one-third, lending confidence to our approximation for $\beta$. Nevertheless, there is a substantial spread in the literature for $\beta$ depending on the wind shear profile, the wave parameters and the approximation method.

3.3. Algebraic wind shear profile

Next we consider the algebraic profile, chosen here because it can eliminate the necessity for the approximation (3.25) of the integral term in (3.24). These are a family of profiles:

(3.35) \begin{align} \left. \begin{gathered} W (\kern0.7pt y) = W_r \bigg[\bigg(1 +\frac{y}{y_s} \bigg)^{1/n} - 1\bigg] , \quad 1 +\frac{y}{y_s} = \left[\frac{(W + W_{r})}{W_r }\right]^n ,\\ W_{y} = \frac{W_{r}}{ny_s}\left(1 + \frac{y}{y_s}\right)^{{-}1 +1/n } , \quad \mathcal{S} = \frac{ny_s}{W_r}\left[\frac{(W + W_{r})}{W_r }\right]^{n-1},\\ \mathcal{K} ={-}\frac{n(n-1)y_s}{W_{r}^2}\left[\frac{(W + W_{r})}{W_r }\right]^{n-2}, \quad \mathcal{K}_{W} ={-}\frac{n(n-1)(n-2)y_s}{W_{r}^3}\left[\frac{(W + W_{r})}{W_r }\right]^{n-3}. \end{gathered} \right\} \end{align}

Here the index $n = 2,3,\ldots$ is an integer, and again $y_s$ is a suitable length scale. The expression (3.26) for $E$ becomes

(3.36)\begin{equation} E = \frac{ \mathcal{S}(W_0)}{(W_0 -c_r)} + \frac{ \mathcal{S}(0)}{c_r} - \frac{1}{k(W_0 - c_r)^2} + \mathcal{K}_{Wc}W_0 + \mathcal{K}_c\ln{\left[\frac{(W_0 - c_r)}{c_r } \right]} . \end{equation}

The main cases of interest here are $n=2, 3, 7$:

(3.37)\begin{align} \left. \begin{aligned} n & =2: \quad \mathcal{S} =\frac{2y_s(W + W_{r})}{W_{r}^2}, \quad \mathcal{K} ={-} \frac{2y_s }{W_{r}^2} , \quad\mathcal{K}_W =0,\\ n & =3:\quad \mathcal{S} =\frac{3y_s(W + W_{r})^2}{W_{r}^3}, \quad \mathcal{K} ={-}\frac{6y_s (W + W_r )}{ W_{r}^3 }, \quad \mathcal{K}_W ={-}\frac{6y_s }{ W_{r}^3 },\\ n & = 7:\quad \mathcal{S} =\frac{7y_s(W + W_{r})^6}{W_{r}^7}, \quad \mathcal{K} ={-}\frac{42y_s (W + W_r )^5}{ W_{r}^7 }, \quad \mathcal{K}_W ={-}\frac{210y_s (W + W_r)^4}{W_{r}^7}. \end{aligned} \right\} \end{align}

In the first two cases, $n =2,3$, $\mathcal {K}_W$ is a constant, and so the integral approximation (3.25) is exact. The case $n=7$ has been used as a representation for very-near-surface turbulent wind (see e.g. Counihan Reference Counihan1975; Hsu, Meindl & Gilhousen Reference Hsu, Meindl and Gilhousen1994), though not nearly as commonly as the logarithmic wind shear profile. We again set $W_r = 0.9\ \textrm {m}\ \textrm {s}^{-1}$ but the value of $y_s$ is more difficult to determine here. From (3.30) $\beta _M = -1/({\rm \pi} \mathcal {K}_c )(k W_{r}^2)$ is proportional to $y_s^{-1}$ and if $y_s$ is too small, $\beta$ is too large. To achieve some consistency with the logarithmic profile of § 3.2 here we choose $y_s$ so that in the deep-water limit $\beta _M$ has a comparable value to that for the logarithmic profile. Hence from (3.33) a factor proportional to $\exp {(c_r/W_{r})}$ is needed to adjust $y_s$ and so here we set $y_s = 1 \ \textrm {m}$. Plots of $\alpha, \beta$ versus $W_0$ are shown in figure 4 for $n=2$, from which we choose the benchmark value $W_{0} = W _{0B} = 39.17 \ \textrm {m}\ \textrm {s}^{-1}$. As expected in the deep-water limit $\beta$ has a maximum of $\beta _M = 0.99$ at $E=0$, and $\alpha$ has maximum and minimum values of $\pm \beta _M/2$ at $E = \mp {\rm \pi}\mathcal {K}_c = \pm 7.76$, $W_0 = 12.97 \ \textrm {m}\ \textrm {s}^{-1}, 567.57 \ \textrm {m}\ \textrm {s}^{-1}$; these values are not seen in figure 4 as the relationship between $E$ and $W_0$ is highly nonlinear, but are confirmed in a more detailed plot with larger range of $W_0$ (not shown here). The case $n=3$ is not shown here, but is similar with much smaller values for $\alpha, \beta$, for instance $\beta _{M} = 0.034$. The case $n=7$ produced values for $\alpha, \beta$ much smaller even than the case $n=3$, and so these are not shown here.

Figure 4. Plots of $\alpha, \beta$ for the algebraic wind shear profile (3.35) when $n=2$, with $y_s=1 \ \textrm {m}$, for a $5$ s wave in deep water. (a) Plot of $\beta$ versus $W_0$; maximum of $\beta = 0.99$ at $W_0=39.17 \ \textrm {m}\ \textrm {s}^{-1}$. (b) Plot of $\alpha$ versus $W_0$; maximum/minimum of $\alpha = \pm 0.49$ for $W_0=12.97, 567.57\ \textrm {m}\ \textrm {s}^{-1}$ (not shown here as out of range).

3.4. Exponential wind shear profile

Our third choice is an exponential profile, similar to, but simpler than, a hyperbolic tangent profile. It was shown by Miles in the appendix to the paper by Morland & Saffman (Reference Morland and Saffman1993) that the modal equation (3.8) could be solved exactly in terms of a hypergeometric function, leading to an exact analytic expression for $\beta$. But we will not use that here as it is quite complicated to express, and our purpose is to evaluate the approximate expressions (3.27). This profile was also examined by Young & Wolfe (Reference Young and Wolfe2014), Bonfils et al. (Reference Bonfils, Mitra, Moon and Wettlaufer2022) and Abid & Kharif (Reference Abid and Kharif2023), sometimes with an extension into the water to account for a wind-driven surface drift layer:

(3.38) \begin{equation} \left. \begin{gathered} W(\kern0.7pt y) = W_{\infty}(1 - \exp{({-}y/y_s)}) , \quad \frac{y}{y_s} = \ln{\Big[\frac{W_{\infty} }{(W_{\infty} - W)}\Big]} ,\\ W_y = \frac{W_{\infty}}{y_s}\exp{({-}y/y_s)}, \quad \mathcal{S} = \frac{y_s }{(W_{\infty} - W)} ,\\ \mathcal{K} ={-} \frac{y_s}{(W_{\infty} - W )^2} , \quad \mathcal{K}_W ={-} \frac{2y_s}{(W_{\infty} - W )^3} . \end{gathered} \right\} \end{equation}

The expression (3.26) for $E$ becomes

(3.39) \begin{align} \left. \begin{gathered} E = \frac{y_s}{W_{\infty}}\left[\frac{W_{\infty}}{(W_0 - c_r)(W_{\infty} - W_{0})} + \frac{1}{c_r }\right] - \frac{1}{k(W_0 - c_r)^2} + \mathcal{K}_{Wc}W_0 + \mathcal{K}_c \ln{\left[\frac{(W_0 - c_r)}{c_r } \right]},\\ \mathcal{K}_c ={-} \frac{y_s}{(W_{\infty} - c_{r})^2 } , \quad \mathcal{K}_{Wc} ={-} \frac{2y_s}{(W_{\infty} - c_{r} )^3} . \end{gathered} \right\} \end{align}

Here $W(\kern0.7pt y) \to W_{\infty }$ as $y \to \infty$ and it turns out that the choice of $W_{\infty }$ and $y_s$ is quite sensitive. For instance if $W_{\infty }= 20 \ \textrm {m}\ \textrm {s}^{-1}$, $y_s = 10 \ \textrm {m}$ led to $W_{0B} = 13.15\ \textrm {m}\ \textrm {s}^{-1}$, $\beta _M= 36.30$, which seems too large by an order of magnitude. Decreasing $y_s$ or increasing $W_{\infty }$ led to an even larger value. As a guide we note that Morland & Saffman (Reference Morland and Saffman1993) put $W_{\infty } = 20 u_{*}$ and showed that then the growth rate depended on the single parameter $\varpi =2 y_s g/W_{\infty }^2$. With our $u_{*} = 0.36 \ \textrm {m}\ \textrm {s}^{-1}$, this leads to $W_{\infty } = 7.2\ \textrm {m}\ \textrm {s}^{-1}$, much smaller than we have used for the two other profiles, but in the range used by Young & Wolfe (Reference Young and Wolfe2014). Then $\varpi = 0.1$ yielded $y_s = 0.26\ \textrm {m}$, $W_{0B} = 3.35\ \textrm {m}\ \textrm {s}^{-1}$, $\beta _M = 3.45$ which is very consistent with the calculation of Morland & Saffman (Reference Morland and Saffman1993). However, decreasing $\varpi$ and so decreasing $y_s$, or increasing $W_{\infty }$ led to unsatisfactorily large values of $\beta _M$, as can be seen in Morland & Saffman (Reference Morland and Saffman1993). We conclude that although this wind shear profile looks promising because it is simple and the modal equation has an analytical solution, it is unsuitable in practice as it is very sensitive to the profile parameters.

5.1. Reduction to shallow-water equations

In this subsection we reduce the full system (2.10)–(2.13) to a modified form of the well-known shallow-water equations, as in Dutykh & Dias (Reference Dutykh and Dias2007) and Dutykh (Reference Dutykh2009). In the shallow-water asymptotic limit $H \to 0$, we make the approximation that $\partial /\partial x \ll \partial /\partial y$. Then at leading order $\phi \sim \varPhi (x,t)$ becomes independent of $y$ and Laplace's equation (2.10) yields

(5.1)\begin{equation} \phi \sim \varPhi(x,t) + (\kern0.7pt y +H)\phi_{y}(\kern0.7pt y={-}H) - \frac{(\kern0.7pt y+H)^2}{2}\varPhi_{xx} + \cdots . \end{equation}

Here $\phi _{y}(\kern0.7pt y=-H)$ is given by the boundary condition (2.30). As in Dutykh & Dias (Reference Dutykh and Dias2007) and Dutykh (Reference Dutykh2009), it is assumed that the frictional and forcing effects are weak and in particular this requires that the frictional boundary layers are thin and much smaller than the depth $H$, that is, $(\nu /\omega )^{1/2} \ll H$. As the depth decreases this requires a much weaker frictional effect. This is in addition to the requirement that the frictional boundary layers have small amplitude $(\nu /\omega )^{1/2} \ll |A|$.

The free-surface boundary condition (2.12) yields the equation

(5.2)\begin{equation} \varPhi_t +\frac{1}{2}\varPhi_x^2 + \eta ={-} \frac{P_a}{\rho_w}+ 2\nu \varPhi_{xx}. \end{equation}

The pressure forcing term is given by (2.36) using the shallow-water limits of $\alpha,\beta$. The second equation in this shallow-waterlimit is conservation of mass (2.28) given here by

(5.3)\begin{equation} \eta_t +((H + \eta )\varPhi_x )_x = 2\nu \eta_{xx} - \int_{0}^{t} \left(\frac{\tau_a}{\rho_w }\right)_x \,{\rm d}t \,+ \int_{0}^{t} \left(\frac{\tau_b }{\rho_w }\right)_x \,{\rm d}t . \end{equation}

Putting $U=\varPhi _x$ as the leading-order horizontal velocity, these take the form (see Dutykh & Dias Reference Dutykh and Dias2007)

(5.4)\begin{gather} U_t + UU_x + \eta_x = F_1 ={-} \frac{P_{ax}}{\rho_w} + 2\nu U_{xx} , \end{gather}

(5.5)\begin{gather} \eta_t +((H + \eta )U)_x = F_2 = 2\nu \eta_{xx} - \int_{0}^{t}\left(\frac{\tau_a}{\rho_w }\right)_x \,{\rm d}t + \int_{0}^{t} \left(\frac{\tau_b }{\rho_w }\right)_x \,{\rm d}t . \end{gather}

The pressure term $P_{a}$ is given by (2.36) using the shallow-water limits of $\alpha, \beta$, and the surface wave stress $\tau _a$ is given by the shallow-water limit of (2.44). Similarly the bottom stress $\tau _b$ is given by the shallow-water limits of either the laminar expression (2.29) or the turbulent expression (2.33):

(5.6)\begin{equation} \frac{\tau_{b[s]} }{\rho_w } = \left(\frac{\nu }{{\rm \pi} }\right)^{1/2}\,\int^{t}_{0} \frac{U_{t}(x,t-\sigma )}{\sigma^{1/2}}\,{\rm d}\sigma, \quad \frac{\tau_{b[t]} }{\rho_w } = C_D |U | U . \end{equation}

In the absence of the frictional effects these are the conventional shallow-water equations. This can be clearly seen if we replace $U$ by $\mathcal {U}$ where

(5.7)\begin{equation} \mathcal{U} = U + \int_{0}^{t}\left(\frac{\tau_a}{H\rho_w }\right) \,{\rm d}t - \int_{0}^{t} \left(\frac{\tau_b }{H\rho_w }\right) \,{\rm d}t - \frac{2\nu \eta_{x}}{H} , \end{equation}

which incorporates weak frictional and forcing effects. Then (5.4) and (5.5) become, keeping just leading-order terms,

(5.8)\begin{gather} \mathcal{U}_t +\mathcal{U}\mathcal{U}_x + \eta_x ={-} \frac{P_{ax}}{\rho_w} + 4\nu \mathcal{U}_{xx} + \frac{\tau_a}{H\rho_w } - \frac{\tau_b}{H\rho_w } , \end{gather}

(5.9)\begin{gather} \eta_t +\left((H + \eta )\mathcal{U}\right)_x = 0 . \end{gather}

In the variables $\mathcal {U}, \eta$ (5.4) and (5.5) are the conventional shallow-water equations with a frictional modification in the momentum equation (5.8).

5.2. Linearised shallow-water equations

Before proceeding with solutions of the system (5.4) and (5.5) it is useful to examine the linearised system given by

(5.10)\begin{gather} U_t + \eta_x ={-} \frac{P_{ax}}{\rho_w} + 2\nu U_{xx} , \end{gather}

(5.11)\begin{gather} \eta_t + H U_x = 2\nu \eta_{xx} + \int_{0}^{t} \left(\frac{\tau_b }{\rho_w }\right)_x \,{\rm d}t . \end{gather}

The turbulent term (2.44) for $\tau _a$ is omitted as it is nonlinear. Similarly for bottom stress we use only the laminar expression (5.6) in these linearised equations.

Eliminating $U$ from the system (5.10) and (5.11), substituting for $P_a$ from (2.36) and for $\tau _b$ from (5.6) yields

(5.12)\begin{equation} \eta_{tt}- H\eta_{xx} ={-}Hk^2(\alpha k\eta + \beta \eta_{x}) W_{r}^2 + 4 \nu \eta_{txx} - \left(\frac{\nu }{{\rm \pi} }\right)^{1/2}\,\int^{t}_{0} \frac{\eta_{xx}(x,t-\sigma )}{\sigma^{1/2}}\,{\rm d}\sigma. \end{equation}

In the absence of pressure forcing and friction ($\nu \to 0$) the system (5.10) and (5.11) is the well-known d'Alembert wave equation and has solutions describing two non-dispersive waves, moving with speeds $\pm H^{1/2}$ to the right and left, respectively. The dispersion relation (2.16) correspondingly reduces to

(5.13)\begin{equation} \omega^2 = Hk^2 . \end{equation}

This can easily be derived from the linearised equations (5.4) and (5.5) or by taking the limit $Q \to 0$ in the full dispersion relation (2.16); the latter requires that $\tanh {Q} \approx Q$, valid to an error of $1\,\%$ when $Q < 0.25$. Retaining the pressure forcing and friction terms leads to the expanded dispersion relation; see (2.37):

(5.14)\begin{equation} \omega^2 - Hk^2 = Hk^3 (\alpha + {\rm i}\beta ) W_{r}^2 - 4{\rm i}\nu \omega k^2 - k^2\left(\frac{\nu }{\omega {\rm \pi}}\right)^{1/2}\int^{\omega t}_{0} \frac{\exp{({\rm i}\varTheta)}}{\varTheta^{1/2}}\,{\rm d}\varTheta. \end{equation}

The integral term is the same $I(t)$ as in (2.31) and is evaluated as $({\rm \pi} /2)^{1/2}(1+\textrm {i})$ for the limit $\omega t \to \infty$.

We seek solutions of the full system (5.4) and (5.5) and of the linearised system (5.12) when the initial condition is that for a periodic wave (2.7). Solutions of the full system (5.4) and (5.5) are described in § 5.3 using Riemann invariants. For this linearised equation substitution of (2.7) into (5.12) for a slowly varying amplitude $A$ yields at leading order

(5.15)\begin{equation} 2{\rm i}\omega A_{t} + 2{\rm i}kH A_{x} = Hk^3 (\alpha + i\beta) W_{r}^2 A - 4{\rm i} \omega k ^2 \nu A - k^2\left(\frac{\nu }{2\omega }\right)^{1/2}(1 + {\rm i})A, \end{equation}

where from (2.16) $\omega = H^{1/2} k$ is here real-valued. Alternatively (5.15) follows from (5.14) by replacing $\omega, k$ with $\omega + \textrm {i}\partial /\partial t , k - \textrm {i}\partial /\partial x$.

The general solution of (5.15) is

(5.16) \begin{equation} \left. \begin{gathered} A = A_{0} (\chi)\exp{(\varDelta_{L} t )} , \quad \chi = x- H^{1/2}t ,\\ \varDelta_{L} = \frac{H^{1/2}k^2 }{2}( \beta -{\rm i}\alpha ) W_{r}^2 - 2k ^2 \nu - \frac{k^2}{2\omega}\left(\frac{\nu }{2\omega }\right)^{1/2}(1 - {\rm i}). \end{gathered} \right\} \end{equation}

The function $A_{0}(x)$ is determined by the initial conditions. Modulo the exponential growth the wave moves with the group velocity which in shallow water is $H^{1/2}$. As expected $\textrm {Re} \, \varDelta _{L} = \varDelta _{1} + \varDelta _{2[s]} +\varDelta _{3[s]}$ after evaluating each $\varDelta _{i} , i =1,2,3$, in the shallow-water limit; see (2.39), (2.40) and (2.32).

5.3. Riemann invariants

The system (5.4) and (5.5) is considered with an initial condition which generates at leading order a wave propagating in the positive $x$ direction. In the linearised system of § 5.2 this is a small-amplitude periodic wave; see (2.7):

(5.17)\begin{equation} \eta = A_{0}(x)\exp{({\rm i}kx)} + \text{c.c.}, \quad U = H^{{-}1/2} A_{0}(x)\exp{({\rm i}kx)} + \text{c.c.} \end{equation}

The initial condition for $U$ comes from the linearised equations (5.10) and (5.11) since to leading order for an unforced, friction-free linear wave moving to the right, $U =H^{-1/2}\eta$. The slowly varying envelope $A_{0}(x)$ can be specified arbitrarily, as in the linear solution (5.16), but we note that the choice $A_{0}(x) = M \operatorname {sech}(\gamma x),$ where $M$ is the amplitude of the carrier wave and $\gamma \ll k$, corresponds to the case 3 initial condition of Maleewong & Grimshaw (Reference Maleewong and Grimshaw2022a,Reference Maleewong and Grimshawb) in simulations of the fNLS equation and modified Euler equations.

The shallow-water system (5.4) and (5.5) has Riemann invariants $R_{\pm }$, where

(5.18)\begin{equation} \left. \begin{aligned} R_{{\pm}} & = U \pm 2\mathcal{D}, \quad \mathcal{D} = (H + \eta)^{1/2} - H^{1/2} , \\ \text{and so} & \quad U = \frac{1}{2}(R_{+} + R_{-}), \quad \mathcal{D} = \frac{1}{4}(R_{+} - R_{-}). \end{aligned} \right\} \end{equation}

Using these as new dependent variables, the system (5.4) and (5.5) becomes

(5.19)\begin{equation} \left. \begin{aligned} & \frac{\partial R_{{\pm}}}{\partial t} +V_{{\pm}} \frac{\partial R_{{\pm}}}{\partial x } =\mathcal{F}_{{\pm}} = F_1 \pm (H + \eta)^{{-}1/2}F_2 ,\\ & V_{{\pm}} ={\pm} H^{1/2} + U \pm \mathcal{D} ={\pm} H^{1/2} + \frac{3}{4}R_{{\pm}} + \frac{1}{4}R_{{\mp}} . \end{aligned} \right\} \end{equation}

In the frictional term $\mathcal {F}_{\pm }$, $F_{1,2}$ are the right-hand sides of (5.4) and (5.5), respectively. For a wave propagating in the positive $x$ direction in the unforced friction-free case $R_{-}=0$ is a constant and here we use that as a leading-order approximation with an error of order $\mathcal {F}_{-}$. Then the Riemann invariant $R_{+} \approx 2U$ and the speed $V_{+} =H^{1/2} +3U/2$ implying eventual wave breaking. More precisely, in order to take account of the frictional and forcing terms $\mathcal {F}_{\pm }$ we change the dependent variables from $x,t$ to $\chi, t$, where $\chi = x- H^{1/2}t$ is the dominant phase variable in the linear solution (5.16). In these variables the Riemann system (5.19) for a wave propagating in the positive $x$ direction uncouples and becomes

(5.20)\begin{equation} \frac{\partial R_{+}}{\partial t} + \frac{3}{4}R_{+}\frac{\partial R_{+}}{\partial \chi} =\mathcal{F}_{+}, \quad -2H^{1/2}\frac{\partial R_{-}}{\partial \chi} =\mathcal{F}_{-} . \end{equation}

The initial condition is (5.17) where we note that $x =\chi \ \textrm {at} \ t=0$. The friction terms are given by (5.19) and (2.36), (2.44) and (5.6) are all evaluated in the shallow-water limit. To leading order $\mathcal {F}_{\pm }$ are functions of $\chi$ with a weak $t$ dependence which enables an approximate evaluation of the terms involving $\tau _{a,b}$, where only the turbulent expressions are shown here. For the bottom stress the laminar term was used in the linearised shallow-water equations of § 5.2. Here we also omit the term in $\alpha$ in (2.36) so that the pressure condition is $P_{a} = \rho _a W_{r}^2 \beta \eta _x$. The outcome is

(5.21)\begin{equation} \mathcal{F}_{+} = \frac{\rho_{a}}{\rho_w}\beta k^2 W_{r}^2 \eta - 4\nu k^2 U + H^{{-}1} \frac{\rho_a}{\rho_w} c_d U_{a}^2 k^2 |\eta| \eta - H^{{-}1}C_D |U | U . \end{equation}

Since to leading order $R_{+} =2U =H^{-1/2} 2\eta$, the nonlinear terms $|\eta | \eta, |U|U$ in $\mathcal {F}_{+}$ are proportional to $|R_{+}| R_{+}$. Then (5.21) becomes

(5.22)\begin{align} \mathcal{F}_{+} = \varDelta_{N} R_{+}, \quad \varDelta_{N} = \frac{H^{1/2}k^2 }{2}\beta W_{r}^2 - 2k ^2 \nu + \frac{1}{2H^{3/2}} \frac{\rho_a}{\rho_w} c_d U_{a}^2 k^2 |A| - \frac{1}{2H^{3/2}} C_D |A | . \end{align}

The first two terms agree with those in the linearised expression $\textrm {Re}(\varDelta _{L})$ in (5.16) and the last two terms have been constructed from (2.45) and (2.35) using the approximation $R_{+} =2U =2H^{-1/2}\eta$. A plot of $\varDelta _{N}$ versus $H$ for the logarithmic wind shear profile and the benchmark parameters for a $5$ s wave and several amplitudes $A$ is shown in figure 11. For these parameters $\varDelta _{N}$ is positive indicating wave growth for $H > \tilde {H} \approx 1.5$ increasing as the wave amplitude increases. As the depth decreases further $\varDelta _{N} < 0$ indicating wave decay. This can be compared with figure 9 for $\varDelta$ and while there is qualitative similarity, some numerical differences appear because $\varDelta _N$ is based on the shallow-water approximation a priori, most obvious as $H \to 0$.

Figure 11. Plot of $\varDelta _{N}$ versus $H$ for a logarithmic wind shear profile (3.33) when $W_{0B} = 11.25\ \textrm {m}\ \textrm {s}^{-1}$ for several wave amplitudes $A$ (in units of $\textrm {m}$) when $0< H<1$.

The Riemann invariant equation (5.20) with friction term (5.21) can be solved exactly using characteristics, combined with the variable changes, $R_{+}=R\exp {(\varDelta _{N}t)}$ and $\textrm {d}T /\textrm {d}t = (3/4)\exp {(\varDelta _{N}t)}$. The outcome is

(5.23)\begin{equation} R = \text{constant when}\ \frac{{\rm d}\chi }{{\rm d}T} = R , \quad T = \frac{3}{4}\frac{\exp{(\varDelta_{N}t)} - 1}{\varDelta_{N}} . \end{equation}

The transformed time variable $T$ increases as $t$ increases, varying from $0$ to ($\infty, (3/4)|\varDelta _{N}|^{-1}$) when $\varDelta _{N} \gtrless 0$; as $\varDelta _{N} \to 0, T \to 3t/4$. The characteristics are labelled by the initial value $\chi _0$ so that $\chi - R_0 (\chi _0) T = \chi _0$ where the initial condition is that $R(\chi,T=0) = R_{0}(\chi )$. This produces wave steepening which is enhanced or opposed according as $\varDelta _{N} \gtrless 0$. Wave breaking may then occur and is defined here as intersecting characteristics when there is an infinite slope in $R$, which first occurs when $1 + R_{0\chi _0} T = 0$. With the initial condition (5.17) the explicit solution is

(5.24)\begin{equation} R = R_{0}(\chi- R T) , \quad R_{0}(\chi) = 4H^{{-}1/2}A_{0}(\chi )\cos{(k\chi )} , \end{equation}

where the smooth slowly varying wave envelope $A_{0}(\chi ) = M \operatorname {sech}(\gamma \chi )$, $\gamma \ll k$. The wave-breaking criterion $1 + R_{0\chi _0} T = 0$ is met when $\sin {(k\chi _0)} = 1$ and $4H^{-1/2}kM = 1/T$. Since $k= \omega H^{-1/2}$ in this shallow-water approximation the breaking time is $T= H/\omega M$ and just $H/M$ for a $5$ s wave as then $\omega =1$. Wave breaking always occurs for $\varDelta _{N}>0$ but would be prevented for $\varDelta _{N}<0$ as here, if the initial wave amplitude is small enough, $3M/H < |\varDelta _{N}|$. A plot of (5.24) is shown in figure 12 for two depths, $H=0.2, 0.5$, that is, $h = 1.25, 3.125 \ \textrm {m}$, for a $5$ s wave. The parameters in the initial condition $A_0 (\chi )$ are $M= 0.02, \gamma =0.5$ for depth $H=0.2$ (figure 12a) and $M= 0.1, \gamma =0.5$ for depth $H=0.5$ (figure 12b). The amplitude $M$ is chosen to ensure that the trough-to-crest height $2M$ is much less than $H$ in keeping with our small-amplitude assumption. The wavenumber $k= 2.236, 1.414$ corresponding to $H = 0.2 , 0.5$. These values of $k$ correspond to $Q = kH = 0.447, 0.707$, respectively, which puts them just outside the shallow-water regime, which strictly requires that $\tanh {Q} \approx Q$ valid to an error of $1\,\%$ when $Q < 0.25$. Nevertheless we continue to use these values of $H$ as even smaller values are physically unrealistic. The decay rate $\varDelta _{N} = -0.0075 , -0.0015$ for $H = 0.2, 0.5$, respectively. The predicted time for wave breaking is $T = 2.5, 1.25$, respectively, corresponding to $t = 3.38, 1.67$. This is much smaller than the decay rate as here $3M/H \gg |\varDelta _{N}|$. However, caution is needed here as this wave breaking is an artefact of the shallow-water approximation, as higher-order dispersive terms would inhibit the occurrence of an infinite slope. Indeed at the next order the shallow-water equations are replaced by a modified Boussinesq system (see Dutykh & Dias Reference Dutykh and Dias2007; Dutykh Reference Dutykh2009), which can in turn be further reduced to a modified Korteweg–de Vries equation with solitary wave solutions which grow or decay depending on the balance between forcing and friction. We are currently looking at this development. Also we note that in a fully nonlinear setting periodic water waves are unstable either when $|A|/H > 0.4$, or in deep water when the wave steepness $|A|>0.4$, expressed here in our non-dimensional units. Since $M/H = 0.1, 0.2$ these limits are avoided here.

Figure 12. Evolution of the Riemann invariant solution (5.24) for a $5$ s wave and the benchmark parameters for a logarithmic wind shear profile: (a) $H=0.2$, $M=0.02$; (b), $H=0.5$, $M=0.1$.