Evolution of wind-induced wave groups in water of finite depth

Abstract

The generation of water waves by wind is a fundamental and much-studied problem of scientific and operational concern. One mechanism that has obtained a wide degree of acceptance and use is the Miles air shear flow instability theory in which water waves grow due to energy transfer from the air at a critical level where the wave phase speed matches the wind speed. In this paper we revisit and extend this theory by examining how the predicted wave growth rate depends on the water depth and the wind shear parameters. At the same time we include frictional effects due to water bottom stress and at the water surface, and add an additional driving term due to wave stress in the air near the sea surface, using a turbulent parametrisation. The theory is developed using a frictional modification of the usual potential flow theory for water waves, the modified Euler equations, coupled to linearised air-flow equations. To assist our analysis we develop a long-wave approximation for the key Miles growth parameter, which explicitly shows how this depends on the water depth and the wind shear parameters. Since our focus is on the development of wave groups, we present a forced nonlinear Schrödinger equation in which the forcing term describes wave growth or decay. For very shallow water the modified Euler equations are reduced to their shallow-water counterpart, which are briefly discussed using Riemann invariants.

1. Introduction

The generation and evolution of water waves due to wind action at an air–water interface 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 ocean–atmosphere 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 is still not well understood. 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 contrasting opinions 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, water waves form into characteristic wave groups, and what their essential properties are, depending on the local atmospheric and oceanic conditions.

Over time, several mechanisms have been invoked to describe the generation and evolution of water waves by wind. One that is very well known is a classical shear flow instability mechanism introduced by John Miles in 1957 (Miles 1957), developed further by him and many others, and has been adapted for routine use in wave forecasting models (Miles 1959, 1993; Janssen 2004; Cavaleri et al. 2007; Grimshaw et al. 2018). The theory is based on linear sinusoidal waves with a real-valued wavenumber and a complex-valued frequency so that waves may have a temporal growth rate. There is a significant transfer of energy from the wind to the waves at the critical level in the air near the interface where the wave phase speed matches the wind speed. Independently, also in 1957, Owen Phillips developed a theory for water wave generation due to the flow of a turbulent wind over the sea surface (Phillips 1957). This mechanism is based on a resonance between a fluctuating pressure field in the air boundary layer and water waves due to a match between the water wave wavelength and the length scale of the pressure fluctuations. This leads to a linear growth in the water wave amplitude, contrasting with the exponential growth of the Miles theory. It is now widely believed that the Phillips mechanism applies in the initial stages of wave growth and that the Miles mechanism describes the later stage of wave evolution (Miles 1957; Phillips 1957, 1981). Another quite different mechanism is a steady-state theory, developed by Jeffreys (1925) for separated flow over large-amplitude waves, and later adapted for non-separated flow over low-amplitude waves (Belcher & Hunt 1998). Asymmetry in the free-surface profile is induced by an eddy viscosity closure scheme, which in the air flow allows for an energy flux to the waves.

No single or combined theory has been found completely satisfactory, and most fail to take account of wave transience and the tendency of waves to develop into wave groups; see Zakharov et al. (2015), Zakharov, Resio & Pushkarev (2017) and Zakharov (2018) amongst many similar criticisms. That is the issue we have been looking at (see Grimshaw 2018, 2019a,b; Maleewong & Grimshaw 2022a,b). Like Miles (1957), our analysis is based on linear shear flow instability theory, but incorporates from the outset that the waves will have a wave group structure with both temporal and spatial dependence. The key feature is that the wave group moves with a real-valued group velocity even for unstable waves when the wave frequency and the wavenumber are both complex-valued. In the absence of wind forcing it is well known that the nonlinear Schrödinger (NLS) equation 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 e.g.  Grimshaw 2007; Osborne 2010). Recently it has been proposed that the effect of wind forcing can be captured by the addition of a linear growth term leading to a forced NLS (fNLS) equation (see Leblanc 2007; Touboul et al. 2008; Kharif et al. 2010; Onorato & Proment 2012; Montalvo et al. 2013b; Brunetti et al. 2014; Slunyaev, Sergeeva & Pelinovsky 2015; Grimshaw 2018, 2019a,b; Maleewong & Grimshaw 2022a,b).

In contrast, most of the literature on wind-generated water waves has focused on the development and analysis of the statistical spectrum using the well-known Hasselman equation which describes the evolution of the water wave spectrum represented by the wave action under the influence of nonlinearity due to a field of resonant quartet interactions, a wind forcing source term and wave dissipation, mainly due to wave breaking (see e.g. Janssen (2004); Grimshaw et al. (2018)). While there are many associated analytical and numerical studies of the fully nonlinear Euler equations for water waves, there are comparatively few studies of a fully nonlinear two-fluid air–water system. Much of these have focused on modelling turbulence in the air flow; see for instance the articles by Sajjadi, Drullion & Hunt (2018), Sullivan et al. (2018), Wang, Yan & Ma (2018) and Hao et al. (2018) in the 2017 IUTAM proceedings (Grimshaw et al. 2018). One reason for this is because an inviscid fully nonlinear two-fluid system such as air over water is subject to short-scale Kelvin–Helmholtz instability. While this is a real physical effect, it is not the explanation for the growth of wind waves, which is due to wind shear and critical-level instability in the near-surface air flow. This difficulty can be avoided by replacing the two-fluid system with a fully nonlinear inviscid Euler system for the water, driven by a pressure term at the free surface which directly links the free-surface pressure with the free-surface slope. This modified Euler system is based on the pioneering work of Miles (1957) and was later developed by Kharif et al. (2010) amongst others.

Our interest is in the growth of wave groups, which in the presence of wind forcing can be described by a fNLS equation, which is a weakly nonlinear asymptotic reduction of the modified Euler system (Maleewong & Grimshaw 2022a,b). In this fNLS model, the key new feature imposed on the well-known NLS equation is a linear growth term which causes exponential growth of the wave amplitude at a rate described by a key parameter, $\Delta s^{-1}$, where asymptotically $\varDelta$ is much smaller than a typical wind wave period. Parameter $\varDelta$ depends on the wave parameters, the fluid depth and the wind shear profile parameters and crucially on a non-dimensional parameter, $\beta$, introduced by Miles (1957), who showed that in a linear air-flow analysis the surface air pressure contains a term $\rho _a \beta \eta _x W_{r}^2$, where $\rho _a$ is the air density, $\eta _x$ is the wave slope and $W_r$ is a reference scaling velocity. It is this term which generates the critical-level instability. We note that this theory is unidirectional in the horizontal. A two-dimensional counterpart was introduced by Benney & Roskes (1969) and a wind-forced version analysed by Grimshaw (2019b) and Maleewong & Grimshaw (2023).

In this paper we examine in § 3 in detail how $\beta$ and, hence, $\varDelta$ depend on the wave parameters, the water depth and the parameters defining the wind shear profile. For this purpose we use a long-wave approximation described in § 3.1. Miles (1957) and many others used a logarithmic wind shear profile (see § 3.2), but here we extend that to examine two other similar wind shear profiles, in § 3.3 for an algebraic wind shear profile and in § 3.4 for an exponential wind shear profile. As well as the critical-level instability we examine frictional and wave stress effects at the air–water interface, and at the bottom of the water column. This leads to the sum $\varDelta = \varDelta _1 + \varDelta _2 +\varDelta _3$, where $\varDelta _1$ is a linear combination of wave growth due to critical-level instability and wave decay due to a laminar frictional boundary layer at the water surface, $\varDelta _2$ is a wave growth term due to a parametrisation of turbulent wave stress in the air flow near the water surface and $\varDelta _3$ is a wave decay term due either to a laminar frictional boundary layer at the water bottom or to a parametrisation of turbulent stress at the water bottom. To implement our analysis we develop a long-wave approximation for the key Miles growth parameter $\beta$, which explicitly shows how this depends on the water depth and the wind shear parameters. Our analysis is based on the weak frictional modification of potential flow introduced by Dutykh & Dias (2007) and Dias, Dyachenko & Zakharov (2008) (see also Longuet-Higgins 1992; Dutykh 2009; Kharif et al. 2010) and described here in §§ 2 and 2.1. In § 2.2 we describe wave stresses at the air–water interface and at the bottom, and in § 2.3 we describe the critical-level instability theory of Miles (1957) for wind forcing. In § 3 we present the linearised air-flow equations, and develop the aforementioned long-wave approximation for the parameter $\beta$. In § 4 we briefly describe the fNLS equation, and present there in analytical and graphical form a description of the dependence of $\varDelta$ on water depth, the wave parameters and some key parameters of the wind shear profile, using the logarithmic profile as the main example. In § 5 we present the fully nonlinear modified Euler system in a shallow-water limit, following Dutykh & Dias (2007). In § 5.2 we solve the linearised shallow-water equations for an initial condition of a linear periodic wave with a slowly varying envelope amplitude, and then in § 5.3 we extend this to the fully nonlinear shallow-water system using Riemann invariants. As this inevitably leads to wave steepening, higher-order dispersion is reintroduced in a forced Korteweg–de Vries model, which is the subject of ongoing research. We conclude with a summary and discussion in § 6.

2. Formulation

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 (1992), Dutykh & Dias (2007) and Dias et al. (2008) to account for thin laminar frictional boundary layers at the free surface and at the water bottom, and modified by Kharif et al. (2010) 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 (2007) 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. (1996) 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 1968)

(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 (1957) 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 (1957) 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 (1957) 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 (1957) 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 (1957) 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 1959). We re-examine this again in § 3.

As in Miles (1957), Janssen (2004), Kharif et al. (2010) 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. 2010), 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 2018, 2019b). 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 (2022a,b). Note that the critical layer instability mechanism of Miles (1957) 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 (2007) and Kharif et al. (2010). 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 (1957), 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. 1996). We note here that Branger et al. (2022) 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. Linearised air flow equations

3.1. Formulation and long-wave approximation

As in Miles (1957) 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 (2018, 2019a). 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 (1957) 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 2018, 2019b). 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 1957; Janssen 2004). 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 (1957) 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 (1957) 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 (1957) 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 (1959) 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. 2013a,b; Latifi et al. 2017; Branger et al. 2022) 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 (2004) and Grimshaw (2018, 2019a). 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 (1957). 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 (1957, 1959), 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 (1957) 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. (2022) for instance). The form (3.33), also used by Janssen (2004), 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 (1959) 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 (1959) 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 1975; Hsu, Meindl & Gilhousen 1994), 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 (1993) 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 (2014), Bonfils et al. (2022) and Abid & Kharif (2023), 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 (1993) 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 (2014). 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 (1993). 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 (1993). 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.

4. Forced NLS equation

In the absence of wind forcing, the full Euler equations can be reduced to a NLS equation for the description of a weakly nonlinear wave packet (see Benney & Newell 1967; Zakharov 1968; Hasimoto & Ono 1972; Grimshaw 2007). In the presence of wind forcing the outcome is a fNLS equation (see e.g. Leblanc 2007; Touboul et al. 2008; Kharif et al. 2010; Montalvo et al. 2013b; Onorato & Proment 2012; Brunetti et al. 2014; Slunyaev et al. 2015; Grimshaw 2018, 2019a,b; Maleewong & Grimshaw 2022b, a). Here we give a brief summary to indicate the forcing term, which contains the growth rate $\varDelta = \varDelta _1 +\varDelta _2 + \varDelta _3$ as described in § 2. We then examine each $\varDelta _i$ as a function of depth $H$, the wave and shear profile parameters. We impose a weakly nonlinear asymptotic expansion, expressed here in the scaled non-dimensional variables:

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

with a corresponding expression for $\phi (x, y, t)$. Here $\textrm {c.c.}$ is the complex conjugate and $\cdots$ are higher-order terms in the asymptotic expansion. Term $A(x,t)$ is a slowly varying amplitude and the expansion is jointly with respect to the amplitude and this slow variation. Note that the wave height is twice the crest (trough) height above the undisturbed level.

At leading order one gets the linear dispersion relation (2.16). At the next order one finds that the amplitude $A$ propagates with the group velocity $c_g$ (2.18). Second-order terms yield the second harmonic and a mean flow term. At the third order a compatibility condition yields the fNLS equation:

(4.2)\begin{equation} {\rm i}(A_t + c_g A_x) + \lambda A_{xx} + \mu |A|^2 A= {\rm i}\,\Delta A . \end{equation}

The coefficients $\mu$ and $\lambda$ are given by (see e.g. Hasimoto & Ono 1972)

(4.3)\begin{align} \mu &={-}\frac{\omega k^2}{4 \tanh^4{(Q)} } (9\tanh^4{(Q)} -10\tanh^2{(Q)} +9) \nonumber\\ &\quad +\frac{\omega^3}{2 \tanh^3{(Q)}(H- c_{g}^{2})}(2\tanh{(Q)} (3-\tanh^2{(Q)} ) +3Q(1-\tanh^2{(Q)} )^2 ), \end{align}

(4.4)\begin{equation} \lambda = \frac{\omega_{kk}}{2} , \quad Q =kH . \end{equation}

Coefficient $\lambda < 0$ for all $Q$ while $\mu < 0$ $( > 0)$ according as $Q> Q_c$ $( Q< Q_c)$, where $Q_c = 1.363$. In deep water ($Q \to \infty$), $\mu \to -2\omega k^2$ and $\lambda \to -\omega /8k^2$. In the shallow-water limit $H \to 0$ discussed in § 5, $\lambda \to 0$ and $\mu \to \infty$ thus invalidating the fNLS model. For a fixed but small frequency $\omega$, $k \sim \omega H^{-1/2}$, $Q=kH\sim \omega H^{1/2}$ as $H \to 0$, and then $\lambda \sim -\omega H^2 /2$, $\mu \sim 9/4\omega H^3$.

In the absence of wind forcing, modulation instability occurs when $\mu \lambda > 0$, that is, when $\mu < 0$, $Q > Q_c$, leading to the formation of solitons and breathers as models of wave packet envelopes. In Maleewong & Grimshaw (2022a,b) we showed that inclusion of wind forcing represented by $\varDelta$ in (4.2) modulation instability again occurs but the wave envelope grows exponentially at a rate $2\varDelta$, twice the linear rate. This follows from the energy law (2.21) which continues to hold for (4.2) with $\mathcal {E}$ represented by $\int ^{\infty }_{-\infty } |A|^2 \, {\textrm {d} x}$.

The total growth rate $\varDelta =\varDelta _1 +\varDelta _2 +\varDelta _3$. Term $\varDelta _1$ (2.39) arises from the modified pressure term (2.36) and is due to critical-level instability in the air flow. Term $\varDelta _2 = \varDelta _{2[s]} + \varDelta _{2[t]}$ ((2.40) and (2.45)) arises from friction in the near-surface water boundary layer and turbulent wind stress in the near-surface air flow, respectively. Term $\varDelta _3$ arises from the bottom friction term, either laminar (2.32) or turbulent (2.35). For convenience each of these are repeated here, but using dimensional variables in units of $\textrm {s}^{-1}$, instead of non-dimensional variables. At the same time we also express the growth rate $\omega _i$ for $\varDelta _1$ (2.38) in dimensional variables:

(4.5) \begin{equation} \left. \begin{gathered} \varDelta_1 =\frac{\rho_a}{2 \rho_w} \beta \omega_r \tanh{Q} \frac{W_{r}^{2}}{c_{r}^2} , \quad \frac{\omega_i }{\omega_r } = \frac{\rho_a}{\rho_w}\frac{1}{2g}\beta k W_{r}^2 ,\\ \varDelta_{2[s]} ={-}2 k^2 \kappa , \quad \varDelta_{2[t]} = \frac{\rho_a}{\rho_w} c_d U_{a}^2 \frac{16k^2 |A|}{3c_{r}} ,\\ \varDelta_{3[s]} ={-} \left(\frac{\kappa}{2 \omega_r}\right)^{1/2} \frac{gk^2}{2\omega_r}\frac{1}{\cosh^2{Q}} , \quad \varDelta_{3[t]} ={-} \frac{16 C_D \omega_r k|A|}{3 \cosh{Q} \sinh^2{Q}} . \end{gathered} \right\} \end{equation}

Each of these is plotted as a function of depth $H$ in figures 5–8 for the logarithmic wind shear profile (3.33) of § 3.2. In these plots we first choose the scaling parameter $\varOmega = 2{\rm \pi} /5$ corresponding to a $5$ s wave in deep water. Then $K =\varOmega ^2 /g =0.161 \ \textrm {m}^{-1}$, $\varOmega K^{-1}= g \varOmega ^{-1} = 7.81 \ \textrm {m}\ \textrm {s}^{-1}$. Each $\varDelta$ depends on the wave parameters $\omega, k$ and the depth parameter $H$, where the linear dispersion relation (2.16) reduces this to just two of these three wave parameters. In practice, we fix $\omega = 1.257 \ \textrm {s}^{-1}\,$ corresponding to a $5$ s carrier wave, and then $k$ is found from (2.16) as a function of $H$, where $k$ increases as $H$ decreases; see (2.17a,b). In addition, each $\varDelta$ depends on other physical parameters associated with the water and air properties, and on several other parameters associated with the choice of wind shear profile, importantly $\alpha, \beta, W_{0B}$, where for the logarithmic wind shear profile (3.33) $W_{0B} = 11.25\ \textrm {m}\ \textrm {s}^{-1}$. In particular we set the non-dimensional drag coefficients $C_D, c_d \sim 10^{-3}$ and the water kinematic viscosity $\kappa \sim 10^{-6} \ \textrm {m}^2 \ \textrm {s}^{-1}$. The surface wind $u_{a}$ is specified by (2.43).

Figure 5. Plot of $\varDelta _1$ (in units of $\textrm {s}^{-1}$) versus $H$, for a logarithmic wind shear profile (3.33) for several $W_{0}$ around $W_{0B}$ (in units of $\textrm {m}\ \textrm {s}^{-1}$).

Figure 6. Plots of (a) $\varDelta _{2[s]}$ and (b) $\varDelta _{2[t]}$ (in units of $\textrm {s}^{-1}$) versus $H$ for a logarithmic wind shear profile (3.33) when $W_{0B} = 11.25\ \textrm {m}\ \textrm {s}^{-1}$ and for several wave amplitudes $A$ (in units of $\textrm {m}$).

Figure 7. Plot of $\varDelta _{2} =\varDelta _{2[s]} +\varDelta _{2[t]}$ (in units of $\textrm {s}^{-1}$) 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}$).

Figure 8. Plots of (a) $\varDelta _{3[s]}$ and (b) $\varDelta _{3[t]}$ (in units of $\textrm {s}^{-1}$) 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}$).

Figure 5 for $\varDelta _{1} > 0$ ((2.39) and (4.5)) shows a significant increase in magnitude as $H$ decreases, most marked for very shallow water where $H = Kh \ll 2$, $h \ll 12.5 \ \textrm {m}$. For the most part this can be traced to the corresponding change in $\beta$ shown in figure 3, and described in § 3.2. For $W_0 =W_{0B} = 11.25\ \textrm {m}\ \textrm {s}^{-1}$, $\varDelta _1$ approaches a constant value in the deep-water limit $H \to \infty$, and then as $H$ decreases there is a very slight increase with a barely discernible maximum before decreasing to zero as the depth $H \to 0$. For $W_0 > W_{0B}$, $\varDelta _1$ is smaller and decreases as the depth decreases. When $W_0 < W_{0B}$, $\varDelta _1$ is again smaller in the deep-water limit, but as the depth decreases there is now a marked increase in magnitude to a maximum, more marked the further $W_0$ is decreased. For these parameters $\varDelta _1$ is quite small generating an $e$-folding time scale of order $10^5 \ \textrm {s}$, about one day. This is because we have evaluated it for a $5$ s wave, and the factor $\exp {(c_r/W_{r})}$ in $\beta$, see (3.34), varies exponentially with $c_r$. Hence we tested this by re-evaluating for a $2.5$ s wave, when $c_r$ is reduced from $7.8 \ \textrm {m}\ \textrm {s}^{-1}$ in deep water to $c_r = 3.9 \ \textrm {m}\ \textrm {s}^{-1}$. The outcome is that $\beta$ is increased tenfold and $\varDelta _1$ is reduced to order $10^{-3} \ \textrm {s}^{-1}$, with a dramatic decrease in the $e$-folding time scale to order $10^3 \ \textrm {s}$, less than one hour. This illustrates the extreme sensitivity to the choice of parameters.

Figure 6 for $\varDelta _{2[s]} <0, \varDelta _{2[t]}>0$ ((2.40), (2.45) and (4.5)) and figure 7 for the sum $\varDelta _{2}= \varDelta _{2[s]} + \varDelta _{2[t]}$ show a substantial increase in magnitude as the depth decreases, due to the decrease in $c_r$ as $H \to 0$. This is most marked in very shallow water, $H \ll 0.2$ for $\varDelta _{2[s]}$ and $H \ll 0.02$ for $\varDelta _{2[t]}$; these shallow depths are too small to be of much physical interest. In deep water as $H \to \infty$, $\varDelta _{2[s]} \to -4 \times 10^{-8}\ \textrm {s}^{-1}$ and $\varDelta _{2[t]} \to [-0.9,-3,-5,-7]\times 10^{-11}\ \textrm {s}^{-1}$ for $|A| = [ 0.01, 0.03, 0.05, 0.07]\ \textrm {m}$; $\varDelta _{2[s]}$ dominates over $\varDelta _{2[t]}$ for this deep-water limit. Here $\varDelta _{2[s]}$ has an inverse time scale comparable with $\varDelta _{1}$ and they are often compared together in estimating a time scale for critical-level instability (see Miles 1957; Kharif et al. 2010). For larger depths, $H >1$, $\varDelta _1$ in figure 5 is of order $10^{-5}\ \textrm {s}^{-1}$ but as $H$ decreases to the very shallow limit $H \ll 0.02$, $\varDelta _{2}$ in figure 7 of order $10^{-6}\ \textrm {s}^{-1}$ is dominant. Since the expression (2.45) for $\varDelta _{2[t]}$ has a dependence on $A$, the consequent growth in $A$ is algebraic and not exponential; the energy $\mathcal {E} = \mathcal {E}(t=0) \{1 - C_0\, t \}^{-2}$ and becomes singular on a long time scale as the constant $C_0 \propto |A|(t=0)$. As $H$ approaches zero $\varDelta _{2[s]}$ and $\varDelta _{2[t]}$ are of the same order of magnitude, $10^{-6}\ \textrm {s}^{-1}$. For larger depths, $H >1$, $\varDelta _1$ in figure 5 is of order $10^{-5}\ \textrm {s}^{-1}$ but as $H$ decreases to the very shallow limit $H \ll 0.02$, $\varDelta _{2[s]}$ in figure 6 is of order $10^{-6}\ \textrm {s}^{-1}$ and becomes dominant. Figure 7 for the sum $\varDelta _2$ shows that this is negative in deep water but can become positive for a large enough wave amplitude $A$ as $H$ decreases to a very small value.

For bottom friction, plots of laminar friction $\varDelta _{3[s]} <0$ and turbulent friction $\varDelta _{3[t]} < 0$ ((2.32), (2.35) and (4.5)) are shown in figures 8(a) and 8(b), respectively. Both are negative; in the deep-water limit, $\varDelta _{3[s]}$ has the larger magnitude for all $A$ while $\varDelta _{3[t]}$ of order $10^{-1}\ \textrm {s}^{-1}$ has much the larger magnitude in the shallow-water limit. In the deep-water limit, $\varDelta _{3[s]} \to 0$, $\varDelta _{3[t]} \to 0$ for all $A$ as expected due to the exponential factors in their respective formulae, (2.32) and (2.35). Both $\varDelta _{3[s]}$ and $\varDelta _{3[t]}$ are effectively zero of order $10^{-13}$ and $10^{-17}\ \textrm {s}^{-1}$ for $H > 10$.

A plot of the combined $\varDelta$ (2.22), with the bottom stress given by the turbulent expression, $\varDelta _3 = \varDelta _{3[t]}$, is shown in figure 9 for the logarithmic wind shear profile (3.33) and the benchmark parameters $W_0=11.25\ \textrm {m}\ \textrm {s}^{-1}$, $y_s=0.0002\ \textrm {m}$ for a $5$ s wave. There is a critical depth $H^{*} \approx 3.5$ such that when $H > H^{*}$, $\varDelta$ approaches a constant, the same value as in the full deep-water limit. The corresponding depth $h^{*} \approx 22\ \textrm {m}$. In this case $\varDelta$ is not very sensitive to the amplitudes in the range $A = 0.01\unicode{x2013} 0.1$. The value of $\varDelta$ is positive for $H > \tilde {H}$ indicating wave growth, where for these parameters $\tilde {H} \approx 1.75$ ($\tilde {h} \approx 11.0 \ \textrm {m}$), increasing as the wave amplitude increases.

Figure 9. Plot of $\varDelta$ (in units of $\textrm {s}^{-1}$) versus $H$ for a logarithmic wind shear profile (3.33) for the benchmark values $W_0=11.25\ \textrm {m}\ \textrm {s}^{-1}$, $y_s = 0.0002 \ \textrm {m}$.

The plots (mostly not shown here) were repeated for the logarithmic wind shear profile (3.33) for a $5$ s wave, varying $W_0, y_s \,$. Overall the results are similar. For fixed $y_s = 0.0002\ \textrm {m}$ and $W_0 = 13 > W_{0B} =11.25\ \textrm {m}\ \textrm {s}^{-1}$, $\varDelta$ is reduced for all $H$. But when $W_0=10 < W_{0B}=11.25\ \textrm {m}\ \textrm {s}^{-1}$ and again for fixed $y_s=0.0002 \ \textrm {m}$, $\varDelta$ decreases slightly for $H < H^{*}$ but then exhibits anomalous behaviour as $H \to 0$, with substantial growth to a positive value and then eventual decay (see figure 10). On the other hand, as $y_s$ is varied for a fixed $W_0 = W_{0B} =11.25\ \textrm {m}\ \textrm {s}^{-1}$, we found that for a larger $y_s = 0.002\ \textrm {m}$ $\varDelta$ is reduced for all $H$, but for a smaller $y_s = 0.00002\ \textrm {m}$ $\varDelta$ increases and the anomalous behaviour at very small $H$ reappears. We also carried out similar but less detailed examination for $\varDelta$ as regards $H$ for the algebraic wind shear profile (3.35) and the exponential wind shear profile (3.38). The results, not shown here, were broadly similar, albeit with different values for $W_0, y_s$.

Figure 10. Plot of $\varDelta$ (in units of $\textrm {s}^{-1}$) versus $H$ for a logarithmic wind shear profile (3.33) for $W_0=10\ \textrm {m}\ \textrm {s}^{-1} < W_{0B}$, $y_s = 0.0002 \ \textrm {m}$.

5. Shallow water

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 (2007) and Dutykh (2009). 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 (2007) and Dutykh (2009), 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 2007)

(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 (2022a,b) 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 2007; Dutykh 2009), 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$.

6. Summary and discussion

In this paper we have used the weak frictional and forcing modification of potential flow developed by Dutykh & Dias (2007), Dutykh (2009) and Kharif et al. (2010), presented here in § 2.1, to describe wind wave growth due to critical-layer instability, and the modifications due to wave stresses, laminar and turbulent, at the air–water interface and at the water bottom (see §§ 2.2 and 2.3). Our focus is on the evolution of wave groups described by the fNLS equation, (4.2) in § 4, for the slowly varying amplitude $A$ of a wave packet with carrier frequency and wavenumber $\omega, k$; see (4.1). Our aim is to determine how the growth rate $\varDelta$ in that equation depends on the wave parameters, the water depth, the frictional coefficients and the parameters defining the wind shear profile. For this purpose in § 3 we re-examine the linearised air-flow equations used in the pioneering work of Miles (1957). In § 3.1 we introduce a long-wave approximation (3.14) to obtain an explicit formula for the key parameter $\beta$ introduced by Miles (1957) which determines the growth rate due to critical-level instability (2.38). As well as re-examining the well-known logarithmic wind shear profile in § 3.2, we also examine two other similar wind shear profiles, algebraic and exponential, in §§ 3.3 and 3.4. The outcome for $\varDelta$ is described graphically in § 4. The main effects due to water depth occur at very small depths $H =Kh$, the non-dimensional depth; here $K$ is the scaling wavenumber, see (2.9), used to cast the equations in non-dimensional variables and is also the wavenumber of the underlying carrier wave in deep water, while $h$ is the dimensional depth. In § 5 we briefly describe the reduction of the full system in the shallow-water limit, leading to a frictional and forcing modification of the well-known shallow-water equations.

The total growth rate $\varDelta = \varDelta _1 + \varDelta _2 +\varDelta _3$ is composed of three parts, (2.22): (1) $\varDelta _1$ is positive and is due to the critical-level instability term; (2) $\varDelta _2 = \varDelta _{2[s]} + \varDelta _{2[t]}$ where the first term is negative and is due to laminar friction at the water surface and the second term is positive and is due to turbulence in the air near the air–water interface; and (3) either $\varDelta _3 = \varDelta _{3[s]}$ or $\varDelta _3 = \varDelta _{3[t]}$, both being negative and are respectively due to either laminar or turbulent bottom friction. The dimensional expressions for each of these are summarised in (4.5). The combination $\varDelta$ can be positive expressing wave growth, or negative expressing wave decay, depending inter alia on the water depth. In the deep-water limit the bottom stress is zero, effectively, so for $H > 10, h > 62.5 \ \textrm {m}$ and then $\varDelta = \varDelta _1 + \varDelta _2$ depends only on the carrier wave and wind shear profile parameters. Parameter $\varDelta _1 >0$ depends mainly on the value of the wind shear curvature $\mathcal {K}_c$ ((3.17) and (3.20)) at the critical level, while $\varDelta _{2} = \varDelta _{2[s]} + \varDelta _{2[t]}$ can be positive or negative depending on the relative values of $\varDelta _{2[s]} <0, \varDelta _{2[t]}>0$, that is, on the water laminar boundary friction at the water surface as regards the turbulence in the air near the air–water interface. In practice, it is usually assumed that in deep water $\varDelta > 0$. Each $\varDelta _{i}, i =1,2,3$, increases in magnitude as $H$ decreases, as seen in figures 5–8. However, although the largest magnitudes occur for such small $H$ that any physical significance can be discounted, nevertheless the tendency to increase in magnitude as the depth decreases is potentially important and we note that it is due in part to a decrease in the phase velocity $c_r = \omega / k$ as the depth decreases. Eventually, as the depth decreases to a very small value where bottom stress is the dominant feature, $\varDelta < 0$.

The linearised air-flow equations used by Miles (1957) and many others are re-examined in § 3 for three wind shear profiles $W(\kern0.7pt y)$. First, in § 3.1 we present the main result of this paper which is a long-wave approximation (3.14) designed to provide an explicit analytic expression for the key parameter $\beta$ introduced by Miles (1957) to determine the growth rate due to critical-level instability, (2.38). In § 3.2 we examine the well-known logarithmic wind shear profile used by Miles (1957) and many others, and then more briefly we examine an algebraic and an exponential profile in §§ 3.3 and 3.4. For each, the key parameters are a reference velocity $W_r$, a length scale $y_s$ and an upper limiting velocity $W_0$, where $W(\kern0.7pt y)$ increases monotonically from $0$ to $W_0$ as $y$ increases. For the logarithmic wind shear profile we set $W_r = 0.9\ \textrm {m}\ \textrm {s}^{-1}$ and $y_s = y_{*} = 0.0002\ \textrm {m}$, a roughness length scale, using empirical expressions as in the seminal paper by Miles (1957), although we also varied $y_s$ around this value. For the algebraic and exponential wind shear profiles we kept $W_r = 0.9\ \textrm {m}\ \textrm {s}^{-1}$ but of necessity choose different values for $y_s$ and $W_0$. We found that while the choice of wind shear profile and associated parameters led to similar qualitative outcomes, there were considerable quantitative differences. In particular the precise value of $W_0$ was important, even though it was constrained so that there is a critical level at all depths, and could not be too large. For the logarithmic wind shear profile of § 3.2 we chose $W_0 =W_{0B} =11.25 \ \textrm {m}\ \textrm {s}^{-1}$ since this maximised the value of $\beta = 1.745$ and hence also of $\varDelta _1 >0$ as described in § 3.2. In § 3.3 we followed a similar strategy for an algebraic profile; when $W(\kern0.7pt y)$ increases quadratically with $y$ which leads to $W_0 =W_{0B} = 39.17 \ \textrm {m}\ \textrm {s}^{-1}$ and $\beta = 0.99$. Then in § 3.4 we examined an exponential profile which yielded $W_0 =W_{0B} = 3.35 \ \textrm {m}\ \textrm {s}^{-1}$ and $\beta = 3.45$ albeit for a carefully chosen set of parameters.

This sensitivity to the choice of wind shear profile and to $W_0$ in particular is a major concern for the implementation of analytical theories such as this one into operational wind wave models; see Janssen (2004) and Grimshaw et al. (2018) for instance, where the Miles critical-level theory is adapted to produce a growth term in the Hasselman equation for the evolution of the water wave spectrum represented by the wave action density.

In order to examine the shallow-water limit $H \to 0$ in slightly more detail we considered the shallow-water reduction of the full system in § 5.1, keeping the frictional and forcing modifications as in Dutykh & Dias (2007) and Dutykh (2009). First in § 5.2 we solved the linearised shallow-water equations for a wave packet which propagates with the linear group velocity, $H^{1/2}$ in the shallow-water limit, with an exponential growth at rate $\varDelta _L$, (5.16), the linearised shallow-water approximation of $\varDelta$. Then in § 5.3 we examined the full nonlinear shallow-water system using Riemann invariants to describe a wave packet moving in the positive $x$ direction, again with an exponential growth rate, now $\varDelta _N$, (5.22), the shallow water reduction of $\varDelta$. The main new feature to emerge is that wave breaking indicated by an infinite slope will occur unless $\varDelta _N$ is sufficiently negative. The breaking time is quite short, much less than the time scale $|\varDelta _{N}|^{-1}$, and is an artefact of the shallow-water approximation. Retention of the next-order dispersive terms leads to a modified Boussinesq system (Dutykh & Dias 2007; Dutykh 2009). In ongoing work we take that step slightly further and are examining a modified Korteweg–de Vries equation to describe wind waves in shallow water. We note that the theory based on the fNLS equation is unidirectional in the horizontal. A two-dimensional counterpart was introduced by Benney & Roskes (1969) and a wind-forced version analysed by Grimshaw (2019b).

Our two main aims in this work were to examine how the key parameters in the laminar air-flow model of Miles (1957) vary with the wave parameters and the fluid depth, and how sensitive is the predicted wave growth rate $\varDelta$ to these parameters and to the choice of wind shear profile. Our focus is on the emergence of wave groups modelled by the fNLS equation. In assessing the growth rate, we used the weak frictional modification of potential flow of Dutykh & Dias (2007) and Dutykh (2009) to include laminar frictional effects in the water surface boundary layer, and at the water bottom, where we also examined the alternative of a turbulent parametrisation of bottom stress. At the same time we inserted a turbulent parametrisation of the near-surface wind wave stress as an additional driving term to the Miles critical-level instability theory. Overall the growth rate increases in magnitude as the depth decreases but a significant and anomalous increase only occurs at very small depths, too small in the context of this model to be of practical importance. However, there was considerable sensitivity to the parameters of the laminar air-flow model, where we mainly examined the well-known logarithmic wind shear profile, but found a similar outcome for two other wind shear profiles, algebraic and exponential. We conclude that using a laminar air-flow model to predict a wave growth rate for the Hasselman equation for the evolution of the water wave spectrum is potentially difficult. As is well known, although the growth rate predicted by the Miles critical-level instability theory is widely used in practice, it needs adjustment to an increased value to take account of observations, and to put operational wind wave forecasting on a secure level. A related and deeper issue is whether a laminar air-flow model can be used at all, when usually the observed winds are turbulent. For many years the answer to this was to choose the laminar wind shear profile to be an averaged flow, hence the standard choice of a logarithmic wind shear profile. Even with the advent of more and more powerful computers, we suggest that a solely computational approach to the air–water system such as that recently put forward by Pizzo, Deike & Ayet (2021) and Wu, Popinet & Deike (2022) is still beyond practical use. Hence the need for analytical models, such as that described here, even though its parameters will need careful examination.