3.1. Overview of wind turbulence field evolution

In this section, we first give an overview of the wind field evolution with a focus on the turbulent motion ($u', v', w'$) defined in (2.4). For the length consideration of the paper, we present results in the growing wind of case CU8 when the flow dynamics are similar among the different cases in table 1. Figure 2 shows $u'$ and $w'$ at the start and end of the growing wind gust. The magnitude of $u'$ is generally stronger than $w'$ at both $t=0$ and $t=T_g$, as shown by their contours on the $x$–$z$ and $y$–$z$ planes. The streaks near the wave surface, denoted by the isosurfaces of $u'$, are enhanced significantly after the impact of the wind gust. On the other hand, the isosurfaces of $w'$ at $t=0$ and $t=T_g$ are qualitatively similar. Explanations for these differences are provided from the analyses in § 3.2.

Figure 2. Turbulent motions at (a,b) $t=0$ and (c,d) $t=T_g$ for the growing wind in case CU8. Plotted are the contours of $u'$ and $w'$ on the $x$–$z$ and $y$–$z$ planes, and the isosurfaces corresponding to $u'^{+0}=\pm 7$ in (a,b) and $w'^{+0}=\pm 3$ in (c,d), where the superscript ‘$+0$’ denotes the normalization based on the wall unit at the initial state (see table 1).

The mean wind profiles before, during and shortly after the gust event are plotted in figure 3. During the gust event $0< t/T_g<1$, the mean flow away from the wave surface ($\zeta >0.01\lambda$) experiences a similar change and the velocity profiles are shifted by a roughly constant value. The dominant mechanism can be understood through the averaged streamwise momentum equation, which reduces to $\partial U/\partial t=B(t)$ to the leading order. Since $B(t)$ is nearly a constant, the time variation in the mean profile is linear and uniform. While a logarithmic region remains in the mean profile at the end of the gust $t=T_g$, the wind turbulence is in a non-equilibrium state. Therefore, the slope of the mean profile is no longer determined only by the friction velocity, and the mean profile cannot be used to calculate the wave growth rate in the same way as when the wind turbulence is steady (Miles Reference Miles1957; Belcher & Hunt Reference Belcher and Hunt1993) or quasi-stationary (Miles & Ierley Reference Miles and Ierley1998). In the region near the wave surface $\zeta <0.01\lambda$, the velocity gradient (i.e. the mean shear) changes significantly during the gust because of the boundary condition constraint imposed by the wave surface (see Appendix B). After the gust event, the wind field continues to evolve (comparing the profiles at $t=T_g$ and $t=1.5T_g$ in figure 3) and a new equilibrium state can be established after a sufficiently long time. Compared with the wind gust event, the temporal variations in the wind field after ${O}(2T_g\unicode{x2013}4T_g)$ are less drastic in the return to equilibrium and therefore excluded from the following analysis.

Figure 3. Evolution of the mean wind velocity profiles in (a) the growing wind and (b) the decaying wind, in case CU8. The profiles are normalized in the initial wall unit and plotted during $0\leq t/T_g\leq 1.5$. The black dashed line denotes $U=c$.

The response of turbulence to the wind gust is in sharp contrast to that of the mean velocity profile. Figure 4 shows the profiles of Reynolds stress components at several time instants during and after the wind gust event in case CU8. Here, $\langle v'v'\rangle$ is not plotted because its evolution is qualitatively similar to $\langle w'w'\rangle$. When the growing wind gust occurs ($0< t< T_g$), we observe a rapid increase in the streamwise velocity fluctuations $\langle u'u'\rangle$ (figure 4a), while $\langle w'w'\rangle$ (figure 4b) has relatively stable values, and the evolution speed of $\langle u'w'\rangle$ (figure 4c) is between those of $\langle u'u'\rangle$ and $\langle w'w'\rangle$. At $t=1.5T_g$, the magnitudes of all these Reynolds stress components become significantly larger as the turbulence field continues to evolve. In the decaying wind (figure 4d–f), the responses of the turbulent stress first appear near the wave surface and then extend to the outer region. These results are similar to the pipe flow experiment of accelerating and decelerating turbulence (He & Jackson Reference He and Jackson2000). The behaviours of the turbulent stresses are related closely to the gust-induced change in the mean shear. As shown above, the mean velocity profile exhibits a linear and uniform time variation such that the mean shear remains constant in the outer region and changes drastically only near the wave surface.

Figure 4. Vertical distribution of the Reynolds stresses for (a–c) the growing wind and (d–f) the decaying wind, in case CU8.

The evolution of the Reynolds stresses is also demonstrated by tracking their peak values in time, shown in figure 5. Here, we also plot the time variation of the mean flow energy $U^{2}(t)/U^{2}(0)$ in the region $\zeta >0.01\lambda$. As expected, the mean flow sees a real-time response with the onset of the gust. The peak streamwise turbulence variation $\langle u'u'\rangle _{p}=\max {\langle u'u'\rangle }$ changes rapidly after a short delay. When the wind is growing, an overshoot appears in $\langle u'u'\rangle _{p}$ after the gust ends, and its temporal evolution is qualitatively similar to the channel flow result (He & Seddighi Reference He and Seddighi2015). The other two components, $\langle w'w'\rangle _{p}$ and $\langle -u'w'\rangle _{p}$, experience relatively less change in the growing and decaying wind gust. Compared with the evolution of the mean flow energy, the turbulence peak values show a delayed response (also known as the memory effect) to the gust. Because of this phase delay, the turbulent kinetic energy (TKE) relative to the mean energy, i.e. $\langle u_i'u_i'\rangle /U^{2}$, decreases in the growing wind and increases in the decaying wind. Detailed explanations are provided by the turbulence budget analysis in the next subsection.

Figure 5. Time variation of the peak values of the turbulence variances in (a) the growing wind and (b) the decaying wind, in case CU8. Also plotted is the mean flow energy $U^{2}(t)$ at $\zeta =0.051\lambda$. All terms are normalized by their values at $t=0$.

3.2. Transient behaviours of a Reynolds stress budget and wave-induced effect

We perform a Reynolds stress budget analysis in this subsection to investigate the mechanism responsible for the turbulent stress evolution shown in § 3.1. In a previous study by Yang & Shen (Reference Yang and Shen2010), the budget of the TKE under the steady wind–wave condition was obtained. The detailed budget for the Reynolds stress components in unsteady wind fields over waves, however, was rarely reported. Following the work of Xuan, Deng & Shen (Reference Xuan, Deng and Shen2020) on Langmuir turbulence under waves, we can write the Reynolds stress budget equation of wind turbulence in the frame of reference moving with the wave as

(3.1)\begin{align} \frac{\partial\langle u_i' u_j'\rangle}{\partial t}&=\underbrace{-(U-c)\boldsymbol{\cdot}\boldsymbol{\nabla}\langle u_i' u_j'\rangle}_{A} +\underbrace{\nu\,\nabla^{2}\langle u_i'u_j'\rangle}_{D^{\nu}}\nonumber\\ &\quad \underbrace{{}-\langle u_i'u_k'\rangle\,\frac{\partial U_j}{\partial x_k}-\langle u_j'u_k'\rangle\,\frac{\partial U_i}{\partial x_k}}_{P^{m}} \underbrace{{}-\langle u_i'u_k'\rangle\,\frac{\partial \tilde{u}_j}{\partial x_k}-\langle u_j'u_k'\rangle\,\frac{\partial \tilde{u}_i}{\partial x_k}}_{P^{w}} \nonumber\\ &\quad \underbrace{{}-\frac{1}{\rho}\left(\frac{\partial \langle p'u_i'\rangle}{\partial x_j}+\frac{\partial \langle p'u_j'\rangle}{\partial x_i}\right)}_{T^{p}} +\underbrace{\left\langle\frac{p'}{\rho}\left(\frac{\partial u_i'}{\partial x_j}+\frac{\partial u_j'}{\partial x_i}\right)\right\rangle}_{R}\nonumber\\ &\quad \underbrace{{}-\frac{\partial \langle u_i'u_j'u_k'\rangle}{\partial x_k}-\frac{\partial \langle u_i'\tau_{jk}^{d}\rangle}{\partial x_k}-\frac{\partial \langle u_j'\tau_{ik}^{d}\rangle}{\partial x_k}}_{T^{t}} \underbrace{{}-2\nu\left\langle\frac{\partial u_i'}{\partial x_k}\,\frac{\partial u_j'}{\partial x_k}\right\rangle}_{\epsilon}, \end{align}

where $A$ is the advection, $D^{\nu }$ is the viscous diffusion, $P^{m}$ is the production by the mean shear, $P^{w}$ is the production by the wave-coherent motions, $T^{p}$ is the pressure transport, $R$ is the pressure–strain correlation, $T^{t}$ is the spatial flux associated with turbulence and SGS stress, and $\epsilon$ is the dissipation. For the non-stationary turbulence field with the wind gust over the wave, each budget term is a function of two spatial coordinates (i.e. $\xi$ and $\zeta$) and the time.

To investigate the time evolution of a budget term, we perform the integration along both the $\xi$ and $\zeta$ directions to evaluate its net contribution in the entire domain. For example, with the dissipation, we define

(3.2)\begin{equation} \breve{\epsilon}_{ij}(t)=\frac{1}{L_xL_z}\int \epsilon_{ij}(\xi,\zeta)\,\mathrm{d}\xi\,\mathrm{d}\zeta, \end{equation}

where the integral, denoted by $\breve {(\ )}$, is performed in the whole computational domain (He & Seddighi Reference He and Seddighi2013); alternatively, it can be computed separately in part of the domain (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2020), e.g. in the buffer layer and the logarithmic layer. The net contributions of the other budget terms ($\breve {P}_{ij}^{m}$, $\breve {R}_{ij}$, etc.) are calculated following (3.2).

The time history of the net contributions for $\langle u'u'\rangle$ is shown in figure 6, where the dominant budget terms are plotted and the secondary terms are neglected for simplicity. The rapid change in the production by the mean shear $\breve {P}_{11}^{m}$ is a key contributor to the time variations of $\langle u'u'\rangle$. In both growing wind and decaying wind, there is a clear phase lag in the pressure–strain correlation $\breve {R}_{11}$, compared with the production $\breve {P}_{11}^{m}$ and the dissipation $\breve {\epsilon }_{11}$. The net contribution of the normal components of the pressure–strain correlation terms, $\breve {R}_{11}$, $\breve {R}_{22}$ and $\breve {R}_{33}$, to the TKE vanishes because their role is to redistribute the energy among the different normal Reynolds stress terms through the pressure. The direction of the energy transfer is predominantly from $\langle u'u'\rangle$ to $\langle v'v'\rangle$ and $\langle w'w'\rangle$, as indicated by the negative $\breve {R}_{11}$ in figure 6. Therefore, the delay of $\breve {R}_{11}$ compared with $\breve {P}_{11}^{m}$ explains the phase lag of $\langle w'w'\rangle$ with respect to $\langle u'u'\rangle$ (figure 5).

Figure 6. Time variation of the integrated budget terms for $\langle u'u'\rangle$ for (a) the growing wind and (b) the decaying wind in case CU8. All terms are normalized by $\breve {P}_{11}^{m}(t=0)$.

To evaluate the wave-induced effect, we now examine the spatial distributions of the turbulence production associated with the mean shear, $P_{11}^{m}$, and the wave-coherent shear, $P_{11}^{w}$, defined in (3.1). As shown from comparison between figures 7(a) and 7(b), the magnitude of $P_{11}^{m}$ increases drastically during the growing wind gust in case CU8, while the spatial distributions are unaffected. In particular, the positive peak values of $P_{11}^{m}$ appear on the leeward side of the wave crest, which is similar to the steady wind turbulence in previous numerical studies (e.g. Yang & Shen Reference Yang and Shen2010). A similar feature was observed for the total turbulence production in the tank experiment by Yousefi, Veron & Buckley (Reference Yousefi, Veron and Buckley2021). Since the change in the turbulent stress is relatively small during the gust, the time variation of $P_{11}^{m}$ is caused primarily by the mean shear evolution near the wave surface (figure 22). In the decaying wind case, the qualitative features of $P_{11}^{m}$ are similar except that its magnitude decreases during the gust (results not plotted).

Figure 7. Turbulence production by the mean shear for the growing wind at (a) $t=0$ and (b) $t=T_g$, in case CU8. Note that the contour magnitudes are different in (a) and (b).

Figure 8 shows the spatial distribution of the turbulence production by the wave-coherent shear, $P_{11}^{w}$, in the growing wind. The wind gust induces a change in both the magnitude and spatial patterns of $P_{11}^{w}$, because of the variation in the wave-coherent shear $\partial {\tilde {u}_i}/\partial {x_j}$. The contours of $P_{11}^{w}$ then show a strong dependency on the wave age. While the net contribution of $P_{11}^{w}$ integrated over a wavelength at a constant $\zeta$ is small, it has a local magnitude of the same order as the mean-shear production $P_{11}^{m}$ (see e.g. figures 7b and 8b), which imposes a strong modulation on the local turbulence field by the wave. For the decaying wind, the evolution of $P_{11}^{w}$ is qualitatively a time reversal of figure 8 (not plotted). These transient features of $P_{11}^{w}$ are associated closely with those of the wave-coherent velocity field, with more detailed analyses presented in § 3.3 below.

Figure 8. Turbulence production by the wave-coherent motions for the growing wind in (a,b) case CU8, (c,d) case CU19, and (e, f) case CU42. Plotted are the normalized values $P_{11}^{w}/ku_{*,0}^{3}$ at $t=0$ and $t=T_g$.

The locally axisymmetric turbulence hypothesis is an important aspect in the characterization of turbulent flows. According to the theory developed by George & Hussein (Reference George and Hussein1991), one can validate the hypothesis by examining the relation between certain components of the mean square velocity gradient

(3.3)\begin{equation} \overline{\left(\frac{\partial u'_i}{\partial x_j}\right)^{2}}=\overline{\left(\frac{\partial u'_l}{\partial x_m}\right)^{2}}, \end{equation}

where $\overline {(\ )}$ denotes the (combined) ensemble and plane averaging, and $(i,j,l,m)=(1,2,1,3)$, $(2,1,3,1)$, $(2,3,3,2)$, $(2,2,3,3)$. Here, we assume that $x_1$ is the preferred local axis. To facilitate analysis, we first perform an interpolation in the vertical direction to transform the velocity field from the original boundary-fitted grid to a Cartesian grid where $x_1$, $x_2$ and $x_3$ correspond to the streamwise, spanwise and vertical directions. The velocity gradients are then calculated in regions above the wave crest. The ratios $\overline {({\partial u'_i}/{\partial x_j})^{2}}/\overline {({\partial u'_l}/{\partial x_m})^{2}}$ are qualitatively the same in all cases, and we present the results in the growing wind for case CU8 as an example in figure 9. For an exact axisymmetric turbulent flow, all these ratios should be unity. In our simulation, their profiles show strong deviations from one below $z\sim O(0.1\lambda )$ (figure 9a). The time history of the ratios at $z\approx 0.1\lambda$ is shown in figure 9(b). Because of the memory effect in the turbulence response to the wind gust event, most $\overline {({\partial u'_i}/{\partial x_j})^{2}}/\overline {({\partial u'_l}/{\partial x_m})^{2}}$ are nearly unchanged at this height, except for the one with $(i,j,l,m)=(1,2,1,3)$. We find that this phenomenon is caused mainly by a rapid change in $\overline {({\partial u'_1}/{\partial x_3})^{2}}$ that is associated closely with the streamwise streak variation during the wind gust (see figure 2). The ratios corresponding to $(i,j,l,m)=(1,2,1,3)$ and $(2,3,3,2)$ satisfy the locally axisymmetric turbulence hypothesis with a 5 % error, while for $(i,j,l,m)=(2,1,3,1)$ and $(2,2,3,3)$, the errors are larger (20 %–40 %). In general, the wind turbulence field shows a roughly local axisymmetry away from the wave surface. We remark that for wind turbulence over waves, the calculation of $\overline {({\partial u'_i}/{\partial x_j})^{2}}$ is non-trivial because of the wave-coherent motions (in the wind field) and the fluctuating surface, and therefore more systematic study is needed to examine whether (3.3) is a good indicator of a streaky turbulence structure.

Figure 9. Ratios of the mean square velocity gradient in the growing wind for case CU8: (a) vertical profiles at $t=0$; (b) time variations at height $z=0.1\lambda$.

We next consider the ratio between the turbulent stress components in the inner layer, which, according to Belcher & Hunt (Reference Belcher and Hunt1993), can be estimated from $T_D(z_i)=T_L(z_i)$. Here, $z_i$ is the top of the inner layer, $T_D(z_i)=z_i/u_*$ is the eddy turnover time scale, and $T_L(z_i)=k/(U(z_i)-c)$ is the advection time scale for wind flow past a wavelength. During the wind gust, the inner layer thickness decreases from $0.076\lambda$ to $0.036\lambda$ for the growing wind, and increases from $0.043\lambda$ to $0.087\lambda$ for the decaying wind, in case CU8. To ensure that the analysis is performed in the inner layer, here we choose the fixed height at $z=0.036\lambda$, which is close to the inner layer thickness found for wind turbulence at the wave age $6.5$ in a tank experiment (Buckley & Veron Reference Buckley and Veron2016). Figure 10 shows the time variation of the ratio of the turbulent normal stress to the turbulent shear stress, which was assumed to be constant in the inner layer by Belcher & Hunt (Reference Belcher and Hunt1993) in their study of steady wind over waves. Our result, at the equilibrium state ($t=0$), shows that the stress ratio is only slightly less than their recommended value. Similar to the turbulent stress, the stress ratio has a delayed response to the gust. While the time variation of $-\langle w'^{2}\rangle /\langle u'w'\rangle$ is relatively small, the wind gust impact is significant for $-\langle u'^{2}\rangle /\langle u'w'\rangle$ (note the different scalings of the vertical coordinates in figures 10a,b). In summary, we find that compared to statistically steady wind, the gust induces appreciable differences to the turbulence, and the delayed response to the gust event is a key feature of the transient behaviours of the turbulence field. The properties of the mean square velocity gradient indicate that the local axisymmetry may exist in the wind turbulence field away from the wave surface, while the (more restrictive) conditions for local isotropy are less likely to be satisfied.

Figure 10. Time variations of the ratio of the turbulent normal stress to the turbulent shear stress at height $z=0.036\lambda$ in case CU8. Also plotted are the values corresponding to an equilibrium boundary layer over a wave (Belcher & Hunt Reference Belcher and Hunt1993): $-\langle u'^{2}\rangle /\langle u'w'\rangle =6.3$ and $-\langle w'^{2}\rangle /\langle u'w'\rangle =1.7$, denoted by the black dashed line.

Figure 11. Wave-coherent streamwise velocity and its phase angle relative to the surface elevation for the growing wind in (a,b) case CU8, (c,d) case CU19, and (e, f) case CU42. The results are plotted at (a,c,e) $t=0$ and (b,d, f) $t=T_g$. Also listed is the instantaneous wave age $c/u_*$. Here, $\tilde {u}$ is normalized by $u_{*,0}$.

Figure 12. Same as figure 11, but for the vertical wave-coherent velocity $\tilde {w}$.

3.3. Evolution of wave-coherent motions

This subsection focuses on the transient features of the wave-coherent streamwise velocity $\tilde {u}$, vertical velocity $\tilde {w}$ and pressure $\tilde {p}$, which are useful for the quantitative analyses in the following subsections. We first illustrate $\tilde {u}$ and $\tilde {w}$ for the growing wind cases at the initial steady state $t=0$ and the end of the gust $t=T_g$, as plotted in figure 11. For all three cases, the contour patterns of the wave-coherent motions under the unsteady wind condition are qualitatively similar to those under the steady wind condition, as long as their wave ages are in the same range. In the small wave age range ($c/u_*\sim 10$, see figures 11a,b,d), $\tilde {u}$ shows the non-separated sheltering feature, which corresponds to a phase shift of ${\rm \pi} /2$ near the wave surface. The distribution of $\tilde {u}$ varies drastically with height near the wave surface. Because the wave ages throughout the gust event in case CU8 are in the same range, the contour shapes of $\tilde {u}$ are qualitatively unchanged during the gust despite a significant increase in its magnitude (comparing figures 11a,b). In case CU19, the wave-coherent motions at $t=0$ decrease rapidly with height, and overall their magnitudes are relatively small in most of the domain. For the wave age in this case, the gust event has a profound effect on both the magnitude and the contour pattern of the wave-coherent motions. Most notably, the contours of $\tilde {u}$ near the wave surface are tilted towards the $-x$ direction at $t=0$, and the $+x$ direction at $t=T_g$ (figures 11c,d). In the intermediate wave age, i.e. $c/u_*\sim O(20)$, the contour patterns correspond to a phase shift of $2{\rm \pi}$ with the height from $z=0$ to $z=0.05\lambda$ (see figures 11c, f). Note that the discontinuity near $z=0.01\lambda \unicode{x2013}0.02\lambda$ is due to the definition of the phase angle, which limits $\phi _{\tilde {u}}$ between $-{\rm \pi}$ and ${\rm \pi}$. For the fast propagating wave at $t=0$ in case CU42 (see figure 11e), $\tilde {u}$ is nearly symmetric with respect to the wave surface except for a thin layer near the wave surface (for a similar pattern, see Cao & Shen Reference Cao and Shen2021).

The distributions of $\tilde {w}$ in the growing wind are plotted in figure 12. Compared with the streamwise velocity $\tilde {u}$ with complex contour patterns, $\tilde {w}$ exhibits relatively simple spatial distributions. Regardless of the wave age and the time instant, the imaginary part of $\tilde {w}$ always dominates, i.e. $\phi _{\tilde {w}}\approx \pm {\rm \pi}/2$ in the entire domain. When the wave age is not large (figures 12a,b,d), $\phi _{\tilde {w}}$ changes from $-{\rm \pi} /2$ to ${\rm \pi} /2$, while for the relatively large wave ages in figures 12(c,e, f), $\phi _{\tilde {w}}$ has a constant value close to $-{\rm \pi} /2$ in the region $z<0.1\lambda$. This feature of $\tilde {w}$ agrees with the Miles (Reference Miles1957) theory based on the Rayleigh equation, which predicts the phase change across the critical layer. A brief review of the critical layer at the steady state can be found in Appendix C.

The evolution of the amplitudes of the wave-coherent motions in the growing wind is shown in figure 13. The vertical coordinates for $|\tilde {u}|$ and $|\tilde {w}|$ are plotted on a logarithmic scale to highlight the near-surface structures. In case CU8, $|\tilde {u}|$ and $|\tilde {w}|$ reach their peak values around $\zeta \sim O(10^{-2})\lambda$. In case CU42, $|\tilde {u}|$ and $|\tilde {w}|$ decrease nearly monotonically with the height except for $|\tilde {u}|$ at early times. In case CU19, the vertical distribution of the wave-coherent velocity amplitudes sees a transition from a monotonic function of height to a more complex pattern during the gust. The pressure amplitudes $|\tilde {p}|$ on the logarithmic scale are nearly straight lines with the same slope in most cases, following an $\exp (-k\zeta )$ decay with the height, especially when $\zeta$ is large. The only exceptions occur at $t=0$ in case CU19 and $t=T_g$ in case CU42, where $|\tilde {p}|$ is prominent only at the wave surface. Away from the surface, the values of $|\tilde {p}|$ are too small and may be subject to fluctuations in the triple decomposition statistics. The velocity amplitudes show a similar exponential decay above $\zeta \sim O(10^{-2})\lambda$. The exponential decay of the wave-coherent motions is a well-known feature of the asymptotic inviscid solutions at large $z$ (Belcher & Hunt Reference Belcher and Hunt1993). In contrast, a noticeable deviation from $\exp (-k\zeta )$ may be observed in the distribution of the pressure amplitude near the wave surface because of the wave geometry and the strong gradient in the vertical velocity, as suggested by previous studies (e.g. Yang & Shen Reference Yang and Shen2010; Grare et al. Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013). As the mean velocity increases rapidly during the gust, the viscous effect is enhanced mostly in the thin region near the wave surface, as shown in figure 3. But the outer region remains approximately inviscid, and as a result, the exponential decay still holds there.

Figure 13. Vertical distributions of the amplitudes of the wave-coherent motions $\tilde {u}$, $\tilde {w}$ and $\tilde {p}$ for the growing wind in (a–c) case CU8, (d–f) case CU19, and (g–i) case CU42. The exponential decay $\exp (-k\zeta )$ is denoted by the black dashed lines. The results are plotted at $t/T_g=0, 0.5, 1.0$. All terms are normalized using the wall unit at $t=0$.

To identify the wave boundary layer, we plot separately the profiles of the turbulent and wave-coherent parts of the Reynolds shear stress (figure 14). The sign of $-\langle \tilde {u}\tilde {w}\rangle$ depends on the instantaneous wave age. For slow waves in case CU8 (figures 14a,d), the wave-coherent stress has an opposite sign to the turbulent part, which agrees with the findings of Hara & Sullivan (Reference Hara and Sullivan2015). For fast waves in case CU42 (figures 14c, f), their signs are the same. The transition between these two patterns occurs at the wave age $c/u_*\sim 20$, as indicated by the results in case CU19 (figures 14b,e). When $-\langle \tilde {u}\tilde {w}\rangle$ is sufficiently large compared with $-\langle u'w'\rangle$ (figures 14c,d), the wave boundary layer is most distinguishable, with thickness $0.008\lambda \unicode{x2013}0.02\lambda$ (see also the tank experiment by Buckley & Veron Reference Buckley and Veron2016).

Figure 14. Turbulent and wave-coherent parts of the Reynolds shear stress for the growing wind at (a–c) $t=0$ and (d–f) $t=T_g$. All terms are normalized using the wall unit at $t=0$. The instantaneous wave ages for (a–f) are $c/u_*=8.3$, $19.2$, $41.6$, $4.5$, $10.4$, $22.3$, respectively.

To summarize, the wave-coherent motions can be identified clearly throughout our simulations, which justifies the use of the triple decomposition in the transient wind–wave process. Owing to page limits for the paper, the time series of the two-dimensional distributions of all wave-coherent components in both growing wind and decaying wind are depicted in the supplementary movies available at https://doi.org/10.1017/jfm.2022.577. A key feature that we see from figures 11–13 is that while the turbulent stresses exhibit the memory effect associated with the gust events (§ 3.2), the wave-coherent motions have a quasi-stationary response to the wind gust and they are strongly dependent on the instantaneous wave age. Furthermore, the growing wind gust has a counterintuitive phenomenon in case CU42: the wave-coherent motions are suppressed, i.e. their magnitudes at $t=T_g$ (figures 11fand 12f) are smaller than their values at $t=0$ (figures 11e and 12e), while in a decaying wind their amplitudes increase. This phenomenon is explained in the following analysis based on the momentum equation.

3.4. Analysis of the wave-coherent momentum equation

To understand further the mechanisms that govern the evolution of the wave-coherent motions during wind gust, we perform more detailed analyses by resorting to the momentum equation. While an energy budget equation for the wave-coherent motions can also be derived in a way similar to the TKE equation (3.1), the phases of the source terms would be complicated because of the product of trigonometric functions. The momentum equation, on the other hand, has a simpler form and preserves directly the phases of the wave-coherent motions. By performing the triple decomposition, the vertical wave-coherent momentum equation can be written as

(3.4)\begin{equation} \frac{\partial\tilde{w}}{\partial t}=\underbrace{-(U-c)\,\frac{\partial\tilde{w}} {\partial \xi}}_{\mathcal{A}} -\underbrace{\frac{\partial\widetilde{u'w'}}{\partial \xi}}_{\mathcal{T}} -\underbrace{\frac{1}{\rho}\,\frac{\partial \tilde{p}}{\partial \zeta}}_{\mathcal{P}} +\underbrace{\nu\,\nabla^{2}\tilde{w}}_{\mathcal{V}}, \end{equation}

where $\mathcal {A}$, $\mathcal {T}$, $\mathcal {P}$ and $\mathcal {V}$ denote the effects of advection, turbulence, pressure gradient and viscosity, respectively. The effect of the wave-coherent stresses $\tilde {u}_i\tilde {u}_j$ ($i,j=1,3$) is omitted because they have a direct impact on the second harmonic wavenumber $4{\rm \pi} /\lambda$ instead of the primary wavenumber $2{\rm \pi} /\lambda$, and thus have little effect on the evolution of $\tilde {w}$. Apparently, the advection term $\mathcal {A}$ and the pressure gradient term $\mathcal {P}$ have rapid responses to the wind gust because they are determined by the mean velocity and the wave-coherent motions. The turbulence term $\mathcal {T}$, on the other hand, shares a similar memory effect with $\widetilde {u'w'}$.

Similar to the decomposition of the wave-coherent motions shown in (2.6), the real part and the imaginary part of each term in the vertical momentum equation (3.4) can be calculated. Figure 15 shows the results for the real parts at $t=0.5T_g$. Note that the viscous term is found to be negligibly small and has been omitted. In all three cases, the pressure gradient and the advection terms are significant in balancing the real part of (3.4). The turbulence effect nearly vanishes in case CU42. For cases CU8 and CU19, the effect of the turbulence term $\mathrm {Re}[\mathcal {T}]$ is relatively strong in the decaying wind, while in the growing wind it is relatively weak. This phenomenon is similar to the time variation of the turbulence intensity discussed in § 3.1, and is also associated with the memory effect of the turbulent stress, which leads to a phase lag in $\mathrm {Re}[\mathcal {T}]$ with respect to the rapidly varying $\mathrm {Re}[\mathcal {A}]$ and $\mathrm {Re}[\mathcal {P}]$. The temporal growth rate $\mathrm {Re}[\partial \tilde {w}/\partial t]$ is negligibly small in all cases, suggesting that the real parts are largely balanced regardless of the gust event.

Figure 15. Real part of the terms in the vertical wave-coherent momentum equation for (a–c) the growing wind and (d–f) the decaying wind, at $t=0.5T_g$. The instantaneous wave ages are $c/u_*=5.1, 11.7, 25.2, 7.1, 16.6, 35.2$ in (a–f), respectively. All terms are normalized by $ku_{*,0}^{2}$.

In contrast to the relatively simple distributions of the real parts of the terms in (3.4), the distributions of the imaginary parts are more complex, as shown in figure 16. The imaginary part of the viscous term, compared with its real part, is important in the vertical momentum equation near the wave surface, because $\mathrm {Im}[\mathcal {V}]=\nu \,\nabla ^{2}\mathrm {Im}[\tilde {w}]$, and $\mathrm {Im}[\tilde {w}]$ dominates over $\mathrm {Re}[\tilde {w}]$ (figure 12). In general, for the imaginary part, the values of different terms in (3.4) are of the same order, especially in case CU8. Unlike $\mathrm {Re}[\partial \tilde {w}/\partial t]$, which is close to zero, $\mathrm {Im}[\partial \tilde {w}/\partial t]$ becomes positive in the growing wind and negative in the decaying wind. Hence $\mathrm {Im}[\partial \tilde {w}/\partial t]$ is the main contributor to the contour pattern changes presented in figure 12. In particular, for the growing wind in case CU42, $\mathrm {Im}[\tilde {w}]$ at $t=0$ is negative, noticing that $\mathrm {Im}[\tilde {w}]=|\tilde {w}|\sin \phi _{\tilde {w}}<0$ when $\phi _{\tilde {w}}\approx -{\rm \pi} /2$ (see figure 12e). The positive $\partial \mathrm {Im}[\tilde {w}]/\partial t$ then causes a magnitude reduction in $\mathrm {Im}[\tilde {w}]$ during the growing wind event, which is confirmed in figures 12(e, f).

Figure 16. Same as figure 15, but for the imaginary parts.

While the time variation in $\tilde {w}$ is dominated by $\mathrm {Im}[\partial \tilde {w}/\partial t]$, the real and imaginary parts of the momentum equation (3.4) are not independent of each other but dynamically coupled through the differential operator $\partial /\partial \xi$. For the advection term and the turbulence term, we have

(3.5)$$\begin{gather} \mathrm{Re}[\mathcal{A}]=k(U-c)\,\mathrm{Im}[\tilde{w}], \end{gather}$$

(3.6)$$\begin{gather}\mathrm{Im}[\mathcal{A}]={-}k(U-c)\,\mathrm{Re}[\tilde{w}], \end{gather}$$

(3.7)$$\begin{gather}\mathrm{Re}[\mathcal{T}]= k\,\mathrm{Im}[\widetilde{u'w'}], \end{gather}$$

(3.8)$$\begin{gather}\mathrm{Im}[\mathcal{T}]={-} k\,\mathrm{Re}[\widetilde{u'w'}]. \end{gather}$$

In case CU42, for instance, $\mathrm {Re}[\mathcal {A}]$ (figures 15c, f) is much larger than $\mathrm {Im}[\mathcal {A}]$ (figures 16c, f) because the wave-coherent vertical velocity is dominated by the out-of-phase imaginary part $\mathrm {Im}[\tilde {w}]$ (see figure 12).

Next, we perform a quantitative analysis of the wave-coherent motions by utilizing and comparing with the viscous curvilinear model proposed by Cao et al. (Reference Cao, Deng and Shen2020). The model originates from the triple decomposition of the Navier–Stokes equations for wind turbulence over a monochromatic wave under the steady wind–wave condition. In their derivation, the nonlinear terms associated with the turbulent stress and the second-order terms of $O((ka)^{2})$ related to the grid transformation are neglected. The governing equation for the vertical wave-coherent velocity $\tilde {w}$ reads

(3.9)\begin{align} &\displaystyle {-\frac{\nu}{ik}}\left(\frac{\mathrm{d}^{4}}{\mathrm{d}\zeta^{4}} -2k^{2}\,\frac{\mathrm{d}^{2}}{\mathrm{d}\zeta^{2}}+k^{4}\right)\mathcal{F}[\tilde{w}] +\left[(U-c)\left(\frac{\mathrm{d}^{2}}{\mathrm{d}\zeta^{2}}-k^{2}\right) -\frac{\mathrm{d}^{2}U}{\mathrm{d}\zeta^{2}}\right]\mathcal{F}[\tilde{w}] \nonumber\\ &\quad \displaystyle =\mathcal{F}\left[\frac{\partial }{\partial t}\left(\frac{\partial\tilde{u}}{\partial\zeta}-\frac{\partial\tilde{w}}{\partial\xi}\right)\right] +\nu\,\mathcal{F}[\eta]\,\frac{\mathrm{d}^{2}}{\mathrm{d}\zeta^{2}} \left[\left(\frac{\zeta}{h}-1\right)\frac{\mathrm{d}^{2}U}{\mathrm{d}\zeta^{2}}\right], \end{align}

where $\mathcal {F}$ denotes the Fourier transform operator. For completeness, we keep the time derivative term that was neglected by Cao et al. (Reference Cao, Deng and Shen2020). But its actual value is also found to be insignificant compared with the other terms in the present wind gust cases. The mean profile $U(\zeta )$ is provided by the simulation, and the boundary conditions are

(3.10a,b)$$\begin{gather} \tilde{w}(\zeta=0)=w_s,\quad \tilde{w}(\zeta\rightarrow\infty)=0, \end{gather}$$

(3.11a,b)$$\begin{gather}\left.\frac{\mathrm{d}\tilde{w}}{\mathrm{d}\zeta}\right|_{\zeta=0} ={-}\frac{\mathrm{d}u_s}{\mathrm{d}\xi}+\frac{\mathrm{d}U}{\mathrm{d}\zeta}\,\frac{\mathrm{d}\eta}{\mathrm{d}\xi},\quad \left.\frac{\mathrm{d}\tilde{w}}{\mathrm{d}\zeta}\right|_{\zeta\rightarrow\infty}=0, \end{gather}$$

where $u_s$ and $w_s$ are the streamwise and vertical wave orbital velocities, respectively (see Appendix A).

Here, the Dirichlet boundary condition for $\tilde {w}$ at $\zeta =0$ is prescribed by the wave orbital velocity, and its derivative ${\mathrm {d}\tilde {w}}/{\mathrm {d}\zeta }$ is obtained from the continuity equation. In practice, the top boundary condition ($\zeta \rightarrow \infty$) can be evaluated at a height where the magnitude of $\tilde {w}$ is sufficiently small, say $\zeta =0.5\lambda$ (Cao et al. Reference Cao, Deng and Shen2020). Near the wave surface, the imaginary part of the vertical wave-coherent velocity dominates over the real part because of the ${\rm \pi} /2$ phase difference between the vertical wave orbital velocity and the surface elevation (see (A6)).

In figure 17, we present $\mathrm {Im}[\tilde {w}]$ at three instants in the growing wind gust. The evolution of $\mathrm {Im}[\tilde {w}]$ dominates the qualitative features shown in figure 12. For example, the phase change of ${\rm \pi}$ in the vicinity of the wave surface in case CU8 (figure 12a) corresponds to the sign reversal at small $\zeta /\lambda$ in figure 17(a). Consistent with the positive $\partial \mathrm {Im}[\tilde {w}]/\partial t$ shown in figure 17, $\mathrm {Im}[\tilde {w}]$ increases with the growing wind in all cases (figure 12). In general, the LES result agrees well quantitatively with the model prediction at both the steady state ($t=0$) and the non-equilibrium state ($t=0.5T_g$ and $t=T_g$). Hence the quasi-linear mechanism (Cao et al. Reference Cao, Deng and Shen2020) that governs the response of $\mathrm {Im}[\tilde {w}]$ to the mean velocity profile $U$ is effective even under the gusty condition.

Figure 17. Imaginary part of the vertical wave-coherent velocity for the growing wind in (a) case CU8, (b) case CU19, and (c) case CU42. Also plotted are the corresponding results computed from the viscous curvilinear model, denoted by the black solid lines.

For $\mathrm {Re}[\tilde {w}]$, on the other hand, the difference between the model prediction and the LES result is noticeable. In figure 18, we plot the instantaneous result at $t=0.5T_g$ as an example. Qualitatively, the result predicted by the model has a shape similar to the LES result, in which $\mathrm {Re}[\tilde {w}]$ first increases with the height and then decreases after reaching a positive peak value. The modelling of $\mathrm {Re}[\tilde {w}]$ is challenging because it has a small magnitude compared with $\mathrm {Im}[\tilde {w}]$. At $t=0.5T_g$, the peak value of $\mathrm {Re}[\tilde {w}]$ is approximately one-fifth of the corresponding value of $\mathrm {Im}[\tilde {w}]$ in case CU8 (comparing figures 17a and 18a). At larger wave ages in case CU42, the difference between the magnitudes of $\mathrm {Im}[\tilde {w}]$ and $\mathrm {Re}[\tilde {w}]$ can be as large as two orders of magnitude (figures 17c and 18c). The viscous curvilinear model fails to predict the weak $\mathrm {Re}[\tilde {w}]$ because the impacts of the turbulence and second-order terms omitted in the derivation of the model are no longer negligible.

Figure 18. Real part of the vertical wave-coherent velocity for the growing wind in (a) case CU8, (b) case CU19, and (c) case CU42. The results computed from the viscous curvilinear model are denoted by the black solid lines.

From the analyses in this subsection, we obtain a clear picture of the wave-coherent motions. First, compared with the turbulent stress, the wave-coherent motions respond instantly to the gust event and show strong dependency on the wave age. Besides, the gust event leads to a drastic change in the out-of-phase component $\mathrm {Im}[\tilde {w}]$, and has little impact on the in-phase component $\mathrm {Re}[\tilde {w}]$. In the next subsection, we investigate whether these phenomena are responsible for producing the drag hysteresis that has been observed in laboratory and field experiments (Waseda et al. Reference Waseda, Toba and Tulin2001; Hwang et al. Reference Hwang, García Nava and Ocampo-Torres2011) and how the change in $\mathrm {Im}[\tilde {w}]$ contributes to the variation in wind–wave energy transfer.

3.5. Wind gust impact on wave growth

It is well recognized that the wave-coherent motions play a critical role in the wind–wave momentum and energy transfer. In this subsection, we investigate how their time variations affect these processes. We first examine the evolution of the form drag $\tau _p$ and viscous drag $\tau _{\nu }$. In the present study, we use the sign convention such that the drags are positive when the direction of the momentum flux is from the wind to the wave, such as at small wave ages. As shown in figure 19, the viscous drag $\tau _{\nu }$ changes continuously during the gust and remains positive. In all cases, $\tau _{\nu }$ at the end of the growing wind gust is slightly larger than that at the start of the decaying wind gust, while at the end of the decaying gust, $\tau _{\nu }$ nearly vanishes and is significantly smaller than its initial value in a growing wind. This hysteresis effect exhibited by $\tau _{\nu }$ is associated primarily with the evolution of the mean velocity profile and the mean shear. As discussed in § 3.1, the uniform change in the mean velocity caused by the external force occurs mainly in the outer region where the mean shear is expected to be roughly unchanged during the gust. The near-wall region is dominated by viscosity such that the mean velocity and its gradient are affected by their past states.

Figure 19. Evolution of form drag and viscous drag during the wind gust ($0< t< T_g$). The start ($t=0$) and end ($t=T_g$) of the gust are denoted by dash marks (–) and circles, respectively, and the arrows denote the direction of time increase.

The form drag $\tau _p$ sees an abrupt change at the start and end of the gust. During the wind gust, the variations in $\tau _p$ are relatively small at first, and this feature is shared in all cases. Near the end of the gust, there is a substantial change in $\tau _p$ for cases CU8 and CU19, while in case CU42 the change in $\tau _p$ is trivial. To understand the behaviours of $\tau _p$, we consider its definition:

(3.12)\begin{equation} \tau_p=\frac{1}{L_x}\int_0^{L_x}\frac{p\eta_x}{\rho}\,\mathrm{d}x =\frac{ka\left.\mathrm{Im}[\tilde{p}]\right|_{\zeta=0}}{\rho}. \end{equation}

The pressure at the wave surface can be obtained by integrating the imaginary part of the vertical momentum equation (3.4). Keeping only the dominant terms, we can rewrite (3.12) as

(3.13)\begin{equation} \tau_p=ka\int_0^{\lambda} k(U-c)\,\mathrm{Re}[\tilde{w}]\,\mathrm{d}\zeta +ka\int_0^{\lambda}\frac{\partial}{\partial t}\mathrm{Im}[\tilde{w}]\,\mathrm{d}\zeta. \end{equation}

The contribution to the form drag is then decomposed into two parts, the advection of $\mathrm {Re}[\tilde {w}]$ and the time variation of $\mathrm {Im}[\tilde {w}]$ caused by the wind gust. For a wind turbulence field at the steady state, only the advection term needs to be included in the calculation, while the gustiness term reduces to zero (see e.g. Cao et al. Reference Cao, Deng and Shen2020).

When the wind is unsteady, the gustiness term becomes important. In figure 20, we plot the time-averaged contributions from the advection and gustiness terms to the form drag during $0< t< T_g$. For the advection term, the contribution varies with the wave age, which determines the advection velocity $U-c$. For the contribution by the gustiness, however, the wave age has less influence because the magnitudes of $\partial \mathrm {Im}[\tilde {w}]/\partial t$ are comparable in different cases (see figure 16). Figure 20 shows that in all cases, the gustiness term is comparable to (cases CU8 and CU19) or stronger than (case CU42) the advection term. Therefore, the out-of-phase component $\mathrm {Im}[\tilde {w}]$ that changes rapidly during the gust plays an at least equally important role as the in-phase component $\mathrm {Re}[\tilde {w}]$ in the wind–wave momentum transfer.

Figure 20. Contributions of the advection and gustiness terms to the form drag in (a) the growing wind and (b) the decaying wind. All terms are calculated by taking the time average during $0< t< T_g$.

We now revisit figure 19 to explain the behaviour of $\tau _p$ based on the decomposition (3.13). The abrupt change immediately after $t=0$ and before $t=T_g$ is caused by the contribution of the gustiness term. During $0< t< T_g$, cases CU8 and CU19 have a noticeable change in $\tau _p$ because the advection term itself is changing with the mean velocity, and its contribution is comparable to that of the gustiness term (figure 20). For case CU42, although the mean velocity and the advection term change substantially during the wind gust, the magnitude of the advection term is much smaller than that of the gustiness term, as shown in figure 20. Therefore, the value of $\tau _p$ is relatively invariant during this process.

To quantify the amount of wind energy transferred to wave, we calculate the wave growth rate $\beta$ (Miles Reference Miles1957) and the temporal energy growth rate $\gamma$ (Donelan & Pierson Reference Donelan and Pierson1987):

(3.14)$$\begin{gather} \beta=\frac{2\tau_p}{(aku_*)^{2}}, \end{gather}$$

(3.15)$$\begin{gather}\gamma=\left(\frac{u_*}{c}\right)^{2}\beta=\frac{2\tau_p}{(akc)^{2}}. \end{gather}$$

In table 3, we list the time-averaged results of $\beta$ and $\gamma$ during the wind gust ($0< t< T_g$). We follow Miles & Ierley (Reference Miles and Ierley1998) to obtain $\bar {\beta }$ by calculating separately the mean form drag and the average of the square of the friction velocity, and then using them to calculate $\bar {\beta }$. For both $\bar {\beta }$ and $\bar {\gamma }$, there is a distinct difference between the growing wind result and the decaying wind result. The primary reason for this difference is due to the hysteresis phenomenon in the variation of the form drag and the viscous drag (figure 19).

Table 3. Time-averaged $\beta$ and $\gamma$ during $0< t< T_g$. Here, $\overline {(\ )}$ is the time-averaging operator. Also listed is the strength of the gust event $|\Delta U/U|/\omega T_g$, where $\Delta U/U=[U(0.02\lambda,T_g)-U(0.02\lambda,0)]/U(0.02\lambda,0.5T_g)$ is the relative change in the mean profile, and $\omega$ is the wave frequency.

We plot $\beta$ as a function of the wave age $c/u_*$ in figure 21(a). The results denoted by the light grey and dark grey symbols are instantaneous values at $t=0$, while those during the gust, denoted by coloured symbols, are calculated as the time-averaged values (see table 3). At the steady state before the gust event occurs, $\beta$ varies with the wave age and is in general consistent with the results in the literature (Sullivan et al. Reference Sullivan, Mcwilliams and Moeng2000; Kihara et al. Reference Kihara, Hanazaki, Mizuya and Ueda2007; Yang & Shen Reference Yang and Shen2010; Druzhinin et al. Reference Druzhinin, Troitskaya and Zilitinkevich2012; Åkervik & Vartdal Reference Åkervik and Vartdal2019). When the wind gust starts, there is an immediate change in the wave age associated with the time-varying friction velocity $u_*$ (see the evolution of the friction drag in figure 19). The shift in the wave age has been incorporated into the wind input source terms in phase-averaged wave models to account for the gustiness effect on the wave energy growth (Janssen Reference Janssen2004). For atmospheric variability with a time scale much larger than the wave period and the wind turbulence in a quasi-equilibrium state, this simple approximation is likely appropriate. But more generally, especially for a rapidly changing wind field, this is not the case. Our result in figure 21(a) shows a noticeable deviation of $\bar {\beta }$ from the steady state value. The strength of the wind gust can be measured roughly by $|\Delta U/U|/\omega T_g$ (table 3), or, more generally, $(\partial U/\partial t)/\omega |U|$ (see Miles & Ierley Reference Miles and Ierley1998). Note that while by this definition the wind gust event is the weakest in the high wave age case CU42, the deviation of $\beta$ from the steady-state values is most significant in this case because the gustiness effect has the largest relative contribution to the form drag (see figure 20).

Figure 21. (a) Wave growth rate $\beta$ scaled with the friction velocity $u_*$. (b) Temporal growth rate $\gamma$ scaled with the mean wind velocity at the height of half-wavelength $U_{\lambda /2}$. In (a), the results in cases CU8, CU19 and CU42 are denoted by triangles, squares and pentagons, respectively. The symbols with the filled colours light grey and dark grey denote the initial steady-state results, and those with the colours gold and blue are the time-averaged values $\bar {\beta }$ during the growing wind and the decaying gust, respectively. Also plotted are the steady-state results with similar wave steepness ($ak<0.1$) obtained using direct numerical simulation (Sullivan et al. Reference Sullivan, Mcwilliams and Moeng2000; Kihara et al. Reference Kihara, Hanazaki, Mizuya and Ueda2007; Yang & Shen Reference Yang and Shen2010; Druzhinin et al. Reference Druzhinin, Troitskaya and Zilitinkevich2012) and LES (Åkervik & Vartdal Reference Åkervik and Vartdal2019), denoted by circles. In (b), only those results with positive values of $\gamma$ and $U_{\lambda /2}/c-1$ are plotted. The black dashed line denotes the parametrization proposed by Donelan et al. (Reference Donelan, Babanin, Young and Banner2006): $\gamma =0.17(U_{\lambda /2}/c-1)|U_{\lambda /2}/c-1|$.

In figure 21(b), we plot $\gamma$ using an alternative scaling based on $U_{\lambda /2}$, the characteristic mean velocity at $\zeta =\lambda /2$ (Donelan et al. Reference Donelan, Babanin, Young and Banner2006). The steady-state values are consistent with the parametrization proposed by Donelan et al. (Reference Donelan, Babanin, Young and Banner2006), $\gamma =0.17$. While this parametrization is effective in reducing the scattering compared with the scaling based on $u_*$ shown in figure 21(a), it underestimates $\gamma$ in the growing wind and overestimates $\gamma$ in the decaying wind. The deviation, which can be up to several orders of magnitude, is a consequence of the change in the mean velocity $U_{\lambda /2}$ and $\gamma$ itself. Therefore, the non-stationary effect associated with the wind gust cannot be incorporated into such parametrized wind input models by simply using a shift in the wave age $c/U_{\lambda /2}$. More research on the model development, which is beyond the scope of this work, is called for in future studies.