2.1. Problem definition

In the early stage of the wind-wave generation process, the surface elevation is not substantially high, such that nonlinear wave-wave interactions can be ignored (Phillips Reference Phillips1957). Therefore, the linearised water wave equation is adequate to model the evolution of surface elevation. The linearised governing equation for the surface wave elevation $\eta$ can be expressed in wavenumber space as follows:

(2.1) \begin{align} \frac {\partial ^2\hat \eta }{\partial t^2}+\Lambda ^2(\boldsymbol {k})\hat \eta =-\frac {k}{\rho _w}\hat {p}(\boldsymbol {k},t), \end{align}

where $(\hat \cdot )$ represents the Fourier transform. Here, $\hat \eta (\boldsymbol {k},t)$ and $\hat p(\boldsymbol {k},t)$ denote the complex-valued wave amplitude and pressure, respectively, at the water’s surface for wavenumber $\boldsymbol {k}$ and time $t$ . The angular frequency $\Lambda (\boldsymbol {k})$ indicates the dispersion of the surface waves and is $\Lambda (\boldsymbol {k})=\sqrt {g|\boldsymbol {k}|}$ for gravity waves and $\Lambda (\boldsymbol {k})=\sqrt {g|\boldsymbol {k}|+(\sigma /\rho _w)|\boldsymbol {k}|^3}$ for gravity–capillary waves. Here, $g$ denotes gravitational acceleration, $\sigma$ denotes the surface tension of the air–water interface and $\rho _w$ denotes water density. Hereafter, $\rho_w$ is incorporated into $\hat{p}(\boldsymbol {k},t)$ for simplicity.

In the present study, the wind-induced water shear current is not considered. The influence of water shear flow on wave growth has been investigated in several studies. Nové-Josserand et al. (Reference Nové-Josserand, Perrard, Lozano-Durán, Benzaquen, Rabaud and Moisy2020) demonstrated that the presence of longitudinal water currents generates fine-scale surface wrinkles. Geva & Shemer (Reference Geva and Shemer2022) derived coupled Orr–Sommerfeld equations to include the effects of water shear flow. In the present study, we assume the water is initially calm and the air–water interface is flat, consistent with the framework established by Phillips (Reference Phillips1957). While wind shear can induce drift velocities in the water as waves develop, previous numerical studies have shown that, during the principal stage of wind-wave generation starting from calm water, the amplitude of the drift current velocity is small compared with the wave phase speed and can be neglected (Li & Shen Reference Li and Shen2022 b). Our present study focuses on the initial stage of wave generation, which precedes the principal stage in Phillips’ theory. At this stage, the assumption of zero surface drift velocity is reasonable. The assumption of negligible surface drift velocity has also been adopted in other literatures with a focus on wave growth when wind-induced water currents are not strong (e.g. Perrard et al. Reference Perrard, Lozano-Durán, Rabaud, Benzaquen and Moisy2019; Maleewong & Grimshaw Reference Maleewong and Grimshaw2022; Zhang, Hector & Moisy Reference Zhang, Hector, Rabaud and Moisy2023). The scenario of appreciable shear currents studied by Nové-Josserand et al. (Reference Nové-Josserand, Perrard, Lozano-Durán, Benzaquen, Rabaud and Moisy2020) and Geva & Shemer (Reference Geva and Shemer2022) is also an important problem, but is beyond the scope of the present work and will be a subject of our future research. In the present study, water dissipation is not considered, which is a mild assumption and is consistent with Phillips’ theory. During the initial stage of wind-wave generation, when waves are excited from a calm water surface, the wind input dominates and overcomes water dissipation. The effect of water viscosity can be incorporated into the governing equation (2.1) by introducing an additional term $4\nu k^2 \partial _t \eta$ (Perrard et al. Reference Perrard, Lozano-Durán, Rabaud, Benzaquen and Moisy2019), where $\nu$ is water kinematic viscosity. Our future research will address wind-wave generation considering the effects of water viscosity, which is particularly relevant in the presence of sheared water currents.

In the theory presented in Phillips (Reference Phillips1957), the temporal evolution of air pressure fluctuations is described according to Taylor’s frozen hypothesis (Taylor Reference Taylor1938), which can be formulated as

(2.2) \begin{align} \hat {p}(\boldsymbol {k},t)=\hat {p}_0(\boldsymbol {k})\exp (- \mathrm {i}\boldsymbol {k}\cdot \boldsymbol {U} t), \end{align}

where $\hat {p}_0(\boldsymbol {k})$ denotes the air pressure fluctuations at time $t=0$ , and $\boldsymbol {U}$ denotes the convection velocity of the air pressure fluctuations. In the present study, we assume that $\boldsymbol {U}$ is a constant vector pointing in the $x$ -direction with a magnitude of $U$ . The convection velocity $U$ in our study represents the convection of pressure fluctuations at the air–water interface. As demonstrated by Li & Shen (Reference Li and Shen2022 b, figure 8), the convection velocity of these air pressure fluctuations at the water’s surface exhibits only slight variations with respect to the wavenumbers over the range of scales relevant to wind-wave generation. Therefore, assuming a constant convection velocity of air pressure fluctuations at the air–water interface across different wavenumbers is a reasonable simplification for this analysis. This assumption aligns with previous theoretical studies of wind-wave interactions (e.g. Perrard et al. (Reference Perrard, Lozano-Durán, Rabaud, Benzaquen and Moisy2019)). Under this assumption, the governing equation for $\hat \eta$ becomes

(2.3) \begin{align} \frac {\partial ^2\hat \eta }{\partial t^2}+\Lambda ^2(\boldsymbol {k})\hat \eta =-k\hat {p}_0(\boldsymbol {k})\exp (- \mathrm {i}\boldsymbol {k}\cdot \boldsymbol {U} t). \end{align}

The assumption that turbulent pressure fluctuations retain coherence over short durations is based on Taylor’s frozen hypothesis, which is applicable in the initial stage of wind-wave generation in Phillips’ theory. This hypothesis assumes that turbulence structures are convected by the mean flow without significant temporal evolution, allowing us to treat the pressure fluctuations as coherent over the relevant time scales. Choi & Moin (Reference Choi and Moin1990) evaluated Taylor’s frozen turbulence hypothesis, demonstrating its overall effectiveness, particularly for large turbulence structures. This convective nature of turbulence pressure fluctuations has been extensively studied in the turbulence literature (e.g. Choi & Moin Reference Choi and Moin1990; Luhar, Sharma & MeKoen Reference Luhar, Sharma and McKeon2014; Yang & Yang Reference Yang and Yang2022). Given the initial conditions of $\hat \eta (\boldsymbol {k},t=0)=0$ and $\partial _t\hat \eta (\boldsymbol {k} ,t=0)=0$ , the solution to $\hat \eta$ is calculated as

(2.4) \begin{align} \hat \eta (\boldsymbol {k},t)=\frac {k \hat {p}_0[-\mathrm {i}\sin (\Lambda t)\boldsymbol {k}\cdot \boldsymbol {U}+\Lambda \cos (\Lambda t)-\Lambda \exp (-\mathrm {i}\boldsymbol {k}\cdot \boldsymbol {U} t) ]}{\Lambda [-(\boldsymbol {k}\cdot \boldsymbol {U})^2+\Lambda ^2]}, \end{align}

where $k=|\boldsymbol {k}|$ is the magnitude of the wavenumber vector $\boldsymbol {k}$ . The wave energy component $|\hat \eta |^2$ has the following explicit expression:

(2.5) \begin{align} |\hat \eta |^2(\boldsymbol {k},t) & = \frac {\Phi _p k^2}{\Lambda ^2[(\boldsymbol {k}\cdot \boldsymbol {U})^2-\Lambda ^2]^2}[(\boldsymbol {k}\cdot \boldsymbol {U})^2+\Lambda ^2-2\Lambda ^2\cos (\Lambda t)\cos (\boldsymbol {k}\cdot \boldsymbol {U} t)\notag \\ & \quad -2\Lambda \sin (\Lambda t)\boldsymbol {k}\cdot \boldsymbol {U}\sin (\boldsymbol {k}\cdot \boldsymbol {U} t)-((\boldsymbol {k}\cdot \boldsymbol {U})^2-\Lambda ^2)\cos ^2(\Lambda t)], \end{align}

where $\Phi _p=|\hat {p}_0(\boldsymbol {k})|^2$ is the energy spectrum of air pressure fluctuations at the air–water interface. The surface elevation variance is obtained from the integral of $|\hat \eta |^2(\boldsymbol {k},t)$ over the wavenumber space as follows:

(2.6) \begin{align} \langle \eta ^2\rangle (t)= \iint |\hat \eta |^2(\boldsymbol {k},t) \mathrm {d}\boldsymbol {k}. \end{align}

The goal of the present study is to understand the temporal behaviours of the wave energy, i.e. the surface elevation variance $\langle \eta ^2\rangle (t)$ . Next, we review the original Phillips approach for analysing wave energy growth, discuss its limitations and introduce our new approach based on complex analyses.

2.2. Review of the Phillips approach and its limitations

In this section, we present a brief overview of the Phillips approach for analysing the evolution of wave energy (i.e. (2.6)). The wave energy component increases quadratically over time when it is located on the resonance curve in the spectral space, defined as $\boldsymbol {k}\cdot \boldsymbol {U}-\Lambda (\boldsymbol {k})=0$ . However, the wave growth rate decreases when the wavenumber deviates from this resonance curve. To investigate the temporal behaviour of the wave energy, Phillips (Reference Phillips1957) utilised a new orthogonal curvilinear coordinate system that is associated with the resonance curve $\boldsymbol {k}\cdot \boldsymbol {U}-\Lambda (\boldsymbol {k})=0$ . Figure 1(a) shows the resonance curves for gravity waves and gravity–capillary waves. Phillips introduced a new local orthogonal coordinate system $(\chi ^t,s)$ , where $\chi ^t= (\boldsymbol {k}\cdot \boldsymbol {U}-\Lambda (\boldsymbol {k}) )t$ represents the non-dimensional distance to the resonance curve and $s$ denotes the distance along the curve $\chi ^t=0$ from the intersection of this curve with the $k_x$ -axis. The coordinate system $(\chi ^t,s)$ is illustrated using red arrows in figure 1(b). The integral of the wave energy components over the two-dimensional wavenumber space can hence be calculated in this new $(\chi ^t,s)$ coordinate system by incorporating Jacobian matrices. Phillips (Reference Phillips1957) aimed to demonstrate that the integral of the wave energy components along $\chi ^t$ between $\chi ^t=-2\pi$ and $\chi ^t=2\pi$ is proportional to time, resulting in a linear increase in the wave energy. However, this approach has several limitations and lacks rigour in the following aspects.

1. (i) First, the $(\chi ^t,s)$ coordinate system is not explicitly defined; thus deriving the Jacobian matrix that connects the $(k_x,k_y)$ and $(\chi ^t,s)$ coordinate systems at any arbitrary location is a non-trivial task.

2. (ii) Second, in Phillips’ analyses, the Jacobian matrix was approximated using its value on the resonance curve ( $\chi ^t=0$ ) for various $\chi ^t$ and $s$ values. Phillips (Reference Phillips1957) refers to this approximation as “simple.” This approximation did not lead to a quantitative result.

3. (iii) Third, the integration along $\chi ^t$ was numerically restricted to the interval $(-2\pi ,2\pi )$ , as depicted using the green dashed arrow in figure 1(b). With this constraint, the contributions of the wave energy components outside this range to the total wave energy were neglected. In other words, the calculations were not performed across the entire wavenumber space.

Hence, Phillips (Reference Phillips1957) only pointed out that the integration of the wave energy components along $\chi ^t$ is a linear function of time under the above approximations, but he did not obtain a quantification of the temporal behaviour of the wave energy with rigorous mathematical proof. In his original paper, Phillips described his method as “a simple and useful approximation to the integral”. The fundamental reason why the Phillips approach cannot yield more rigorous results lies in the choice of the resonance curve-based orthogonal curvilinear coordinate system $(\chi ^t,s)$ . However, we find that this curvilinear coordinate system is not necessary for wave evolution analyses. In the present study, we show that the integral of the wave energy components can be directly analysed in a Cartesian coordinate system $(k_x, k_y)$ using a complex analysis method, as described in § 2.3. Based on this new mathematical formulation, we develop a concrete theoretical framework for the wind-wave generation mechanism in the initial stage for gravity waves and gravity–capillary waves, as discussed in §§ 3 and 4, respectively, with the flow physics elucidated.

Figure 1. Panel (a) shows the resonance curves for gravity waves and gravity–capillary waves to illustrate their difference. Panel (b) illustrates the Phillips (Reference Phillips1957) approach for calculating $\langle \eta ^2\rangle$ based on the $(\chi ^t,s)$ coordinate system in the wavenumber space. The content is a rework of figure 2 in Phillips (Reference Phillips1957).

2.3. Our new approach

We find that an analysis performed directly using the Cartesian wavenumber coordinates $(k_x,k_y)$ is feasible for analytically obtaining the linear increase in the leading-order term of the upper-bound solution of the surface elevation variance over time. This approach differs from that of Phillips (Reference Phillips1957). We present the main mathematical framework and results of the present study in this section. Extensive mathematical analyses and investigations into the flow physics are provided in §§ 3 and 4.

We evaluate the temporal increase in the surface elevation variance. Equation (2.6), which describes the surface elevation variance $\langle \eta ^2\rangle$ , can be calculated by integrating first over variable $k_x$ and then over $k_y$ . Let function $E(k_y,t)$ be the integral of the wave energy component $|\hat \eta |^2$ along $k_x=(-\infty ,\infty )$ , i.e.

(2.7) \begin{align} E(k_y,t)=\int _{-\infty }^{\infty } |\hat \eta |^2(k_x,k_y,t) \mathrm {d}k_x. \end{align}

Therefore, the surface elevation variance $\langle \eta ^2\rangle$ can be expressed as

(2.8) \begin{align} \langle \eta ^2\rangle = \int _{-\infty }^{\infty } E(k_y,t)\mathrm {d} k_y. \end{align}

In the present study, we prove that for any $k_y\in \mathbb {R}$ , the leading-order term of the upper-bound solution of $E(k_y,t)$ is a linear function of time $t$ .

We rewrite the wave energy component $|\hat \eta |^2(k_x, k_y, t)$ in (2.5) as

(2.9) \begin{align} |\hat \eta |^2(k_x,k_y,t) & = \frac {\Phi _p (k_x^2+k_y^2)}{\Lambda ^2(B^2-\Lambda ^2)^2} \left \{\frac {1}{2}B^2+\frac {3}{2}\Lambda ^2-\Lambda (\Lambda +B)\cos [(\Lambda -B)t]\right .\notag \\ &\quad \left .-\Lambda (\Lambda -B)\cos [(\Lambda +B)t]-\frac {1}{2}(B^2-\Lambda ^2)\cos (2\Lambda t)\right \}, \end{align}

where $B=\boldsymbol {k}\cdot \boldsymbol {U}=k_xU$ . Because the pressure spectrum $\Phi _p$ is generally not an analytic function, we use the following inequality to move the term $\Phi _p$ outside of the integral in (2.7). Combining (2.7) and (2.9), we obtain the following estimation:

(2.10) \begin{align} E(k_y,t)\leqslant \Vert \Phi _p (k_x^2+k_y^2)\Lambda ^{-3} \Vert _{L^\infty _{x}}\int _{-\infty }^{\infty } F(k_x, k_y,t) \mathrm {d}k_x, \end{align}

where $F(k_x,k_y,t)$ is constructed as

(2.11) \begin{align} F(k_x,k_y,t)=&\frac {\Lambda }{(B^2-\Lambda ^2)^2} \left \{\frac {1}{2}B^2+\frac {3}{2}\Lambda ^2-\Lambda (\Lambda +B)\cos [(B-\Lambda )t]\right .\notag \\ &\left .-\Lambda (\Lambda -B)\cos [(\Lambda +B)t]-\frac {1}{2}(B^2-\Lambda ^2)\cos (2\Lambda t)\right \}. \end{align}

Here, the notation $\Vert \cdot \Vert _{L_x^\infty }$ denotes the $L^\infty$ -norm with respect to $x$ . For a general function $f(x,y)$ , we have

(2.12) \begin{align} \Vert f \Vert _{L_x^\infty }&=\max _{x\in \mathbb {R}} |f(x,y)|. \end{align}

The inequality (2.10) always holds due to the positivity of function $|\hat \eta |^2(k_x,k_y,t)$ . The pressure spectrum $\Phi _p$ can be considered a function that decays quickly when the wavenumber $\boldsymbol {k}$ approaches zero or infinity in the wavenumber space. Therefore, the $L^\infty$ norm of $\Phi _p (k_x^2+k_y^2)\Lambda ^{-3}$ has a finite value. We choose a specific expression for $F(k_x,k_y,t)$ to ensure that the improper integral of $F(k_x,k_y,t)$ with respect to $k_x$ on the right-hand side of (2.10) also has a finite value. For gravity waves, $\Lambda$ and $B$ behave as $\Lambda \sim |\boldsymbol {k}|^{1/2}$ and $B\sim k_x$ , respectively. For gravity–capillary waves, $\Lambda$ and $B$ behave as $\Lambda \sim |\boldsymbol {k}|^{3/2}$ and $B\sim k_x$ , respectively, when $k_x$ approaches infinity (i.e. $k_x\rightarrow \infty$ ). Therefore, for any wavenumber $k_y$ , the function $F$ behaves as $k_x^{-3/2}$ when $k_x\rightarrow \infty$ for both gravity waves and gravity–capillary waves. The $-3/2$ power-law decay ensures that the integral is finite. When $k_y=0$ , $F$ behaves as $F\sim k_x^{-1/2}$ . The integral is also finite when the endpoint of the integration interval is $k_x=0$ . Therefore, the right-hand side of (2.10) has a finite value.

Because $F(k_x,k_y,t)$ is an even function with respect to $k_x$ , we can limit our focus to the integral along the positive $x$ -axis, i.e. the interval $(0,\infty )$ . Hence, we define function $G(k_y,t)$ as

(2.13) \begin{align} G(k_y,t)=\int _0^\infty F(k_x,k_y,t) \mathrm {d}k_x. \end{align}

The function $G(k_y,t)$ represents half of the integral of $F(k_x,k_y,t)$ with respect to $k_x$ over the range $(-\infty ,\infty )$ .

To understand the temporal evolution of the surface elevation variance, we should examine the right-hand side of (2.10), especially the function $G(k_y,t)$ in (2.13). Next, we prove that for gravity waves, the function $G(k_y,t)$ is bounded by a linear function of time plus an oscillatory component with an amplitude that is restricted by a constant. The linear increase in $G(k_y,t)$ suggests that the surface elevation variance $\langle \eta ^2\rangle (t)$ increases linearly with time $t$ . However, the linear increase in the surface elevation variance for gravity–capillary waves depends on whether the convection velocity of the air pressure fluctuations at the water’s surface exceeds a certain threshold $U_{{min}}$ . Specifically, for gravity–capillary waves, when the convection velocity of the air pressure fluctuations at the water’s surface, $U$ , exceeds $U_{{min}}$ , the function $G(k_y,t)$ is restricted by a linear function of time plus an oscillatory component with an amplitude that is limited by a constant. Consequently, the surface elevation variance $\langle \eta ^2\rangle$ exhibits linear growth over time. However, when the magnitude of the convection velocity of the air pressure fluctuations, $U$ , is less than $U_{{min}}$ , the function $G(k_y,t)$ is restricted by a constant that is independent of time $t$ , which suggests that waves do not grow under such weak wind conditions. The detailed proof is presented in § 3 for gravity waves and § 4 for gravity–capillary waves.

3.1. Properties of the resonance curve for gravity waves

The dispersion relation of gravity surface waves relates the angular frequency $\Lambda$ of the waves to their wavenumber $\boldsymbol {k}$ and is given by $\Lambda =\sqrt {g|\boldsymbol {k}|}$ , where $|\boldsymbol {k}|=\sqrt {k_x^2+k_y^2}$ and $g$ represents gravitational acceleration. The denominator in (2.11) is singular when the following resonance condition is satisfied: $B-\Lambda =0$ . For a chosen $k_y\in \mathbb {R}$ , the value of $k_x$ such that the point $(k_x,k_y)$ is located on the resonance curve $\chi =0$ in the wavenumber space must satisfy the following condition:

(3.1) \begin{align} k_x=\pm \frac {g}{U^2}\sqrt {\frac {1}{2}\left (1+\sqrt {1+\frac {4k_y^2U^4}{g^2}}\right )}. \end{align}

This implies that, for any $k_y$ , there exists a positive value of $k_x$ that satisfies the resonance condition $\Lambda - \boldsymbol {k} \cdot U = 0$ . In the first quadrant of the wavenumber space, the resonance curve $\chi$ begins at $(g/U^2,0)$ on the $x$ -axis and extends to infinity, as depicted in figure 1( $a$ ).

3.2. Calculation of $G(k_y,t)$

Next, we calculate the function $G(k_y,t)$ defined in (2.13). The key aspect of our approach is that we employ the residue theorem to calculate the integral in the definition of $G(k_y,t)$ , as shown in (2.13). First, we perform an analytical continuation to convert the function $F(k_x,k_y,t)$ in the integrand on the right-hand side of (2.13) into a complex function with good properties. We calculate the leading-order term using the residue theorem, resulting in a linear function of time. The oscillatory term is obtained based on complex analyses.

Figure 2. Illustration of the indented closed loop for integration, as outlined in (3.4).

3.2.1. Analytical continuation of $F$ in the complex plane

In this subsection, we present a novel approach to evaluate the integral of $F$ using complex analyses. We transform the real-valued function $F$ into a complex-valued function in the complex $z$ -plane. Hereafter, $z$ represents a complex number. We examine a related complex function $F_1(z,k_y,t)$ , which is obtained by analytically continuing the function $F(k_x,k_y,t)$ into the complex $z$ -plane as follows:

(3.2) \begin{align} F_1(z,k_y,t) & = \frac {\Lambda }{(B^2-\Lambda ^2)^2} \left \{\frac {1}{2}B^2+\frac {3}{2}\Lambda ^2-\Lambda (\Lambda +B)\exp [\mathrm {i}(B-\Lambda )t]\right .\notag \\ & \quad \left .-\Lambda (\Lambda -B)\exp [\mathrm {i}(\Lambda +B)t]-\frac {1}{2}(B^2-\Lambda ^2)\exp (2\mathrm {i}\Lambda t)\right \}. \end{align}

Here, $B=Uz$ and $\Lambda =\sqrt {g\sqrt {z^2+k_y^2}}$ , as we have substituted the variable $k_x$ with the complex-valued variable $z$ . It is evident that $F(k_x,k_y,t)$ is equal to the real part of $F_1(k_x,k_y,t)$ for any real values of $k_x$ , $k_y$ and $t$ . The function $F_1$ has a pole at $z=z_*$ on the positive real axis, where $B - \Lambda = 0$

(3.3) \begin{align} z_* =\frac {g}{U^2}\sqrt {\frac {1}{2}\left (1+\sqrt {1+\frac {4k_y^2U^4}{g^2}}\right )}. \end{align}

We consider a closed loop $\Gamma$ in the complex plane, as depicted in figure 2, which is defined as $\Gamma =[0, z_*-r]+C_r+[z_*+r, R]+C_R+L_R$ . Here, $C_r( {\theta })=z_*+re^{\mathrm {i}(\pi -\theta )}$ , with $0\leqslant \theta \leqslant \pi$ a clockwise circular path around $z_*$ with a radius of $r$ ; $C_R(\theta )=Re^{\mathrm {i}\theta }$ , with $0\leqslant \theta \leqslant \pi /4$ , is an anticlockwise circular path around 0 with a radius of $R$ ; and $L_R(\theta )=Re^{\mathrm {i}\pi /4}(1-\theta )$ , with $0\leqslant \theta \leqslant 1$ , is a straight path pointing towards the origin. We choose the angle between the real axis and the ray $-L_R$ to be $\pi /4$ . However, the choice of this angle is arbitrary, and the real part of the integral of the function $F_1$ along $L_R$ is not dependent on this angle. This is ensured by the fact that the contribution of intergral along $C_R$ vanishes as $R\rightarrow \infty$ regardless the angle between $L_R$ and the $x$ -axis, as further discussed in § 3.2.2. The independence of the angle can also be confirmed through numerical integration. Because the function $F_1$ has no singularities inside the closed loop $\Gamma$ , the integral of $F_1$ with respect to $z$ along the contour $\Gamma$ is equal to zero. This can be written as

(3.4) \begin{align} \left (\int _0^{z_*-r}+\int _{C_r}+\int _{z_*+r}^R+\int _{C_R}+\int _{L_R}\right )F_1(z,k_y,t)\mathrm {d}z=0. \end{align}

Therefore, the integration of $F_1$ along the real positive axis can be calculated as

(3.5) \begin{align} \int _0^{\infty }F_1(k_x,k_y,t) \mathrm {d}k_x=&\lim _{R\rightarrow \infty }\int _{0}^{R}F_1(k_x,k_y,t) \mathrm {d}k_x\notag \\ =&\lim _{R\rightarrow \infty }\lim _{r\rightarrow 0}\mathrm {Re}\left [\left (\int _0^{z_*-r}+\int _{z_*+r}^R\right )F_1(z,k_y,t)\mathrm {d}z\right ]\notag \\ =&\lim _{R\rightarrow \infty }\lim _{r\rightarrow 0}\mathrm {Re}\left [\left (-\int _{C_r}-\int _{C_R}-\int _{L_r}\right )F_1(z,k_y,t)\mathrm {d}z\right ]\notag \\ =&-\lim _{R\rightarrow \infty }\mathrm {Re}\left [\int _{C_R}F_1(z,k_y,t) \mathrm {d}z\right ] -\lim _{r\rightarrow 0}\mathrm {Re}\left [\int _{C_{r}}F_1(z,k_y,t) \mathrm {d}z\right ]\notag \\ &-\lim _{R\rightarrow \infty }\mathrm {Re}\bigg [\int _{C_{L_R}}F_1(z,k_y,t) \mathrm {d}z\bigg]. \end{align}

3.2.2. Calculation of the integral of $F_1$ along the path $C_R$

We evaluate $F_1$ along the path $C_R$ . Upon examining the expression of $F_1$ in (3.2), we observe that the essential properties of $F_1$ are governed by three exponential functions, namely, $\exp [\mathrm {i}(B-\Lambda )t]$ , $\exp [\mathrm {i}(B+\Lambda )t]$ and $\exp (2\mathrm {i}\Lambda t)$ . Assuming $z = R\exp (\mathrm {i}\theta )$ with $0\leqslant \theta \leqslant \pi /4$ , for sufficiently large values of $R$ , we have the following asymptotic expansion:

(3.6) \begin{align} B-\Lambda &=UR\mathrm {e}^{\mathrm {i}\theta }-\sqrt {g}( R^2 \mathrm {e}^{2\mathrm {i}\theta }+k_y^2)^{\frac {1 }{4}} \nonumber\\ &=UR\mathrm {e}^{\mathrm {i}\theta }-\sqrt {gR} \mathrm {e}^{\frac {\mathrm {i}\theta }{2}}+O(R^{-\frac {3}{2}}). \end{align}

Therefore, for sufficiently large values of $R$ , we obtain the following expression for $B-\Lambda$ :

(3.7) \begin{align} B-\Lambda =\alpha _1+\mathrm {i}\beta _1, \end{align}

where $\alpha _1$ and $\beta _1$ are real positive functions of $R$ . The leading-order approximation shows that $\alpha _1\approx UR\cos (\theta )$ and $\beta _1\approx UR\sin (\theta )$ for sufficiently large values of $R$ . Hence, the term $\exp [\mathrm {i}(B-\Lambda )t]$ becomes

(3.8) \begin{align} \exp [\mathrm {i}(B-\Lambda )t]=\exp (-\beta _1 t+\mathrm {i}\alpha _1 t). \end{align}

This analysis demonstrates that the real part of $\exp [\mathrm {i}(B-\Lambda )t]$ is a bounded function and it decays as $R$ or $t$ approaches infinity.

Using a similar method, we next study the properties of the function $\exp [\mathrm {i}(B+\Lambda )t]$ . Because $0\leqslant \theta \leqslant \pi /4$ , we have the following expansion for $B+\Lambda$ :

(3.9) \begin{align} B+\Lambda =\alpha _2+\mathrm {i}\beta _2, \end{align}

where $\alpha _2$ and $\beta _2$ are both real positive functions of $R$ , and the leading-order approximation shows that $\alpha _1\approx UR\cos (\theta )$ and $\beta _1\approx UR\sin (\theta )$ for sufficiently large $R$ . The term $\exp [\mathrm {i}(B+\Lambda )t]$ becomes

(3.10) \begin{align} \exp [\mathrm {i}(B+\Lambda )t]=\exp (-\beta _2 t-\mathrm {i}\alpha _2 t). \end{align}

This result indicates that the real part of $\exp [\mathrm {i}(B+\Lambda )t]$ is a bounded function and it decays as $R$ or $t$ approaches infinity.

For the term $\exp (2\mathrm {i}\Lambda t)$ , we calculate the following asymptotic expansion for a sufficiently large value of $R$ :

(3.11) \begin{align} 2\mathrm {i}\Lambda =\alpha _3+\mathrm {i}\beta _3, \end{align}

where $\alpha _3$ and $\beta _3$ are positive functions of $R$ that satisfy $\alpha _3\approx 2\sqrt {gR}\cos (\theta /2)$ and $\beta _3\approx 2\sqrt {gR}\sin (\theta /2)$ , respectively, for sufficiently large values of $R$ . Hence, the term $\exp (2\mathrm {i}\Lambda t)$ becomes

(3.12) \begin{align} \exp (2\mathrm {i}\Lambda t)=\exp (-\beta _3 t-\mathrm {i}\alpha _3 t), \end{align}

of which the real part is also a bounded function of the variables $R$ and $t$ .

The above analyses suggest that, along the curve $C_R$ ,

(3.13) \begin{align} |F_1|\lt \frac {C}{R^2}, \end{align}

where $C$ is a function of $k_y$ and $t$ and does not depend on $R$ . Therefore, we obtain the following estimation:

(3.14) \begin{align} 0\leqslant \left |\lim _{R\rightarrow \infty }\int _{C_R}F_1 \mathrm {d}z\right | \leqslant \lim _{R\rightarrow \infty }\int _{C_R}\left |F_1(z)\right | \mathrm {d}z \leqslant \lim _{R\rightarrow \infty } \int _{C_R} \frac {C}{R^2}\mathrm {d}z=0. \end{align}

We thus conclude that

(3.15) \begin{align} \lim _{R\rightarrow \infty }\int _{C_R}F_1(z) \mathrm {d}z=0. \end{align}

3.2.3. Calculation of the integral of $F_1$ along the path $C_r$

The integral of $F_1$ along the semicircular path $C_r$ can be calculated using the residue theorem, from which we have

(3.16) \begin{align} \lim _{r\rightarrow 0}\int _{C_{r}}F_1(z) \mathrm {d}z&=-\pi \mathrm{iRes}(F_1,z_{*}). \end{align}

Here, the notation $\mathrm {Res}(F_1,z_{*})$ denotes the residue of the function $F_1$ at the point $z=z_*$ . The function $F_1$ has a first-order pole at $z = z_*$ . As a result, we have

(3.17) \begin{align} \mathrm {Res}(F_1,z_*)=\lim _{z\rightarrow z_*} F_1(z)(z-z_*). \end{align}

To determine the residue of the function $F_1$ in the above equation, we employ a perturbation method to calculate the limit near $z=z_*$ . We first perform a Taylor expansion of the term $B-\Lambda$ around $z=z_*$ and obtain

(3.18) \begin{align} B-\Lambda =\frac {(z_*^2+2k_y^2)U}{2z_*^2+2k_y^2}(z-z_*)+O((z-z_*)^2). \end{align}

The same procedure is applied to both the denominator and numerator of the function $F_1(z,k_y,t)$ near $z=z_*$ . For simplicity, we define a dimensionless local Froude number $\gamma _y$ based on the wavenumber $k_y$ and convection velocity of the pressure fluctuations $U$ as

(3.19) \begin{align} \gamma _y=\frac {U^2 k_y}{g}. \end{align}

Hence, we obtain that

(3.20) \begin{align} {\lim _{r\rightarrow 0}\int _{C_{r}}F_1(z,k_y,t) \mathrm {d}z=-\frac {\pi g \big(1+\sqrt {4\gamma _y^2+1}\,\big)^{\frac {3}{2}}}{2\sqrt {2}\sqrt {4\gamma _y^2+1}U^2}t}. \end{align}

Equation (3.20) indicates that this component is linearly proportional to time $t$ .

3.2.4. Calculation of the integral of $F_1$ along the path $L_R$

Next, we consider the integral of $F_1$ along path $L_R$ , where the variable $z$ can be parameterised as a function of a real variable $r$ , such that $z(r) = r\exp (\mathrm {i}\pi /4)$ with $0\leqslant r\leqslant R$ . Therefore, we have

(3.21) \begin{align} \int _{L_R} F_1(z)\mathrm {d}z&=\exp \left (\frac {\mathrm {i}\pi }{4}\right )\int _R^0 F_1(z(r))\mathrm {d}r\notag \\ &=-\frac {1+\mathrm {i}}{\sqrt {2}}\int _0^R\frac {K(r,t)}{J(r)}\mathrm {d}r. \end{align}

In (3.21), the numerator $K(r,t)$ and denominator $J(r)$ are

(3.22) \begin{align} K(r,t) & = \Lambda _0(r)\left \{\frac {1}{2}B_0(r)^2+\frac {3}{2}\Lambda _0(r)^2-\Lambda _0(r)[\Lambda _0(r)+B_0(r)]\exp [\mathrm {i}(B_0(r)-\Lambda _0(r))t]\right .\notag \\ & \quad -\Lambda _0(r)[\Lambda _0(r)-B_0(r)]\exp [\mathrm {i}(B_0(r)+\Lambda _0(r))t]\notag \\ & \quad \left .-\frac {1}{2}[B_0(r)^2-\Lambda _0(r)^2]\exp [2\mathrm {i}\Lambda _0(r) t]\right \}, \end{align}

and

(3.23) \begin{align} J(r)=&\big[B_0(r)^2-\Lambda _0(r)^2\big]^2, \end{align}

respectively. Here, the complex-valued variables $B_0(r)$ and $\Lambda _0(r)$ are functions of the real-valued number $r$ . They can be written as

(3.24) \begin{align} B_0(r)&=\frac {1+\mathrm {i}}{\sqrt {2}}Ur, \end{align}

(3.25) \begin{align} \Lambda _0(r)&=\sqrt {g\sqrt {\mathrm {i}r^2+k_y^2}}. \end{align}

We aim to examine the time evolution behaviour of $F_1(z)$ for any complex variable $z$ located along the path $L_R$ . For large values of $k_y$ , $\mathrm {Im}(B_0(r) - \Lambda _0(r)) \gt 0$ always holds. We can split the complex-valued term $B_0(r) - \Lambda _0(r)$ into a real part $\alpha _4(r)$ and an imaginary part $\beta _4(r)$

(3.26) \begin{align} B_0(r)-\Lambda _0(r)=\alpha _4(r)+i\beta _4(r). \end{align}

Here, $\alpha _4(r)$ and $\beta _4(r)$ are real functions and $\beta _4(r)\geqslant 0$ . Hence we have

(3.27) \begin{align} \exp [\mathrm {i}(B_0(r)-\Lambda _0(r))t]=\exp [\mathrm {i}\alpha _4(r)t]\exp [-\beta _4(r)t]. \end{align}

We observe that the magnitude of $\exp (\mathrm {i}(B_0(r) - \Lambda _0(r))t )$ is a bounded function of time $t$ for all complex variables $z(r)$ located along the path $L_R$

(3.28) \begin{align} \left |\exp \left [\mathrm {i}(B(r)-\Lambda (r))t\right ]\right |\leqslant 1. \end{align}

Similarly, we can obtain that

(3.29) \begin{align} \left |\exp [\mathrm {i}(B(r)+\Lambda (r))t]\right |&\leqslant 1,\end{align}

(3.30) \begin{align} \left |\exp [\mathrm {i}(2\Lambda (r))t]\right |&\leqslant 1. \end{align}

Therefore, by using the triangle inequality, the modulus of the complex-valued function $K(r)$ is bounded by

(3.31) \begin{align} |K(r,t)|\leqslant |B_0(r)|^2|\Lambda _0(r)|+4|\Lambda _0(r)|^3+2|B_0(r)\Lambda _0^2(r)|. \end{align}

According to (3.27), $K(r,t)$ has the following asymptotic behaviour for large values of $t$ :

(3.32) \begin{align} \lim _{t\rightarrow \infty }K(r,t)=\frac {1}{2}B_0^2(r) {\Lambda _0(r)}+\frac {3}{2}\Lambda _0^3(r). \end{align}

We note that the right-hand side of (3.31) is independent of time $t$ . The denominator $J(r)$ on the right-hand side of (3.21) is positive for all $z$ on $L_R$ , indicating that (3.21) does not have singularities along the path $L_R$ . Therefore, from (3.21), we can estimate its modulus as follows:

(3.33) \begin{align} \left |\int _{L_R}F_1(z)\mathrm {d}z\right | &\leqslant \int _0^\infty \frac {|B_0(r)|^2|\Lambda _0(r)|+4|\Lambda _0(r)|^3+2|B_0(r)\Lambda _0^2(r)|}{|B_0(r)^2-\Lambda _0(r)^2|^2}\mathrm {d}r. \end{align}

The integrand on the right-hand side of (3.33) has an asymptotic behaviour of $r^{-1.5}$ as $r$ approaches infinity. Consequently, the integration over the path $L_R$ yields a finite value that depends on the parameters $U$ , $g$ and the local spanwise wavenumber $k_y$ . When $k_y$ is small, the positivity of $\beta _4(r)$ cannot hold true for all values of $r\gt 0$ , although the decaying property of $B_0(r)-\Lambda _0(r)$ ensures that $\beta _4(r)$ is positive for large values of $r$ . As a result, under such circumstances, the pointwise estimation of (3.28) is not valid with respect to $r$ due to the potential exponential growth of time $t$ . However, the function $F(k_x,k_y,t)$ exhibits at most quadratic growth over time, and thus the exponential growth term should cancel out during the integration of $F_1$ along the path $L_R$ . In numerical examples in the next subsection, we confirm that the integral of $F_1$ along the path $L_R$ shows oscillatory behaviour with a decaying amplitude.

3.2.5. Properties of the function $G$

In this subsection, we summarise the calculations of $G(k_y,t)$ as defined in (2.13). Recall that integrating the function $F_1$ with respect to $z$ along the positive real axis yields the same result as the function $G(k_y,t)$ (see (3.5), (3.15), (3.20) and (3.21)). Then, we obtain the following expression for the function $G(k_y,t)$ :

(3.34) \begin{align} G(k_y,t)=M(k_y,t)+N(k_y,t), \end{align}

where the functions $M(k_y,t)$ and $N(k_y,t)$ are

(3.35) \begin{align} M(k_y,t)&= {\frac {\pi g\left (1+\sqrt {4\gamma _y^2+1}\right )^{\frac {3}{2}}}{2\sqrt {2}\sqrt {4\gamma _y^2+1}U^2}t}, \end{align}

(3.36) \begin{align} N(k_y,t)&=\mathrm {Re}\left [\int _0^\infty \frac {1+\mathrm {i}}{\sqrt {2}}\frac {K(r,t)} { {J}(r)} \mathrm {d}r\right ]. \end{align}

Here, $\gamma _y$ is the local Froude number defined in (3.19), $N(k_y,t)$ is the component that deviates from linear growth, and $J(r)$ is defined in (3.23). The function $N(k_y,t)$ oscillates with time $t$ , but its amplitude is bounded by a constant, as defined on the right-hand side of (3.33). Let $q(k_y)$ be the linear growth rate of $M(k_y,t)$

(3.37) \begin{align} q(k_y)=\frac {\pi g\left (1+\sqrt {4\gamma _y^2+1}\,\right )^{\frac {3}{2}}}{2\sqrt {2}\sqrt {4\gamma _y^2+1}U^2}. \end{align}

Because $\gamma _y$ scales as $\gamma _y\sim U^2$ , we obtain that the linear growth rate $q(k_y)$ has the following approximations for small $\gamma _y$ and large $\gamma _y$ :

(3.38) \begin{align} q(k_y)\approx \begin{cases} \dfrac {\pi g}{2\sqrt {2}U^2}, &\gamma _y\ll 1\\[10pt] \dfrac {\pi \sqrt {gk_y}}{2U}, &\gamma _y\gg 1 \end{cases}\ \end{align}

which shows that $q(k_y)$ scales as $q(k_y) \sim U^{-2}$ for small $\gamma _y$ and $q(k_y) \sim U^{-1}$ for large $\gamma _y$ . Additionally, as time $t$ approaches infinity, we have the following limit result:

(3.39) \begin{align} \lim _{t\rightarrow \infty }N(k_y,t)=\mathrm {Re}\left [ \int _0^\infty \frac {(1+\mathrm {i})\sqrt {g}\left (\mathrm {i}U^2r^2+3g\sqrt {\mathrm {i}r^2+k_y^2}\right )\big(\mathrm {i}r^2+k_y^2\big)^{\frac {1}{4}}}{4\left (\mathrm {i}U^2r^2-g\sqrt {\mathrm {i}r^2+k_y^2}\right )^2}\mathrm {d}r \right ]. \end{align}

Figure 3. Numerical examples of the calculation of $G(k_y,t)$ for representative values of the wavenumber $k_y$ , with results for $k_y=0\,\mathrm {m^{-1}}$ shown in the left column (panels a, c and e) and results for $k_y=2\,\mathrm {m^{-1}}$ shown in the right column (panels b, d and f). We set $g=9.8\,\mathrm {m\,s^{-2}}$ and $U=2\,\mathrm {m\,s^{-1}}$ . Panels (a) and (b) show $G(k_y,t)$ , with the blue lines indicating the results obtained from numerical integration, the red lines indicating the linear growth component, $M(k_y,t)$ and the black dots representing the sum of the linear growth component $M(k_y,t)$ and the oscillatory component $N(k_y,t)$ . In panels (c) and (d) the red lines show the oscillatory component $N(k_y,t)$ , and the blue lines show $G(k_y,t)-M(k_y,t)$ . Panels (e) and (f) show the numerical errors in the calculations, with the blue lines representing $G(k_y,t)-M(k_y,t)-N(k_y,t)$ .

We present numerical examples for calculating the function $G(k_y,t)$ for gravity waves, as shown in figure 3. To verify (3.34), we compute the function $G(k_y,t)$ by performing a numerical integration of (2.13) and compare the result with the sum of $M(k_y, t)$ and $N(k_y, t)$ . This comparison is illustrated in panels (a) and (b) of figure 3. The results confirm that the function $G(k_y, t)$ predominantly exhibits linear growth over time and that the linear function $M(k_y, t)$ perfectly captures the leading-order behaviour of $G(k_y, t)$ . The time evolution of the function $N(k_y,t)$ , which is the oscillatory component of $G(k_y,t)$ , is shown in panels (c) and (d) of figure 3. As $N(k_y,t)$ oscillates, its amplitude gradually decreases over time. This phenomenon supports the conclusion that the function $N(k_y,t)$ converges to a finite limit value as time approaches infinity, as demonstrated by (3.39). Panels (e) and (f) of figure 3 show the numerical error for the integration of $G(k_y,t)$ , calculated as $G(k_y,t)-M(k_y,t)-N(k_y,t)$ . The difference between $G(k_y,t)$ and $M(k_y,t)+N(k_y,t)$ remains less than $10^{-4}$ while $G(k_y,t)$ ranges from $0$ to $60$ , thereby validating the theoretical analyses.

3.3. Overall increases in the surface elevation variance $\langle \eta ^2\rangle$ over time

In this subsection, we demonstrate the evolution of the surface elevation variance $\langle \eta ^2\rangle$ based on (2.8). By combining (2.10), (2.13) and (3.34), we obtain the following estimation of the function $E(k_y,t)$ :

(3.40) \begin{align} E(k_y,t)\leqslant 2\Vert \Phi _p (k_x^2+k_y^2)\Lambda ^{-3} \Vert _{L^\infty _{x}}[M(k_y,t)+N(k_y,t)]. \end{align}

Hence, the surface elevation variance $\langle \eta ^2\rangle$ can be calculated from (2.8) as

(3.41) \begin{align} \langle \eta ^2\rangle & \leqslant 2\int _{-\infty }^\infty \Vert \Phi _p(k_x^2+k_y^2)\Lambda ^{-3}\Vert _{L_x^\infty } [M(k_y,t)+N(k_y,t)] \mathrm {d}k_y\notag \\ & \leqslant \Vert \Phi _p(k_x^2+k_y^2)\Lambda ^{-3}q(k_y) \Vert _{L_x^\infty L_y^1} t\notag \\ & \quad +2\int _{-\infty }^\infty \Vert \Phi _p(k_x^2+k_y^2)\Lambda ^{-3}\Vert _{L_x^\infty } N(k_y,t) \mathrm {d}k_y, \end{align}

where $\gamma _y$ and $q(\gamma _y)$ are defined in (3.19) and (3.37), respectively. Here, the notation of $\Vert \cdot \Vert _{L_x^\infty L_y^1}$ denotes the $L^\infty$ -norm with respect to $x$ and $L^1$ -norm with respect to $y$ . For a general function $f(x,y)$ , we have

(3.42) \begin{align} {\Vert f \Vert _{L_x^\infty L_y^1}=\int _{-\infty }^\infty \max _{x\in \mathbb {R}} |f(x,y)|\mathrm {d}y.} \end{align}

Equation (3.41) explicitly demonstrates the upper bound of the temporal increase in the energy of gravity waves; that is, the leading-order term increases linearly over time, and there also exists an oscillatory term representing the fluctuations in the surface elevation variance. The pressure fluctuation spectrum $\Phi _p$ is proportional to $u_\tau ^4$ , where $u_\tau$ is the air friction velocity (Phillips Reference Phillips1957). Considering the scaling of $q(k_y)$ with respect to $U$ (see (3.38)), we can presumably infer that the wave energy $\langle \eta ^2\rangle$ scales as $\langle \eta ^2\rangle \sim U^s$ with $s\in [2,3]$ during the initial stage of gravity wave growth.

In conclusion, we calculate the surface elevation variance $\langle \eta ^2\rangle$ for gravity waves forced by wind by integrating the wave energy component $|\hat \eta (\boldsymbol {k})|^2$ in the wavenumber space. We apply a complex analysis method to transform the integration along the $k_x$ axis to a singular integration within the complex $z$ -plane and resolve the singularity using the residue theorem. Our results indicate that, for gravity waves forced by wind, the leading-order term of the upper-bound solution of the surface elevation variance is a linear function of time, while the sub-leading term shows oscillatory behaviour.

4.1. Existence condition for resonance curves of gravity–capillary waves

We aim to identify the conditions under which the resonance mechanism exists for a given wavenumber $k_y$ . The analyses in § 3 demonstrate that the resonance condition always holds for any value of $k_y$ in the case of gravity waves. However, this is not the case for gravity–capillary waves.

We assume that the resonance condition is satisfied at wavenumber $\boldsymbol {k} = (k_x, k_y)$

(4.5) \begin{align} \boldsymbol {k} \cdot \boldsymbol {U}=\Lambda (\boldsymbol {k}). \end{align}

This expression is equivalent to

(4.6) \begin{align} k_x U=\sqrt {g\sqrt {k_x^2+k_y^2}+\frac {\sigma }{\rho _w}\left (\sqrt {k_x^2+k_y^2}\,\right )^3}. \end{align}

We define a dimensionless Froude number $\gamma _0$ , which is based on the pressure convection velocity $U$ , critical wavenumber $k_0$ and gravitational acceleration $g$ , as follows:

(4.7) \begin{align} \gamma _0=\frac {U^2 k_0}{g} {=\frac {2U^2}{c_{{min}}^2}}. \end{align}

By setting $\kappa _x=k_x/k_0$ and $\kappa _y=k_y/k_0$ , we can simplify (4.6) as

(4.8) \begin{align} \gamma _0 \kappa _x^2=\sqrt {\kappa _x^2+\kappa _y^2}+\left (\sqrt {\kappa _x^2+\kappa _y^2}\,\right )^3. \end{align}

Therefore, $\kappa _x^2$ is a positive root of the following cubic equation for $\xi$

(4.9) \begin{align} \xi ^3+\big(3\kappa _y^2+2-\gamma _0^2\big)\xi ^2+\big(1+\kappa _y^2\big)\big(1+3\kappa _y^2\big)\xi +\kappa _y^2\big(1+\kappa _y^2\big)^2=0. \end{align}

The coefficients in (4.9) are all real numbers, and we can obtain the properties of the roots using the discriminant. Let $\xi _i ( i=1,2,3)$ be the three roots of the cubic equation (4.9). We note that (4.9) has either three real roots or one real root and two complex conjugate roots.

We first consider the case when $\kappa _y\neq 0$ . According to Vieta’s formulas, we have

(4.10) \begin{align} \xi _1\xi _2\xi _3=-\kappa _y^2\big(1+\kappa _y^2\big)\lt 0. \end{align}

Therefore, when $\xi _1$ and $\xi _2$ are complex conjugate roots, the only real root $\xi _3$ is negative, i.e. $\xi _3\lt 0$ . When all the roots are real numbers, $\xi _1\geqslant \xi _2\gt 0$ and $\xi _3\lt 0$ . To determine the properties of the solutions, we next calculate the discriminant $\Delta _1$ for the cubic equation (4.9)

(4.11) \begin{align} \Delta _1 =& \big(1 + \kappa _y^2\big)^2\gamma _0^2 [4 \kappa _y^2 \gamma _0^4+ (1 - 18 \kappa _y^2 - 27 \kappa _y^4)\gamma _0^2 - 4]. \end{align}

The cubic equation (4.9) for $\xi$ has three real roots if and only if the discriminant satisfies $\Delta _1\geqslant 0$ . Equation (4.9) has one negative root and two complex conjugate roots when the discriminant satisfies $\Delta _1\lt 0$ .

Based on (4.11), we can infer that the sign of $\Delta _1$ is the same as that of $F_{\gamma }$

(4.12) \begin{align} F_{\gamma }\big(\gamma _0^2\big)=4 \kappa _y^2 \gamma _0^4+ \big(1 - 18 \kappa _y^2 - 27 \kappa _y^4\big)\gamma _0^2 - 4, \end{align}

which can be interpreted as a quadratic function of $\gamma _0^2$ . The discriminant of the quadratic function $F_{\gamma }$ is always positive; hence, there exist two real roots $\psi _1$ and $\psi _2$ such that $F_{\gamma }(\psi _1)=F_{\gamma }(\psi _2)=0$ . For simplicity, we assume that $\psi _1\lt \psi _2$ . Owing to Vieta’s formula, we have $\psi _1\psi _2=-4\lt 0$ , and we obtain $\psi _1\lt 0\lt \psi _2$ . Therefore, we conclude that $F_\gamma (\psi )\lt 0$ for any $\psi \in (0,\psi _2)$ and that $F_\gamma (\psi )\gt 0$ for any $\psi \in (\psi _2,\infty )$ . By direct calculations, we obtain that

(4.13) \begin{equation} \psi _2=\frac {-1+18\kappa _y^2+27\kappa _y^4+\sqrt {\big(1+\kappa _y^2\big)\big(1+9\kappa _y^2\big)^3}}{8\kappa _y^2}. \end{equation}

Hence, we have precise control over the sign of $\Delta _1$

(4.14) \begin{align} \Delta _1\lt 0,\quad & \mathrm {when}\, \gamma _0\lt \sqrt {\psi _2},\end{align}

(4.15) \begin{align} \Delta _1\gt 0,\quad & \mathrm {when}\, \gamma _0\gt \sqrt {\psi _2}. \end{align}

If the condition $\gamma _0\lt \sqrt {\psi _2}$ is satisfied, the cubic equation (4.9) for $\xi$ has only one negative real root $\xi _1$ . However, for $\xi _1$ to satisfy (4.8), it must also be non-negative, which is contradictory. As a result, for any $k_y$ , the resonance condition (4.6) has no solution when $\gamma _0\lt \sqrt {\psi _2}$ . When $\gamma _0\gt \sqrt {\psi _2}$ , we have the discriminant $\Delta _1\gt 0$ , and there exist two positive roots (i.e. $\xi _1$ and $\xi _2$ ) for the cubic equation (4.9). Therefore, the resonance condition (4.6) holds when the wavenumbers $\kappa _x$ are $\kappa _x=\sqrt {\xi _1}$ and $\kappa _x=\sqrt {\xi _2}$ . The function $\psi _2$ monotonically increases with respect to $\kappa _y$ when $\kappa _y$ is non-negative. Hence, when $\kappa _y$ equals zero, $\psi _2$ obtains its minimum value of $\psi _2=4$ . In conclusion, the resonance curve does not exist in the wavenumber space if the Froude number $\gamma _0\lt \min (\sqrt { \psi _2})=2$ , i.e. $U^2/\sqrt {g\sigma /\rho _w}\lt 2$ .

4.2. Characterisation of the resonance curve

As illustrated in figure 4, for gravity–capillary waves, the dimensionless resonance curve $\chi ^*=0$ in the $(\kappa _x,\kappa _y)$ space can be characterised by three representative points: the two intersections between the resonance curve and the $\kappa _x$ axis and the maximum point $(\kappa _{x,M},\kappa _{y,M})$ , where $\kappa _{y,M}$ is the largest value. To further illustrate the shape of the resonance curve, we estimate the value of $\kappa _{y,M}$

(4.16) \begin{align} \kappa _{y,M}=\max (\kappa _y:\chi ^*(\kappa _x,\kappa _y)=0). \end{align}

We aim to determine an upper bound for $\kappa _y$ that satisfies the resonance condition. Rearranging equation (4.8), we obtain

(4.17) \begin{align} &\big(1+\kappa _y^2\big)\big(1+3\kappa _y^2\big)\kappa _x^2+\kappa _y^2\big(1+\kappa _y^2\big)^2=\kappa _x^4\big(\gamma _0^2-3\kappa _y^2-2-\kappa _x^2\big). \end{align}

The right-hand side of (4.17) has an upper bound due to the inequality of the arithmetic and geometric means. Thus, we have

(4.18) \begin{align} {\kappa _y^6\leqslant \big(1+\kappa _y^2\big)\big(1+3\kappa _y^2\big)\kappa _x^2+\kappa _y^2\big(1+\kappa _y^2\big)^2 \leqslant \frac {4}{27}\big ({\gamma _0^2-3\kappa _y^2-2}\big )^3.} \end{align}

Therefore, for any $\kappa _y$ on the resonance curve $\chi ^*(\kappa _x,\kappa _y)=0$ , we obtain its upper bound as follows:

(4.19) \begin{align} {\kappa _y\leqslant \frac {\sqrt [3]{2}}{\sqrt {3+3\sqrt [3]{4}}}\sqrt {\gamma _0^2-2}.} \end{align}

Here, we use the positivity of $(1+\kappa _y^2)(1+3\kappa _y^2)\kappa _x^2$ to estimate $\kappa _y$ without explicit knowledge of $\kappa _x$ .

Figure 4. Illustration of the resonance curve. The red solid line represents the resonance curve $\chi ^*(\kappa _x,\kappa _y)=0$ based on (4.6). The orange dashed line represents the auxiliary function defined in (4.20). The dashed blue line denotes the theoretical upper bound, as defined in (4.19), for $\gamma _0=3.44$ .

We next determine the maximum value of $\kappa _y$ on the resonance curve $\chi ^*(\kappa _x,\kappa _y)=0$ . We observe that the derivative of $\kappa _y$ with respect to $\kappa _x$ reaches zero at $\kappa _y=\kappa _{y,M}$ . Taking the derivative of (4.8) with respect to $\kappa _x$ and evaluating it at $\kappa _y=\kappa _{y,M}$ , we obtain the following auxiliary function:

(4.20) \begin{align} 3\kappa _x^4 + 2\big(-\gamma _0^2 + 2 + 3\kappa _y^2\big)\kappa _x^2 + \big(1 + \kappa _y^2\big)\big(1 + 3\kappa _y^2\big) = 0, \end{align}

which is plotted in figure 4 for $\gamma _0=3.44$ when the convection velocity is $U=0.3\,\mathrm {m\,s^{-1}}$ . We can obtain the value of $\kappa _{y,M}$ by solving (4.8) and (4.20) simultaneously, which is the interception of the two curves plotted in figure 4. Figure 5 shows the variation in $\kappa _{y,M}$ and its upper bound, as defined in (4.19), for different $\gamma _0$ values. The results show that our estimation of the upper bound reasonably captures the variation in $\kappa _{y,M}$ .

The points at which the resonance curve intersects the $\kappa _x$ -axis can be directly calculated. At $\kappa _y=0$ , we have the following solutions to the resonance curve $\chi ^*(\kappa _x,\kappa _y)=0$ :

(4.21) \begin{align} &\kappa _{x,L} = \frac {1}{2}\gamma _0-\frac {1}{2}\sqrt {\gamma _0^2-4},\\[-10pt]\nonumber \end{align}

(4.22) \begin{align} &\kappa _{x,R} = \frac {1}{2}\gamma _0+\frac {1}{2}\sqrt {\gamma _0^2-4}. \end{align}

The resonance curve exists when $\gamma _0\gt 2$ and the distance between the two intersections, $|\kappa _{x,L}-\kappa _{x,R}|=\sqrt {\gamma _0^2-4}$ , increases as $\gamma _0$ increases. The area bounded by the resonance curve and the $x$ -axis scales as $\gamma _0^2$ , which also increases with increasing $\gamma _0$ .

In summary, we present a quantitative estimation for the shape of the resonance curve for gravity–capillary waves. To accurately compute this curve, one must solve an implicit cubic equation. Our approach simplifies the expression for estimating the resonancecurve height. It also obtains the intersections between the resonance curve and the $\kappa _x$ -axis. The results demonstrate that as the wind speed decreases (i.e. $\gamma _0$ decreases), the size of the domain bounded by the resonance curve and the $\kappa _x$ -axis decreases. When $\gamma _0$ is less than $2$ , the variables $\kappa _{x,L}$ and $\kappa _{x,R}$ , defined in (4.21) and (4.22), are no longer realvalued. This finding implies that when the wind speed falls below a certain threshold, the resonance curve does not exist. In § 4.3, we investigate the wave evolution under such weak wind conditions when the wind speed is too low to trigger the initial stage of wind-wave generation. The strong wind conditions are investigated in § 4.4.

Figure 5. Variation in the maximum value $\kappa _{y,M}$ as a function of $\gamma _0$ . The solid red line represents $\kappa _{y,M}$ solved according to (4.8) and (4.20), and the dashed blue line denotes its theoretical upper bound, as defined in (4.19).

4.3. Weak wind conditions

Following the previous discussion in § 4.2, the resonance mechanism does not occur for gravity–capillary waves under weak wind conditions. In such scenarios, the convection velocity of the air pressure fluctuations at the water’s surface is consistently lower than the phase velocity of the water waves for all wavenumbers. Therefore, the inequality $B\lt \Lambda$ in (3.2) holds for any wavenumber $(k_x,k_y)$ , and there is no singularity in the denominator of $F(k_x,k_y,t)$ , as shown in (2.11). Hence, $F(k_x,k_y,t)$ is bounded by a function that is independent of time $t$

(4.23) \begin{align} F(k_x,k_y,t)&=\frac {\Lambda }{(B^2-\Lambda ^2)^2} \left \{\frac {1}{2}B^2+\frac {3}{2}\Lambda ^2-\Lambda (\Lambda +B)\cos [(B-\Lambda )t]\right .\notag \\ & \quad \left .-\Lambda (\Lambda -B)\cos [(\Lambda +B)t]-\frac {1}{2}(B^2-\Lambda ^2)\cos (2\Lambda t)\right \}\notag \\ &\leqslant \frac {4\Lambda ^3}{(\Lambda ^2-B^2)^2}. \end{align}

Therefore, function $G(k_y,t)$ , i.e. the integral of $F(k_x,k_y,t)$ with respect to $k_x$ , satisfies

(4.24) \begin{align} G(k_y,t)\leqslant \int _0^\infty \frac {4\Lambda ^3}{(\Lambda ^2-B^2)^2} \mathrm {d}k_x. \end{align}

According to the above equation, $G(k_y,t)$ does not increase over time but instead reaches a finite value.

Figure 6 shows examples of the temporal evolution of $G(k_y,t)$ for various values of $U$ and $k_y$ . We set $g=9.8\,\mathrm {m\,s}^{-1}$ , $\rho _w=10^{3}\,\mathrm {kg\,m^{-3}}$ and $\sigma =7\times 10^{-2}\,\mathrm {kg\,s^{-2}}$ . The function $G(k_y,t)$ increases from zero over time and reaches a certain magnitude before beginning to oscillate. The amplitude of the oscillation in $G(k_y,t)$ decreases over time. At a fixed convection velocity of the pressure fluctuations $U$ , the value of the function $G(k_y,t)$ decreases for larger values of $k_y$ . Given a fixed wavenumber $k_y$ , the magnitude of $G(k_y,t)$ increases as the convection velocity $U$ increases. Thus, we can reasonably anticipate that higher wind speeds would induce larger surface elevation deformations. This expectation holds even when the wind speed is not sufficient to trigger wind-wave generation (the latter being a process in which the wave energy increases continuously over time).

Next, we aim to find an analytical expression that provides an upper bound for $G(k_y,t)$ for all values of time $t$ , which will give us an estimation of the wave energy when the wind speed is too low to trigger the initial stage of wind-wave generation.

Figure 6. Time evolution of $G(k_y,t)$ for gravity–capillary waves under weak wind conditions for various parameter values. In all the cases, we set $g=9.8\,\mathrm {m\,s^{-2}}$ , $\rho _w=10^3\,\mathrm {kg\,m^{-3}}$ and $\sigma =7\times 10^{-2}\,\mathrm {kg\,s^{-2}}$ . The blue, red and purple lines correspond to the results obtained for $U=0.05\,\mathrm {m\,s^{-1}}$ , $0.1\,\mathrm {m\,s^{-1}}$ and $0.15\,\mathrm {m\,s^{-1}}$ , respectively. The solid lines, dash-dotted lines and dotted lines represent results for $k_y=0\,\mathrm {m^{-1}}$ , $50\,\mathrm {m^{-1}}$ and $100\,\mathrm {m^{-1}}$ , respectively.

Under weak wind conditions, $\Lambda \gt B$ for any wavenumber $(k_x,k_y)$ . We perform inequality scaling for the denominator of (4.23) and obtain an upper-bound estimation of the integral of $F(k_x,k_y,t)$ with respect to $k_x$ . Our objective is to find a coefficient, denoted as $m\in (0,1)$ , that does not rely on $k_x$ and satisfies the condition

(4.25) \begin{align} m\Lambda ^2 \geqslant B^2, \end{align}

for all values of $k_x$ and $k_y$ . Therefore, we have $\Lambda ^2-B^2\gt (1-m)\Lambda ^2$ . The right-hand side of (4.23) is smaller than a power-law function of $\Lambda$ , whose integral can be estimated by a simple analytical solution. Based on the inequality of arithmetic and geometric means, we have the following estimation:

(4.26) \begin{align} \Lambda ^2(k_x,k_y)=g|\boldsymbol {k}|+\sigma \rho _w^{-1}|\boldsymbol {k}|^3\geqslant 2\sqrt {g\sigma \rho _w^{-1}}(k_x^2+k_y^2). \end{align}

Because $B=Uk_x$ , we combine (4.25) and (4.26) and obtain that (4.25) holds as long as $m$ satisfies

(4.27) \begin{align} m\geqslant \frac {U^2}{2\sqrt {g\sigma \rho _w^{-1}}}= {\frac {U^2}{c_{{min}}^2}}. \end{align}

The existence of $m$ is ensured by the weak wind condition, i.e. $U\lt c_{{min}}$ . Because we have $c_{{min}}=\sqrt {2\sqrt {g\sigma \rho _w^{-1}}}$ , the inequality $U\lt c_{{min}}$ is equivalent to $\ {U^2}/\big({2\sqrt {g\sigma \rho _w^{-1}}}\big)\lt 1$ . Therefore, we can find an $m$ within the interval $(0,1)$ such that $m$ satisfies (4.27). Consequently, we can rescale the denominator of (4.24) to eliminate the function $B$ and obtain the following result:

(4.28) \begin{align} G(k_y,t) &\leqslant \frac {4}{(1-m)^2}\int _0^\infty \frac {1}{\Lambda }\mathrm {d}k_x. \end{align}

However, integrating $\Lambda ^{-1}$ with respect to $k_x$ still does not lead to an analytical solution. To obtain an analytical upper-bound estimation, we partition the integral domain $(0,\infty )$ into two parts: $(0,k_*)$ and $(k_*,\infty )$ , where $k_*$ is an arbitrary positive number. To ensure that the integral remains well defined, we include only the gravity-related term in the partition $(0,k_*)$ and only the surface tension-related term in the partition $(k_*,\infty )$ . This choice preserves sufficient physical information when determining the upper-bound estimation. Hence, we obtain the following estimation for function $G(k_y,t)$ :

(4.29) \begin{align} G(k_y,t)&\leqslant \frac {4}{(1-m)^2}\bigg[\int _0^{k_*} \frac {1}{\sqrt {g\sqrt {k_x^2+k_y^2}}}\mathrm {d}k_x+ \int _{k_*}^\infty \frac {1}{\sqrt {\sigma \rho _w^{-1}\left (\sqrt {k_x^2+k_y^2}\right )^3}}\mathrm {d}k_x\bigg]\notag \\ &\leqslant \frac {4}{(1-m)^2}\left [\int _0^{k_*} \frac {2^{\frac {1}{4}}}{\sqrt {g}(k_x+k_y)^\frac {1}{2}}\mathrm {d}k_x + \int _{k_*}^\infty \frac {2^{\frac {3}{4}}}{\sqrt {\sigma \rho _w^{-1}}(k_x+k_y)^\frac {3}{2}}\mathrm {d}k_x\right ]\notag \\ &\leqslant \frac {2^{\frac {13}{4}}}{(1-m)^2}\left [\frac {\sqrt {k_y+k_*}-\sqrt {k_y}}{\sqrt {g}}+\frac {\sqrt {2}}{\sqrt {\sigma \rho _w^{-1}}\sqrt {k_y+k_*}}\right ]. \end{align}

We use the inequality $k_x^2+k_y^2\geqslant (k_x+k_y)^2/2$ to obtain (4.29), which rescales the integral of both the gravity- and surface tension-related terms to analytic solutions. Note that (4.29) holds for any value of $k_*$ . We choose a specific value of $k_*$ that ensures that the right-hand side of (4.29) reaches a minimum when $k_y$ approaches zero. By calculations, we obtain that $k_*=\sqrt {2}k_0$ .

Figure 7. Evolution of $G(k_y,t)$ for gravity–capillary waves under weak wind conditions for different parameters. The blue solid lines denote the integral representation, i.e. the right-hand side of (4.24), and red dashed lines show the analytical representation of the upper bound, i.e. the right-hand side of (4.30). Panel $(a)$ shows a plot of the representation functions versus the variable $U$ , with $k_y=0$ . Panel $(b)$ shows a plot of the representation functions versus the variable $k_y$ , with $U=0.1\,\mathrm {m\,s^{-1}}$ . In all cases, we choose $g=9.8\,\mathrm {m\,s^{-2}}$ , $\rho _w=10^3\,\mathrm {kg\,m^{-3}}$ and $\sigma =7\times 10^{-2}\,\mathrm {kg\,s^{-2}}$ .

Using the results of our previous analysis on the admissible range of $m$ (see (4.27)), we choose $m=U^2/(2\sqrt {g\sigma \rho _w^{-1}})=\gamma _0/2$ , and obtain

(4.30) \begin{align} G(k_y,t)\leqslant \frac {2^{\frac {21}{4}}}{\left (2-\gamma _0\right )^2}\left [\frac {\sqrt {k_y+\sqrt {2}k_0}-\sqrt {k_y}}{\sqrt {g}}+\frac {\sqrt {2}}{\sqrt {\sigma \rho _w^{-1}}\sqrt {k_y+\sqrt {2}k_0}}\right ]. \end{align}

This estimation has useful physical implications. Equation (4.30) shows that the wave energy is bounded by a function that does not depend on time $t$ when the convection velocity of the air pressure fluctuations is below the critical value, which is the minimum phase velocity of water waves. Additionally, the saturated wave energy increases as the convection velocity increases. This relation between $G$ and $k_y$ can be supported by the upper-bound estimation of $G$ , derived in (4.30). Both terms in parentheses on the right-hand side of (4.30) decrease as $k_y$ increases. Therefore, the upper bound of $G$ decreases as $k_y$ increases. The original analysis by Phillips (Reference Phillips1957) considered only the scenario in which the convection velocity is zero. Our approach here extends the wave energy analysis to weak wind conditions (and the strong wind conditions, which are analysed in § 4.4), enabling a quantitative estimation of the dependency of the saturated wave energy on the convection velocity of the air pressure fluctuations at the water’s surface. Perrard et al. (Reference Perrard, Lozano-Durán, Rabaud, Benzaquen and Moisy2019) investigated quasisteady surface elevation motions influenced by turbulent airflow based on the balance between the energy input from the airflow and the viscous dissipation in the water. An experimental study by Paquier, Moisy & Rabaud (Reference Paquier, Moisy and Rabaud2015) reported that surface deformations can be saturated at low wind speeds. Here, we present a physical mechanism for the generation of quasisteady surface elevation motions, even without the effect of viscous dissipation, when the airflow speed is inadequate to induce the resonance mechanism.

In figure 7, we present numerical results to support the above estimate of $G(k_y,t)$ . The blue solid lines in figure 7 represent the integral representation of the right-hand side of (4.24), while the red dashed lines depict the analytical representation, that is, the right-hand side of (4.30). The results show that the estimation using the analytical representation defined in (4.24) reasonably captures the characteristics of the upper bound of the function $G(k_y,t)$ under weak wind conditions.

Under weak wind conditions, the wave energy eventually reaches a saturation value because the resonance condition for linear growth does not exist. Numerical examples confirm that the saturation energy increases with the wind speed, and our derived upper bound for $G$ under weak wind conditions approaches infinity when $U$ approaches $c_{{min}}$ . This result indicates that, under weak wind conditions, as the wind speed approaches the minimum phase velocity of gravity–capillary water waves, $G$ can be extremely large. This result suggests a continuous transition from weak to strong wind conditions. Furthermore, it is also important to emphasise that our present analysis applies only to the initial stage of wind-wave generation. The governing equations for wave motion in this stage assume Taylor’s frozen turbulence hypothesis. Over longer time scales, time decorrelations of pressure fluctuations occur, and the wave evolution depends on the space-time correlation of air turbulence (Li & Shen Reference Li and Shen2022 b).

4.4. Strong wind conditions

In this subsection, we analyse wave evolution under strong wind conditions, namely, when the convection velocity of the air pressure fluctuations exceeds the minimum phase velocity of the water waves. As stated in the previous discussions, such circumstances are met when the non-dimensional Froude number satisfies $\gamma _0 \gt 2$ ; hence, resonance occurs for wavenumbers within the range of $k_y \leqslant k_{y,M}$ . Here, $k_{y,M}$ represents the maximum wavenumber for which resonance occurs and is expressed as $k_{y,M} = k_0\kappa _{y,M}$ , where $\kappa _{y,M}$ is defined in (4.16). We employ a technique similar to the method proposed in § 3 to analyse the wave energy evolution for gravity–capillary waves next. We decompose the evolution of $G(k_y,t)$ into a linear growth term $M(k_y,t)$ and a residual term $N(k_y,t)$ . Our numerical results show that $M(k_y,t)$ accounts for the leading-order growth of $G(k_y,t)$ in most cases.

4.4.1. Analysis in the complex plane

Similar to the analysis of gravity waves, for gravity–capillary waves we define the analytical continuation of $F$ as the following complex-valued function $F_2(z,k_y,t)$ :

(4.31) \begin{align} F_2(z,k_y,t)=&\frac {\Lambda }{(B^2-\Lambda ^2)^2} \left \{\frac {1}{2}B^2+\frac {3}{2}\Lambda ^2-\Lambda (\Lambda +B)\exp [\mathrm {i}(\Lambda -B)t]\right .\notag \\ &\left .-\Lambda (\Lambda -B)\exp [\mathrm {i}(\Lambda +B)t]-\frac {1}{2}\big(B^2-\Lambda ^2\big)\exp (2\mathrm {i}\Lambda t)\right \}. \end{align}

Here, the complex-valued functions $B$ and $\Lambda$ are defined as

(4.32) \begin{align} B(z)&=Uz, \end{align}

(4.33) \begin{align} \Lambda (z)&=\sqrt {g\sqrt {z^2+k_y^2}+\sigma \rho _w^{-1}\left (\sqrt {z^2+k_y^2}\right )^3}. \end{align}

The term $\exp [\mathrm {i}(\Lambda -B)t]$ in (4.31) ensures that $F_2(z,k_y,t)$ approaches zero for large values of $|z|$ in the first quadrant of the complex plane for gravity–capillary waves.

Figure 8. Illustration of the closed loop for integration, as outlined in (4.34).

When the spanwise wavenumber $k_y$ is less than the maximum wavenumber for which the resonance occurs, i.e. $k_{y,M}$ , two points $(k_{x1},k_y)$ and $(k_{x2},k_y)$ exist that satisfy the resonance condition. We consider a closed loop $\Gamma$ in the complex plane, as depicted in figure 8, which is defined as $\Gamma =[0, z_1-r]+C_{r1}+[z_1+r, z_2-r]+C_{r2}+$ $[z_2+r, R]+C_R+L_R$ . Here, $z_1$ and $z_2$ are the two roots of $B(z)-\Lambda (z)=0$ ; $C_{r1}( {\theta })=z_1+re^{\mathrm {i}(\pi -\theta )}$ , with $0\leqslant \theta \leqslant \pi$ , is a clockwise circular path around $z_1$ with a radius of $r$ ; $C_{r2}( {\theta })=z_2+re^{\mathrm {i}(\pi -\theta )}$ , with $0\leqslant \theta \leqslant \pi$ , is a clockwise circular path around $z_2$ with a radius of $r$ ; $C_R(\theta )=Re^{\mathrm {i}\theta }$ , with $0\leqslant \theta \leqslant \pi /4$ , is an anticlockwise circular path around 0 with a radius of $R$ ; and $L_R(\theta )=Re^{\mathrm {i}\pi /4}(1-\theta )$ , with $0\leqslant \theta \leqslant 1$ , is a straight path pointing towards the origin. Because the function $F_2$ has no singularities inside the closed loop $\Gamma$ , the integral of the function $F_2(z,k_y,t)$ with respect to $z$ along the contour $\Gamma$ is equal to zero, which can be expressed as

(4.34) \begin{align} \left (\int _0^{z_1-r}+\int _{C_{r1}}+\int _{z_1+r}^{z_2-r}+\int _{C_{r2}} +\int _{z_2+r}^R+\int _{C_R}+\int _{L_R}\right )F_2(z,k_y,t)\mathrm {d}z=0. \end{align}

Taking the limits of $R\rightarrow \infty$ and $r\rightarrow 0$ , the integration of $F_2$ along the real positive axis can be calculated as

(4.35) \begin{align} \int _0^{\infty }F(k_x,k_y,t) \mathrm {d}k_x=&\lim _{R\rightarrow \infty }\int _{0}^{R}F(k_x,k_y,t) \mathrm {d}k_x\notag \\ =&\lim _{R\rightarrow \infty }\lim _{r\rightarrow 0}\mathrm {Re}\left [\left (\int _0^{z_1-r}+\int _{z_1+r}^{z_2-r}+\int _{z_2+r}^{R}\right )F_2(z,k_y,t)\mathrm {d}z\right ]\notag \\ =&\lim _{R\rightarrow \infty }\lim _{r\rightarrow 0}\mathrm {Re}\left [\left (-\int _{C_{r1}}-\int _{C_{r2}}-\int _{C_R}-\int _{L_r}\right )F_2(z,k_y,t)\mathrm {d}z\right ]\notag \\ =&-\lim _{R\rightarrow \infty }\mathrm {Re}\left [\int _{C_R}F_2(z) \mathrm {d}z\right ] -\lim _{r\rightarrow 0}\mathrm {Re}\left [\int _{C_{r1}}F_2(z) \mathrm {d}z\right ]\notag \\ &-\lim _{r\rightarrow 0}\mathrm {Re}\left [\int _{C_{r2}}F_2(z) \mathrm {d}z\right ]-\lim _{R\rightarrow \infty }\mathrm {Re}\left [\int _{C_{L_R}}F_2(z) \mathrm {d}z\right ]. \end{align}

Similar to the above analysis for gravity waves, the integration of $F(k_x,k_y,t)$ with respect to $k_x$ along the positive $x$ -axis is equivalent to the integration of the complex function $F_2(z,k_y,t)$ with respect to the complex value $z$ along the paths $C_{r1}$ , $C_{r2}$ , $C_R$ and $L_R$ . Because the magnitudes of $B$ and $\Lambda$ behave asymptotically as $|B|\sim |z|$ and $|\Lambda |\sim |z|^{3/2}$ when $|z|\rightarrow \infty$ , the function $F_2(z,k_y,t)$ , defined in (4.31), behaves as $|F_2(z,k_y,t)|\sim |z|^{-3/2}$ for large values of $|z|$ when $z$ is located in the first quadrant of the complex $z$ -plane. Therefore, the integration of $F_2$ along the path $C_R$ becomes zero as $R$ approaches infinity. Hence, we next focus on the integration of $F_2(z,k_y,t)$ along the paths $C_{r1}$ , $C_{r2}$ and $C_{L_R}$ .

4.4.2. Calculation of the integral of $F_2$ along the paths $C_{r1}$ and $C_{r2}$

The integral of $F_2$ along the semicircular paths $C_{r1}$ and $C_{r2}$ can be calculated using the residue theorem as follows:

(4.36) \begin{align} \lim _{r\rightarrow 0}\int _{C_{r1}}F_2(z) \mathrm {d}z&=-\pi \mathrm{iRes}(F_2,z_{1}),\notag \\ \lim _{r\rightarrow 0}\int _{C_{r2}}F_2(z) \mathrm {d}z&=-\pi \mathrm{iRes}(F_2,z_{2}). \end{align}

The varibales $z_1$ and $z_2$ are two real-valued solutions for $z$ to the following resonance condition equation:

(4.37) \begin{align} Uz=\sqrt {g\sqrt {z^2+k_y^2}+\sigma \rho _w^{-1}\left (\sqrt {z^2+k_y^2}\,\right )^3}. \end{align}

Both singularity points (i.e. $z=z_1$ and $z=z_2$ ) are first-order poles of the complex-valued function $F_2(z,k_y,t)$ . The residues at $z=z_i$ for $i=1,2$ are defined as

(4.38) \begin{align} \mathrm {Res}(F_2,z_i)=\lim _{z\rightarrow z_i} F_2(z)(z-z_i). \end{align}

To compute these residues, we follow the same perturbation method introduced in § 3.2.3. We begin by calculating the Taylor expansions of $B-\Lambda$ around $z-z_1$ and $z-z_2$ , respectively. These expansions are then substituted into (4.38), yielding

(4.39) \begin{align} \mathrm {Res}( {F_2},z_1) &=-\frac { \mathrm {i}U^2 z_1\sqrt {z_1^2+k_y^2}t} {g+3\sigma \rho _w^{-1}(z_1^2+k_y^2)-2U^2\sqrt {z_1^2+k_y^2}},\notag \\ \mathrm {Res}\left( {F_2},z_2\right) &=-\frac { \mathrm {i}U^2 z_2\sqrt {z_2^2+k_y^2}t} {g+3\sigma \rho _w^{-1}(z_2^2+k_y^2)-2U^2\sqrt {z_2^2+k_y^2}}. \end{align}

Hence, the integral of $F_2$ along two paths $C_{r1}$ and $C_{r2}$ is

(4.40) \begin{align} \int _{C_{r1}+C_{r2}}F_2(z,k_y,t)\mathrm{d}z=&-\frac {\pi U^2 z_1\sqrt {z_1^2+k_y^2}t} {g+3\sigma \rho _w^{-1}(z_1^2+k_y^2)-2U^2\sqrt {z_1^2+k_y^2}}\notag \\ &-\frac {\pi U^2 z_2\sqrt {z_2^2+k_y^2}t} {g+3\sigma \rho _w^{-1}(z_2^2+k_y^2)-2U^2\sqrt {z_2^2+k_y^2}}. \end{align}

Notably, this result indicates that the integral of $F_2$ along the paths $C_{r1}$ and $C_{r2}$ is a linear function of time $t$ . At $k_y=k_{y,M}$ , the two roots $z_1$ and $z_2$ coincide and are denoted as $z_M$ . In other words, the integral in (4.40) can then be evaluated with $z_1=z_2=z_M$ . Next, we analyse the integral of $F_2$ along the path $L_R$ .

4.4.3. Calculation of the integral of $F_2$ along path $L_R$

Next, we consider the integral of $F_2$ along path $L_R$ , where the variable $z$ can be parameterised as a function of the real variable $r$ , such that $z(r) = r\exp (\mathrm {i}\pi /4)$ with $0\leqslant r\leqslant R$ . Therefore, we have

(4.41) \begin{align} \int _{L_R} F_2(z)\mathrm {d}z&=\exp \left (\frac {\mathrm {i}\pi }{4}\right )\int _R^0 F_2(z(r))\mathrm {d}r\notag \\ &=-\frac {1+\mathrm {i}}{\sqrt {2}}\int _0^R\frac {K(r,t)}{J(r)}\mathrm {d}r. \end{align}

In (4.41), the numerator $K(r,t)$ and denominator $J(r)$ are

(4.42) \begin{align} K(r,t) & =\Lambda _0(r)\left \{\frac {1}{2}B_0(r)^2+\frac {3}{2}\Lambda _0(r)^2-\Lambda _0(r)[\Lambda _0(r)+B_0(r)]\exp [\mathrm {i}(\Lambda _0(r)-B_0(r))t]\right .\notag \\ & \quad -\Lambda _0(r)[\Lambda _0(r)-B_0(r)]\exp [\mathrm {i}(B_0(r)+\Lambda _0(r))t]\notag \\ &\quad \left .-\frac {1}{2}\big[B_0(r)^2-\Lambda _0(r)^2\big]\exp [2\mathrm {i}\Lambda _0(r) t]\right \}, \end{align}

(4.43) \begin{align} J(r)& =\big[B_0(r)^2-\Lambda _0(r)^2\big]^2, \end{align}

respectively, and $B_0(r)$ and $\Lambda _0(r)$ are defined as

(4.44) \begin{align} B_0(r)&=\frac {1+\mathrm {i}}{\sqrt {2}}Ur,\notag \\ \Lambda _0(r)&=\sqrt {g\sqrt {\mathrm {i}r^2+k_y^2}+\sigma \rho _w^{-1}\left (\sqrt {\mathrm {i}r^2+k_y^2}\right )^3}. \end{align}

We note that the definitions of $B_0(r)$ and $\Lambda _0(r)$ are based on the angle between $L_R$ and the $x$ -axis being $\pi /4$ . However, this choice of angle does not influence the real part of the integral of $F_2$ along $L_R$ because the contribution from the integral along $C_R$ vanishes as $R\rightarrow \infty$ . This phenomenon has been verified through numerical integrations for different angles.

4.4.4. Properties of the function G

Based on (4.35), (4.40) and (4.41) and taking the limit $R\rightarrow \infty$ , we can obtain the following expression for $G(k_y,t)$ , defined in (2.13):

(4.45) \begin{align} G(k_y,t)=M(k_y,t)+N(k_y,t). \end{align}

Here, the functions $M(k_y,t)$ and $N(k_y,t)$ are

(4.46) \begin{align} M(k_y,t)&=\frac {\pi z_1\sqrt {z_1^2+k_y^2}U^2 t} {g+3\sigma \rho _w^{-1}(z_1^2+k_y^2)-2U^2\sqrt {z_1^2+k_y^2}}\nonumber\\& \quad +\frac {\pi z_2\sqrt {z_2^2+k_y^2}U^2 t} {g+3\sigma \rho _w^{-1}(z_2^2+k_y^2)-2U^2\sqrt {z_2^2+k_y^2}}, \end{align}

(4.47) \begin{align} N(k_y,t)=&\mathrm {Re}\left [\frac {1+\mathrm {i}}{\sqrt {2}}\int _0^\infty \frac {K(r,t)}{J(r)}\mathrm {d}r\right ], \end{align}

where the expressions of the numerator $K(r,t)$ and the denominator $J(r)$ on the right-hand side of (4.47) are given in (4.42) and (4.43), respectively. Due to its complexity, we do not explicitly write the complete expression of $N(k_y,t)$ here. The function $M(k_y,t)$ is the leading-order term of $G(k_y,t)$ in most cases, and $N(k_y,t)$ is the residual. The function $N(k_y,t)$ can be computed through numerical integration using relationships in (4.42), (4.43) and (4.47). Unlike the gravity wave case, in (4.46), $M(k_y,t)$ does not exhibit an explicit scaling behaviour with the convection velocity $U$ . Assuming that the two roots are ordered such that $z_1\lt z_2$ , the smaller root $z_1$ (corresponding to the smaller wavenumber) can be treated as the ‘gravity approximation’ of the solution, while the larger root $z_2$ (larger wavenumber) represents the ‘capillary approximation’ of the solution. Preliminary numerical test indicates that the growth rate of $M(k_y,t)$ associated with $z_1$ (the first term on the right-hand side of (4.46)) is negative and scales as $U^n$ with $n\in [-2,-1]$ . The scaling behaviour is similar to the gravity wave case (see (3.38)). On the other hand, the growth rate of $M(k_y,t)$ related to $z_2$ is postive and scales as $U^2$ when $|k_y|\ll z_2$ , reflecting the capillary effect.

In the gravity wave case, we demonstrated that $N(k_y,t)$ is an oscillating function with a bounded amplitude. However, in the gravity–capillary wave case, due to the high-order nonlinearity term in the variable $\Lambda$ , the residual term $N(k_y,t)$ both oscillates and increases linearly over time, but its magnitude is lower than that of the leading-order term $M(k_y,t)$ . A special case in which $N(k_y,t)$ does not increase over time occurs when $g=0$ and $k_y=0$ . Numerical examples are provided below to support the above theoretical findings.

Figure 9 shows numerical examples of calculating the function $G(k_y,t)$ with $M(k_y,t)$ and $N(k_y,t)$ for a gravity–capillary wave with different convection velocity $U$ (i.e. $U=0.3\,\mathrm {m\,s^{-1}}$ and $U=1\,\mathrm {m\,s^{-1}}$ ) and a capillary wave (a limiting case when $g=0\,\mathrm {m\,s^{-2}}$ ). The relationship among $G(k_y,t)$ , $M(k_y,t)$ and $N(k_y,t)$ and their temporal evolution properties are visualised in panels (a), (b) and (c) of figure 9. The results confirm that the function $G(k_y, t)$ predominantly increases linearly over time and that the function $M(k_y, t)$ captures the leading-order behaviour of $G(k_y, t)$ . In panels (a) and (b), the convection velocities are $U=0.3\,\mathrm {m\,s^{-1}}$ and $U=1\,\mathrm {m\,s^{-1}}$ , respectively, and in both cases, the strong wind condition is satisfied: $U\gt c_{{min}}=0.23\,\mathrm {m\,s^{-1}}$ . For the larger convection velocity case (i.e. panel b), the magnitude of $M(k_y,t)$ is closer to that of $G(k_y,t)$ . The time evolution of the function $N(k_y,t)$ , which is the oscillatory component of $G(k_y,t)$ , is shown in panels (d) and (e) of figure 9. In general, $N(k_y,t)$ also increases over time and oscillates, and the magnitude of $N(k_y,t)$ is smaller than that of $M(k_y,t)$ , confirming that $M(k_y,t)$ represents the leading-order behaviour of the temporal evolution of $G(k_y,t)$ . The special case of $g=0$ and $k_y=0$ is shown in the right column of figure 9. In this case, the residual term $N(k_y,t)$ does not increase and instead only oscillates over time, and its amplitude decreases gradually. Panels (g) and (h) of figure 9 show the numerical error for the integration of $G(k_y,t)$ , calculated as $G(k_y,t)-M(k_y,t)-N(k_y,t)$ . The relative error between $G(k_y,t)$ and $M(k_y,t)+N(k_y,t)$ normalised by $G(k_y,t)$ is less than $4\times 10^{-5}$ , which confirms the accuracy of the numerical integrations.

Figure 9. Numerical examples of the calculation of $G(k_y,t)$ , with the results for $k_y=0\,\mathrm {m^{-1}}$ , $U=0.3\,\mathrm {m\,s^{-1}}$ and $g=9.8\,\mathrm {m\,s^{-2}}$ shown in the left column (panels a, d and g), the results for $k_y=10\,\mathrm {m^{-1}}$ , $U=1\,\mathrm {m\,s^{-1}}$ , $U = 0.3\, \mathrm {m s^{-1}}$ and $g=9.8\,\mathrm {m\,s^{-2}}$ shown in the middle column (panels b, e and h) and the results for $k_y=0\,\mathrm {m^{-1}}$ and $g=0\,\mathrm {m\,s^{-2}}$ shown in the right column (panels c, f and i). The surface tension $\sigma$ and water density $\rho _w$ are set to $\sigma =7\times 10 ^{-2}\,\mathrm {kg\,s^{-2}}$ and $\rho _w=10^3\,\mathrm {kg\,m^{-3}}$ , respectively. Panels (a), (b) and (c) show $G(k_y,t)$ , with the blue lines indicating the numerical integration results, the red lines indicating the linear growth component $M(k_y,t)$ and the black dots indicating the sum of the linear growth component $M(k_y,t)$ and the oscillatory component $N(k_y,t)$ . In panels (d), (e) and (f) the red lines represent the oscillatory component $N(k_y,t)$ , and the blue lines represent $G(k_y,t)-M(k_y,t)$ . Panels (g), (h) and (i) show the numerical errors in the calculations, with the blue lines representing $G(k_y,t)-M(k_y,t)-N(k_y,t)$ .

To better understand the relative contributions of $M(k_y,t)$ and $N(k_y,t)$ to $G(k_y,t)$ , we calculate the ratio of $M(k_y,t)/G(k_y,t)$ at $t=1\,\mathrm {s}$ for different $k_y$ and $U$ . Figure 10 illustrates two numerical examples of $M(k_y,t)/G(k_y,t)|_{t=1\,\mathrm {s}}$ varying with the spanwise wavenumber $k_y$ . In both cases, the magnitude of $M(k_y,t)$ is more than one half of the magnitude of $G(k_y,t)$ for $k_y$ in the range of $(0,0.8k_{y,M})$ . For the higher wind velocity case (figure 10 b), $M(k_y,t)/G(k_y,t)|_{t=1\,\mathrm {s}}$ is above $0.998$ for a wide range of $k_y$ . Figure 11 plots the ratio of $M(k_y,t)/G(k_y,t)$ at $t=1\,\mathrm {s}$ , averaged over $k_y\in [0,100]\,\mathrm {m^{-1}}$ , and illustrates that this value approaches one as the convection velocity $U$ increases. Therefore, $M(k_y,t)$ is a reasonable representation of the leading-order linear behaviour of $G(k_y,t)$ in most cases.

Figure 10. Two numerical examples of the ratio $M(k_y,t)/G(k_y,t)$ at $t=1\,\mathrm {s}$ with respect to wavenumber $k_y$ . In each panel, the vertical dashed line indicates the location where $k_y=k_{y,M}$ , and the horizonal dashed line represents the value 1.

Figure 11. Ratio of $M(k_y,t)/G(k_y,t)$ at $t=1\,\mathrm {s}$ , averaged over $k_y$ , as a function of the convection velocity $U$ . The vertical dashed line indicates the location where $U=c_{{min}}$ , and the horizonal dashed line represents the value 1.

4.4.5. Overall increase in the surface elevation variance $\langle \eta ^2\rangle$ over time

When the spanwise wavenumber $k_y$ satisfies $k_y\lt k_{y,M}$ , the resonance condition can be met at two points of $k_x$ for a fixed value of $k_y$ . However, when $k_y$ is larger than the critical value $k_{y,M}$ , no resonance occurs. Hence, we partition the first quadrant of the wavenumber space into two parts: $\Omega _1=\{(k_x,k_y)\in [0,\infty )\times {[}0,k_{y,M})\}$ and $\Omega _2=\{(k_x,k_y)\in [0,\infty )\times (k_{y,M},\infty )\}$ . Based on the above analyses, the wave energy components within region $\Omega _1$ increase over time, whereas those within $\Omega _2$ remain finite. In $\Omega _2$ , it always holds that $\Lambda (\boldsymbol {k})\gt \boldsymbol {k} \cdot \boldsymbol {U}$ , meaning that the resonance condition is not satisfied. The wave dynamics in $\Omega _2$ can be analysed using a similar approach based on the weak wind condition discussed in § 4.3. Hence, the function $E(k_y,t)$ can be expressed as

(4.48) \begin{align} E(k_y,t)&\leqslant 2\Vert \Phi _p (k_x^2+k_y^2)\Lambda ^{-3} \Vert _{L^\infty _{x}}\left [M(k_y,t)+N(k_y,t)\right ], &(k_x,k_y)\in \Omega _1, \end{align}

(4.49) \begin{align} E(k_y,t)&\leqslant 2\Vert \Phi _p (k_x^2+k_y^2)\Lambda ^{-3} \Vert _{L^\infty _{x}}C_1, &(k_x,k_y)\in \Omega _2, \end{align}

where $C_1$ is a constant that does not depend on time $t$ .

Hence, the surface elevation variance $\langle \eta ^2\rangle$ can be calculated from (2.8) as follows:

(4.50) \begin{align} \langle \eta ^2\rangle & \leqslant 2\int _{-k_{y,M}}^{k_{y,M}} \Vert \Phi _p\big(k_x^2+k_y^2\big)\Lambda ^{-3}\Vert _{L_x^\infty } (M(k_y,t)+N(k_y,t) ) \mathrm {d}k_y\notag \\ & \quad +2\bigg(\int _{-\infty }^{-k_{y,M}}+\int _{k_{y,M}}^{\infty } \bigg) \Vert \Phi _p\big(k_x^2+k_y^2\big)\Lambda ^{-3}\Vert _{L_x^\infty } {C_1} \mathrm {d}k_y\notag \\& \quad \leqslant 2\Vert \Phi _p\big(k_x^2+k_y^2\big)\Lambda ^{-3}\unicode {x1D7D9}_{|k_y|\lt k_{y,M}}M(k_y,t) \Vert _{L_x^\infty L_y^1}\notag \\ & \quad + 2\Vert \Phi _p\big(k_x^2+k_y^2\big)\Lambda ^{-3}\unicode {x1D7D9}_{|k_y|\lt k_{y,M}}N(k_y,t) \Vert _{L_x^\infty L_y^1}+C_2, \end{align}

where $\unicode {x1D7D9}_{|k_y|\lt k_{y,M}}$ is the indicator function and $C_2$ is independent of time $t$ . The first term on the right-hand side of (4.50) represents the leading-order temporal increase in the wave energy, with the expression $M(k_y,t)$ explicitly given (see (4.46)), and this term demonstrates the linear increases in the surface elevation variance over time. The residue term related to $N(k_y,t)$ represents the oscillatory behaviour of the wave energy evolution as well as the subdominant linear growth behaviour. The last term on the right-hand side of (4.50) represents to the contribution of the wave energy components that have wavenumbers for which the resonance does not occur.

In summary of § 4, for gravity–capillary waves, we found that resonance occurs when the convection velocity of the air pressure fluctuations at the water’s surface is greater than the critical phase velocity of the water waves, which is equivalent to the condition for the non-dimensional Froude number $\gamma _0\gt 2$ . When wind speeds fall below this threshold, the wind-generated surface elevation does not continually increase over time and instead saturates within a certain period. When the resonance condition is satisfied (i.e. $\gamma _0\gt 2$ ), our analyses demonstrate that the overall wave energy increases over time, and its leading-order component exhibits linear growth behaviour.