Wind-Induced Changes to Surface Gravity Wave Shape in Deep to Intermediate Water

Wave shape (i.e. skewness or asymmetry) plays an important role in beach morphology evolution, remote sensing, and ship safety. Wind's influence on ocean waves has been extensively studied theoretically in the context of growth, but most theories are phase averaged and cannot predict wave shape. Most laboratory and numerical studies similarly focus on wave growth. A few laboratory experiments have demonstrated that wind can change wave shape, and two-phase numerical simulations have also noted wind-induced wave shape changes. However, wind's effect on wave shape is poorly understood, and no theory for it exists. For weakly nonlinear waves, wave shape parameters are the phase of the harmonic relative to the primary frequency (or harmonic phase HP, zero for a Stokes wave) and relative amplitude of the harmonic to the primary. Here, surface pressure profiles (denoted Jeffreys, Miles, and Generalized Miles) are prescribed based on wind-wave generation theories. Theoretical solutions are derived for quasi-periodic progressive waves and the wind-induced changes to HP, relative harmonic amplitude, as well as already known phase speed changes and growth rates. The wave shape parameters depend upon the chosen surface pressure profile, pressure magnitude and phase relative to the wave profile, and the nondimensional depth. Wave asymmetry is linked to the nondimensional growth rate. Atmospheric large eddy simulations constrain pressure profile parameters. HP predictions are qualitatively consistent with laboratory observations. This theory, together with the observables of HP and relative harmonic amplitude, can provide insight into the actual wave surface pressure profile.


Introduction
The shape of surface gravity waves plays a role in many physical phenomena. Wave shape is described by the third-order statistical moments, skewness and asymmetry (e.g. Hasselmann 1962;Elgar et al. 1990). These two parameters are integral in determining sediment transport direction (onshore vs. offshore) and magnitude (e.g. Drake & Calantoni 2001;Hsu & Hanes 2004;Gonzalez-Rodriguez & Madsen 2007), which play key roles in beach morphodynamics (e.g. Hoefel & Elgar 2003;Grasso et al. 2011). Wave shape is also pertinent in remote sensing, where wave skewness modulates the returned waveform in radar altimetry (e.g. Jackson 1979;Hayne 1980;Huang et al. 1983) and wave asymmetry affects the thermal emissions measured in polarimetric radiometry (e.g. Kunkee & Gasiewski 1997;Piepmeier & Gasiewski 2001;Johnson & Cai 2002). Additionally, these wave shape plays a role in determining ship response to wave impacts (e.g. Soares et al. 2008;Oberhagemann et al. 2013). Waves propagating on a flat bottom are ordinarily symmetric, though a number of processes can create asymmetry. While some wave asymmetry inducing-phenomena, such as wave shoaling (e.g. Elgar & Guza 1985, 1986, are well understood, the effect of wind on wave shape is still poorly understood. The influence of wind on ocean waves has been extensively studied, although primarily in the context of wave growth. An initial investigation by Jeffreys (1925) was based on a sheltering hypothesis where separated airflow resulted in reduced pressure on the leeward side of the wave, causing wave growth. While conceptually simple, this mechanism fell out of favour because such separation only occurs near breaking (Banner & Melville 1976) and is therefore unlikely to contribute meaningfully to wave growth (Young 1999). However, Jeffreys's theory has inspired some recent work: Belcher & Hunt (1993) developed a fully turbulent model wherein the sheltering effect causes a thickening of the boundary layer and wave growth, even without separation. Later treatments utilized different physical mechanisms such as resonant forcing by incoherent, turbulent eddies (Phillips 1957), vortex forcing from vertically sheared airflow (e.g. Miles 1957;Lighthill 1962), and non-separated sheltering (e.g. Belcher & Hunt 1993). Janssen (2004) provides an extensive overview of the relevant developments in wind-wave generation theory. When deriving energy and momentum fluxes from air to water, all of these seminal theories on wave growth (e.g. Phillips 1957;Lighthill 1962;Belcher & Hunt 1993) utilized a phase-averaging technique, which removes wave-shape information. Thus, although these theories of wind-wave interaction focused on the wave growth rate, no theoretical work has investigated the effect of wind on wave shape in a physically consistent manner.
Measurements and numerical simulations have also been used to investigate the dependence of wave growth on wind speed. Field measurements (e.g. Longuet-Higgins 1962;Snyder 1966;Hasselmann et al. 1973) and laboratory experiments (e.g. Shemdin & Hsu 1967;Plant & Wright 1977;Mitsuyasu & Honda 1982) have been used to parameterize how quickly intermediate and deep-water waves grow under various wind conditions, including short fetch (e.g. Lamont-Smith & Waseda 2008) and strong wind conditions (e.g. Troitskaya et al. 2012). Note, direct measurements of wave surface pressure (related to growth) are notoriously difficult (e.g. Donelan et al. 2005). Similarly, numerical simulations have also been used to predict wind-induced growth rates. Early atmospheric numerical models used the Reynolds-averaged Navier Stokes equations (RANS, e.g. Gent & Taylor 1976;Al-Zanaidi & Hui 1984) to calculate the energy loss of the wind field, however these early simulations could only approximate turbulence through a Reynolds-averaging process. Recent studies have analysed the turbulence behaviour in detail. Particle image velocimetry (PIV) and laser-induced fluorescence (LIF) have been used for turbulence measurements in laboratory experiments and have revealed turbulent structures above the waves (e.g. Veron et al. 2007;Buckley & Veron 2017, 2019. This turbulent behaviour has also been captured through direct numerical simulations (DNS) of the governing equations (e.g. Yang & Shen 2009 and by parameterizing subgridscale processes in large eddy simulations (LES, e.g. Yang et al. 2013;Hara & Sullivan 2015;Hao et al. 2018). When solving for the atmospheric dynamics, many of these simulations prescribe a static sinusoidal wave shape while focusing on the evolution of the wind field, as well as energy and momentum transfers. Therefore, any wind-induced changes to wave shape were not captured.
While there has been much research regarding wind-induced wave growth, wave shape has seen relatively little work. Two-phase (air and water) numerical simulations have begun incorporating dynamically evolving waves into their analyses (e.g. Liu et al. 2010;Xuan et al. 2016;Deike et al. 2017;Hao & Shen 2019). These studies directly model the evolution of both the air and wave fields in a coupled manner, no longer prescribing a fixed wave shape (e.g. Hara & Sullivan 2015;Hao et al. 2018). Furthermore, some also qualitatively consider how wave shape evolves under the influence of wind (e.g. Yan & Ma 2010;Xie 2014Xie , 2017. However, theses analyses are focused on other parameters and do not quantify precisely how the wave shape changes. Additionally, there have been a small number of field measurements (e.g. Cox & Munk 1956) and laboratory experiments (Feddersen & Veron 2005;Leykin et al. 1995) that have directly investigated how wind affects wave shape. It was found that the skewness and asymmetry depended on wind speed for mechanically-generated waves in relatively deep (Leykin et al. 1995) or intermediate and shallow (Feddersen & Veron 2005) water. In particular, the wave asymmetry (Leykin et al. 1995), skewness (Cox & Munk 1956) and energy ratio of the first harmonic (frequency 2f ) to the primary (frequency f ) (Feddersen & Veron 2005) all increased with wind speed. It would be beneficial to develop a theory that can explain these experimental findings.
In this paper, we develop a theory coupling wind to dynamically evolving intermediate and deep-water waves (kh 1 with k the wavenumber and h the water depth). We consider the fluid domain beneath a periodic, progressive wave that is forced by a prescribed, wave-dependent surface pressure profile. That is, the atmosphere is not treated dynamically. The wind effect on wave shape requires a nonlinear theory. As the surface boundary conditions for gravity waves are nonlinear, the equations are solved using a multiple scales perturbation analysis where the wave steepness ε := a 1 k (with a 1 the amplitude of the primary wave) is small and new, slower timescales are introduced over which the nonlinearities act (see, for example Mei et al. 2005). This formalism has been used to derive the canonical Stokes waves, which are periodic, progressive waves of permanent form in intermediate and deep water (Stokes 1880). By introducing a surface pressure forcing term, we will allow for solutions of the form with the real part is implied. Here η is the wave height, ω is the wave frequency, and a 1 k = ε and a 2 k = O ε 2 are the nondimensional amplitudes of the primary and first harmonic, respectively. We have defined a new parameter, the "harmonic phase" (or HP) β, which is analogous to the biphase, a statistical tool (Elgar & Guza 1985). The wave shape (i.e., skewness and asymmetry) are functions of HP β and relative harmonic amplitude a 2 /(a 2 1 k). For example, both skewness and asymmetry are zero for a 2 /(a 2 1 k) = 0, and asymmetry is zero for β = 0. For a deep-water (kh ≫ 1) Stokes wave without wind forcing, ω = g/k 1 + 1 2 a 2 1 with g the acceleration due to gravity, a 2 /(a 2 1 k) = 1 2 giving a non-zero skewness, and β = 0 yielding no phase difference between the harmonics in (1.1). Indeed, unforced Stokes waves are exactly symmetric at all orders (Toland 2000).
Three surface pressure profiles, derived from the theories of Jeffreys (1925) and Miles (1957), are included in the perturbation expansion. Using the method of multiple-scales, Stokes wave-like solutions are derived, giving the wave shape (via a 2 /(a 2 1 k) and β) dependence on the wind-induced surface pressure profile. Additionally, wave growth will result from the fact that Im{ω} is no longer zero (e.g. Miles 1957). These solutions reduce to the unforced Stokes waves when the pressure forcing vanishes.
In section 2, we set up the equations and define the different pressure profiles used. Section 3 begins the general derivation covering a range of realistic pressure magnitudes, which is continued in appendix A. A key aspect to the derivation (section 3 and appendix A), we include the nondimensional pressure p ′ in the leading-order equations (p ′ = O(ε)) which is the most nuanced, but also allows the formal replacement p ′ → εp ′ or p ′ → ε 2 p ′ , generating the weaker p ′ = O ε 2 and p ′ = O ε 3 solutions (appendix A.5). In section 4, we discuss the results of this analysis. Section 5 details the relevance of these findings and connects them to existing studies. Appendix A extends the general derivation to higher orders in ε to demonstrate a weak amplitude dependence of the shape parameters.

Governing Equations
Here, we specify the equations governing the water dynamics, Homogeneous, incompressible fluids satisfy the incompressible continuity equation within the fluid. We assume irrotational flow, allowing the water velocity u to be written in terms of a velocity potential as u = ∇φ. We define a coordinate system with z = 0 at the mean water level, positive z upwards, and gravity pointing in the −z direction. We assume planar wave propagation in the +x direction and uniform in the y direction. Then, the incompressibility condition becomes Laplace's equation: (2.1) Assuming water of uniform depth with a flat bottom located at z = −h, we impose a no-flow bottom boundary condition Finally, the standard surface boundary conditions (e.g. Whitham 2011) are the kinematic boundary condition, with η(x, t) the surface profile and p(x, t) the surface pressure evaluated at z = η. In section 2.3 we specify the surface pressure profiles.

Assumptions
Our analysis is characterized by a number of nondimensional parameters. The wave slope ε := a 1 k, assumed small, will order our perturbation expansion. Additionally, we will restrict our attention to intermediate and deep water by requiring that the nondimensional depth kh O(1) so that the Ursell parameter is small, ε/(kh) 3 ≪ 1. An additional parameter is the nondimensional surface pressure magnitude induced by the wind discussed in sections 2.3 and 2.4. We seek waves with wavelength λ := 2π/k travelling in the x direction that are periodic in x and quasi-periodic in t: (2.5) with θ defined in (A 66) and analogous to the standard kx − Re{ω}t. Additionally, we neglect the surface tension σ by restricting to wavelengths λ ≫ 2 cm implying a large Bond number (ρg/k 2 σ ≫ 1). Furthermore, we assume no mean Eulerian current. Finally, we seek only primary waves and their bound harmonics as solutions.
In the dynamic boundary condition (2.4), we incorporated the normal stress (surface pressure) but neglected the shear stress, as the normal stress is usually significantly larger than the shear stress (e.g. Kendall 1970;Hara & Sullivan 2015;Husain et al. 2019). Additionally, we note that surface shear stresses cause a slight thickening of the boundary layer, which is equivalent to a pressure phase shift on the remainder of the water column (Longuet-Higgins 1969). Therefore, we can include the effect of shear stresses through a phase shift in the pressure relative to the wave profile. Hence, in this investigation we only consider pressures acting normal to the wave surface. The irrotational assumption was motivated by our assumption that vorticitygenerating wind shear stresses are small; additionally, any such vorticity is constrained to a thin boundary layer just below the surface of the waves (Longuet-Higgins 1969). Finally, viscous forces vanish-necessary for Bernoulli's equation (2.4)-for any flow that is both irrotational and incompressible (with constant viscosity; e.g. Fang 2019). Thus, we will assume irrotational, inviscid flow throughout the fluid interior.

Surface Pressure Profiles
Here, we define the surface pressure profiles used in the analysis. The Jeffreys (1925) theory yields a ("Jeffreys") surface pressure profile, with η(x, t) the surface profile, ρ a the air density, U a characteristic wind speed, and s an empirical, unitless, order-1 constant. In contrast, the Miles (1957) theory of wind-wave growth gives a ("Miles") surface pressure profile of the form withα andβ empirical, unitless, order-1 constants, k the wavenumber, and η a the analytic representation † of η.
Miles's theory was developed for a single sinusoidal wave profile without harmonics. Although naively applying (2.7) to a multi-harmonic (e.g. Stokes) wave is permitted, another suitable generalization, which captures the motivation behind the Miles profile, is specifying the wind as leading or trailing η. Thus, we define another ("Generalized Miles") surface pressure profile as p G (x, t) = rρ a U 2 kη(kx + ψ P , t) , (2.8) with r a new, unitless, order-1 constant, and the new parameter, the "wind phase" ψ P , corresponds to the phase shift between the wave and the pressure profile. As the surface pressure is elevated on the windward (relative to leeward) side of the wave, ψ P > 0 corresponds to wind blowing from the left, assuming ψ P ∈ (−π, π).
To facilitate comparison, the various pressure profiles are written in a common form. Inspired by similarities in (2.6) to (2.8), we define a non-negative pressure magnitude constant, P , which implicitly encodes the wind speed: (2.9) Our analysis will be agnostic regarding the U dependence. The form specified in (2.9) only serves as motivation. Since k η 2 = O(ε) ( · 2 is the L 2 norm), we see can from † The analytic representation of a real function f (x) is f (x) + if (x) withf (x) the Hilbert transform of f (x). For our purposes, only two representations will be relevant: the analytic representation of cos(x) is e ix and that of sin(x) is −ie ix . the various definitions (2.6) to (2.8) that (2.10) For instance, we will define P J for Jeffreys the profile such that with ψ P = ± 1 2 π. Likewise, we will rewrite the Miles and Generalized Miles profiles as p M (x, t) = P M e iψP kη a (x, t) .
(2.12) p G (x, t) = P G kη(kx + ψ P , t) . (2.13) These three surface (z = η) pressure profiles, (2.11) to (2.13), are expanded in a Fourier series to yield simpler equations. Expanding arbitrary functions f (x) in Fourier series as the real part of (2.16) Therefore, we will generically writê withP m = mP J e iψP and ψ P = ± 1 2 π for Jeffreys,P m = P M e i sgn mψP for Miles, and P m = P G e imψP for Generalized Miles profiles. Note, the wind phase ψ P is a free parameter for the Miles and Generalized Miles profiles. Thus, by appropriate choice ofP m , we can represent surface pressures as convolutions of η(x, t) with arbitrary (time-independent) functions of x.
To make these definitions concrete, and to contrast the different forcing types, a deepwater, second-order Stokes wave profile kη = ε cos θ + 1 2 ε 2 cos 2θ (2.18) is shown for ε = 0.2 ( fig. 1a). The Stokes wave profile is used to express the three (unity normalized) surface pressure profiles ( fig. 1b), where ψ P = π/2 is chosen for Miles and Generalized Miles profiles. These three pressure profiles, (2.11) to (2.13), are largely similar to each other, although differences arise due to the Stokes wave harmonics. The derivative in the Jeffreys profile (blue fig. 1b) multiplies each Fourier component by its wavenumber, mk, enhancing higher frequencies. The wind phase ψ P , measured left from θ = 0 to the pressure maximum, shifts the entire pressure waveform relative to the surface waveform η for the Generalized Miles profile (orange fig. 1b). Contrast this to the Miles profile (green fig. 1b), where ψ P shifts the phase of the m-th Fourier component by ψ P /m radians, distorting the waveform η. The LES numerical simulations of Hara & Sullivan (2015) and Husain et al. (2019) show ψ P ≈ 3π/4 for a variety of wind speeds (section 5.1). However, in fig. 1, each surface pressure profile is depicted for ψ P = 1 2 π to facilitate comparison with the Jeffreys case (for which ψ P = ± 1 2 π).

Determination of Pressure Magnitude P
We use existing experimental data to determine the magnitude of P in various contexts. Assuming a logarithmic wind profile, Miles (1957) derived the wave-energy growth rate Normalized surface pressure profiles p(θ) as described in (2.11) to (2.13); see legend. The maximum pressure magnitude is normalized to unity (arbitrary units), and a value of ψ P = 1 2 π was chosen to facilitate comparison with the Jeffreys profile with ψ P positive corresponding to wind blowing to the right. γ, normalized by the (unforced, linear, deep-water) wave frequency f ∞ 0 , for the pressure profile p M (2.12): where, c ∞ 0 = g/k is the linear, deep-water phase speed, ρ w is the water density, and (2.9) is used to define P M . Using the value ψ P = 3π/4 from Hara & Sullivan (2015) and Husain et al. (2019) gives P M k/(ρ w g) = 0.23(γ/f ∞ 0 ). Furthermore, we use empirical data relating wind speed U to growth rate to constrain the P M pressure magnitude constant in deep-water. Figure 2 shows the energy growth rate γ/f ∞ 0 as a function of inverse wave age, u * /c ∞ 0 with u * the friction velocity. The empirical observations of deep-water γ/f ∞ 0 versus u * /c ∞ 0 collapse onto a curve permitting a conversion from u * /c ∞ 0 to γ/f ∞ 0 , yielding P M k/(ρ w g) (2.19). Here, we consider p 2 k/(ρ w g) = O(ε) to O ε 3 , or P k/(ρ w g) = O(1) to O ε 2 . If we assume ε ≈ 0.1 and ψ P ≈ 3π/4, then (2.19) shows we are considering growth rates γ/f ∞ 0 ≈ 4 × 10 −2 to 4. Referring to fig. 2, we see these reside mostly in the laboratory measurement regime, corresponding to u * /c ∞ 0 ≈ 5 × 10 −1 to 5. We can approximate U 10 using the standard logarithmic boundary layer theory (e.g. Monin & Obukhov 1954):  Figure 2. Nondimensional, deep-water wave-energy growth rate γ/f ∞ 0 versus inverse wave age, u * /c ∞ 0 with u * the wind's friction velocity and c ∞ 0 = g/k the linear, deep-water phase speed. The filled symbols represent laboratory measurements while the hollow symbols represent field measurements (from Komen et al. 1994). The solid line represents the fit parameterized by Banner & Song (2002).
It is interesting to examine the pressure forcing magnitude used previously. Phillips (1957) modelled wave growth using a different mechanism than considered here where the pressure forcing was included at the same order as η. That is, p 2 k/(ρ w g) = O(ε), or P k/(ρ w g) = O(1). Referring to fig. 2, this corresponded to strongly forced waves and fast wind (u * /c ∞ 0 = O(1)). Other theoretical works have often chosen p 2 k/(ρ w g) = O ε 2 (e.g. Janssen 1982; Kharif et al. 2010;Leblanc 2007;Onorato & Proment 2012). Thus, the choices of P k/(ρ w g) = O(1) to O ε 2 are all relevant in the literature.

Nondimensionalization
Nondimensional systems are useful in perturbation expansions. Here, a standard nondimensionalization (e.g. Mei et al. 2005) is performed, defining new nondimensional, order-1 primed variables: Notice the factor of ε in the equations for η and φ, since these are assumed small. Unlike in the standard Stokes wave problem, the surface pressure must also be nondimensionalized. As show in (2.10), O( p 2 ) = εO(P ). Thus, we find p and P (as well as their Fourier transforms) are nondimensionalized by with p ′ (x, t) and P ′ (as well as their Fourier transforms) now order-unity and dimensionless. For the remainder of the paper, primes will be dropped and all variables will be assumed nondimensional and order-unity, except where explicitly stated.

Derivation of Wave Shape Parameters
We now couple a prescribed surface pressure profile (2.17) to the nonlinear wave problem (2.1) to (2.4) to derive the wind's effect on wave shape. The end result will be an expression for the nondimensional surface profile of the form Both the HP β and a 2 /a 2 1 encode information about the wave shape. We take the ratio a 2 /a 2 1 because we will find that a 2 ∝ exp(2 Im{ω 0 }t 0 ) while a 1 ∝ exp(Im{ω 0 }t 0 ). As we are mainly interested in the shape, the growth is removed by using the ratio a 2 /a 2 1 . Expanding our nondimensional variables in an asymptotic series of ε, we have Recall that the nondimensional pressure forcing can enter at different orders depending on the wind strength (cf. section 2.4). Here and in appendix A, we include pressure in the leading-order equations (P k/(ρ w g) = O(1), i.e. p 1 = 0) which is the most nuanced, but also allows the formal replacement P → εP or P → ε 2 P generating P k/(ρ w g) = O(ε) and P k/(ρ w g) = O ε 2 solutions (appendix A.5). Here, we solve the equations to O ε 2 , giving the leading order contribution to our shape parameters, β and a 2 /(a 2 1 k). In appendix A, we solve the equations to O ε 4 , yielding the next non-zero correction to the shape parameters (demonstrating a weak amplitude dependence).
Laplace's equation (2.1) is solved via a Fourier transform and, with the bottom boundary condition (2.2), has solution (real part implied) with arbitrary m ∈ N and arbitrary functionφ n,m (t 0 , t 1 , . . .). Furthermore, to express the surface pressure profile p n in terms of the surface height η n (cf. (2.17)), all variables are written as a Fourier series: Aside from the pressure expansion, this follows the standard Stokes expansion methodology (e.g. Mei et al. 2005;Ablowitz 2011). Recall that we previously related (cf. (2.17)) the Fourier-transform of the surface pressure to the surface profile, Thus, the pressure has higher-order corrections because η has higher-order Stokes-like corrections.
We now expand the kinematic (2.3) and dynamic (2.4) boundary conditions in ε and collect terms order-by-order. O(ε) : (3.11) (3.14) We solve these equations to O ε 2 here and leave the solution of O ε 4 equations to appendix A.

O(ε) Equations
Here, we assume the pressure enters the leading order equations with P k/(ρ w g) = O(1). Proceeding to first-order in ε, the linearised boundary conditions are Inserting the Fourier transforms (3.3) to (3.5) and substituting the pressure profile (2.17) givesφ Combining equations to eliminateη 1,1 gives This is the usual, finite depth, linear operator onφ 1,1 modified by the presence ofP 1 , showing thatφ 1,1 (t 0 , t 1 , . . .) is harmonic. Using a bit of foresight to define the constants, we writeφ where The (+) sign corresponds to waves propagating to the right. Inserting this into the surface boundary conditions gives equations for η 1 , This gives with ω 0 given above. It is instructive to consider the real and imaginary parts of ω 0 :
The surface boundary conditions give equations for η 2 : These have the solution Redimensionalizing, we find where we have defined C 2,2 as (3.40) Note that A 2 1 |C 2,2 |/k is the quantity we called A 2 in (3.1). We have now found the primaryη m=1 = εη 1,m + O ε 3 and first harmonicη m=2 = ε 2η 2,2 + O ε 3 . Therefore, the amplitudes of the primary and harmonics are respectively Hence, in order to cancel the t 0 -dependence, we define the relative harmonic amplitude shape parameter as With this definition, (3.39) becomes where we have defined the harmonic phase β as the complex angle ofη m=2 /η 2 m=1 : (3.45) In general, both β and a 2 /(a 2 1 k) will have an expansion in ε sinceη m=2 will have higherorder corrections. For instance, the HP β has expansion β = β 0 + εβ 1 + . . .. Inserting our solution (3.39) into (3.45) gives β 0 , which is just the complex angle of C 2,2 at this order: with a asterisk representing the complex conjugate. Similarly, using (3.43) shows that the leading order term of a 2 /(a 2 1 k) is just |C 2,2 | to this order: (3.47) Without wind (P 1 =P 2 = 0), C 2,2 reduces to 1 4 (2+3 csch 2 (h)) coth(h), or 1 2 in deep water. This reproduces the usual, Stokes waves values of a 2 /(a 2 1 k) = 1 2 in deep water and β → 0. By solving the kinematic and dynamic boundary conditions to O ε 2 we have generated the leading-order terms for β and a 2 /(a 2 1 k). Note that these are also converted to the wave skewness and asymmetry in appendix A.3. We continue this analysis by solving to O ε 4 in appendix A, deriving the first non-trivial correction and demonstrating a weak time-and amplitude-dependence to the shape parameters.

Results
Now, we present the main results of this theory. The harmonic phase β, harmonic coefficients a 1 and a 2 , and complex frequency ω depend on the four nondimensional parameters: the wave steepness ε := a 1 k, the water depth kh, the pressure magnitude constant P k/(ρ w g), and the wind phase ψ P . To reduce the nondimensional parameter range, we keep a fixed ε = 0.1. Recall (section 2.2) the requirement of ε/(kh) 3 1, such that the expansion remains properly ordered, implies kh 0.5, though we keep kh 1. Note that taking kh to ∞ yields solutions on infinite depth. The pressure magnitude constant P is P J , P M , or P G corresponding to the choice of pressure profile. For each solution, taking P → 0 recovers the unforced Stokes wave.
For the remainder of the paper, we will revert to dimensional variables. In particular, the pressure constant P is dimensional again and not necessarily order-1. Replacing the multiple timescales with the true time t in our solution (3.44), we obtain (cf. section 3) a surface height profile η of the form (A 67) kη = (a 1 k)e iθ + (a 1 k) 2 a 2 a 2 1 k e 2iθ+β + . . . , (4.1) with the real part implied. Here, θ is defined in (A 66) and is analogous to the standard kx − Re{ω}t with ω(t) = ω 0 + εω 1 (t) + . . . the full, time-dependent frequency (A 60) and waves propagating to the right for Re{ω} > 0. On the other hand, the wind blows from the left if ψ P > 0, assuming ψ P ∈ (−π, π). Note that the growth of the harmonics means the perturbation expansion is liable to become disordered after a long time. Therefore, these solutions are only valid for finite time.

Harmonic Phase, Relative Harmonic Amplitude, and Wave Shape
The harmonic phase β quantifies the relative phase shift between the primary and first harmonic and was derived (A 52) for all pressure profiles satisfying (2.17) with magnitude P k/(ρ w g) = O(1) to O ε 2 . We now specialize these results to the three pressure profiles of interest.
it for small ε ≪ 1 as (3.47) a 2 a 2 1 k J = 2 + 3 csch 2 (kh) 4 coth(kh) 1 + P 2 k 2 /(ρ 2 w g 2 ) 1 + P 2 k 2 /(ρ 2 w g 2 ) × coth 2 (kh) 2 + coth 2 (kh) (4.8) However, weak wind O(P k/(ρ w g)) ≪ 1 moves the P -dependence to a higher-order than we calculated, simply reverting to the unforced Stokes wave value with A defined in (A 19). Here, the pressure forcing only plays an indirect role by causing a 1 to grow, with direct influence pushed to higher order corrections. In contrast, the Miles profile lacks a direct pressure influence even when P k/(ρ w g) = O(1), again reducing to the unforced Stokes wave value: In general, the relative harmonic amplitude (3.47) a 2 /(a 2 1 k) reduces (to leading order in ε for P k/(ρ w g) = O(1)) to the unforced Stokes case (P 1 =P 2 = 0) precisely when (P 2 −P 1 )/(1 +P 1 ) = 0, which is only satisfied by the Miles profileP 1 =P 2 . Thus, the Miles profile is the unique P k/(ρ w g) = O(1) profile type which has no impact on wave shape (β or a 2 /(a 2 1 k)), at least to leading order in ε. The complete Generalized Miles a 2 /(a 2 1 k) G (A 51) is also plotted in figs. 3b, 4b and 5b. We can simplify a 2 /(a 2 1 k) G by assuming a small wave-steepness ε ≪ 1: (4.11) Instead assuming O(P k/(ρ w g)) ≪ 1 gives instead a 2 a 2 1 k G = 2 + 3 csch 2 (kh) 4 coth(kh) 1 + P k ρ w g sin ψ P (2 cosh ψ P − 1) (4.12) For the chosen values of kh = ∞, ε = 0.1, and P k/(ρ w g) = 1, both the Jeffreys and Generalized Miles profiles induce a harmonic phase magnitude |β| up to π/4 ( fig. 3a). The Jeffreys value of β J = π/4 is placed at ψ P = 1 2 π to correspond with its restriction that ψ P = ± 1 2 π. The Generalized Miles HP β increases from zero at ψ P = 0 ( fig. 3a) to roughly π/16 before decreasing to approximately −π/4 and passing through zero near ψ P = 1 2 π ( fig. 3b). The angle ψ P = 3π/4 is denoted by a dashed line in fig. 3, and this ψ P is utilized hereafter, as suggested by Hara & Sullivan (2015) and Husain et al. (2019). In contrast to the harmonic phase, the relative harmonic amplitude shows opposing behaviour for the two forcing types. The Jeffreys a 2 /(a 2 1 k) J = 0.7 is enhanced relative to the deep-water Stokes value a 2 /(a 2 1 k) = 1 2 , while the Generalized Miles value is suppressed a 2 (a 2 1 k) G 1 2 . Note that fig. 3 only depicts ψ P 0 since β (3.46) is antisymmetric and a 2 /(a 2 1 k) (3.47) is symmetric about ψ P = 0. This is seen by noticing ψ P → −ψ P =⇒P m →P * m . The wave shape parameters show a particularly rich dependence on the pressure magnitude P k/(ρ w g) ( fig. 4). While both Jeffreys and Generalized Miles yield nonzero harmonic phase β for small pressures ( fig. 4a), they have qualitatively different behaviours for large P k/(ρ w g). The Jeffreys profile increase steadily, reaching 3π/8 for P k/(ρ w g) = 3. Instead, the Generalized Miles profile first decreases, reaching a minimum of approximately − 1 4 π at P k/(ρ w g) = 0.6 and then increasing to small, positive values. The relative harmonic amplitude shows ( fig. 4b) virtually no change from the deep-water Stokes value of 1 2 until P k/(ρ w g) = 0.3. Then, the Jeffreys profile increases rapidly, attaining a 2 /(a 2 1 k) J = 1.7 for P k/(ρ w g) = 3. Contrarily, the Generalized Miles profile decreases and asymptotes to a 2 /(a 2 1 k) G ≈ 0.2. Finally, the nondimensional depth kh also modulates the wind effect on wave shape. For the chosen values of P k/(ρ w g) = 1 and ψ P = 3π/4, the Generalized Miles β G ≈ − 1 4 π while Jeffreys β J ≈ + 1 4 π for large kh ( fig. 5a). However, as the kh decreases, both values grow in magnitude with β J increasing faster, nearly reaching β J = 1 2 π at kh = 1. Thus, the shallow depth kh strongly enhances the effect of wind on β. The wind's influence on a 2 /(a 2 1 k) is less pronounced. Notice that the unforced Stokes wave also has a depth dependence for a 2 /(a 2 1 k) (dashed line in fig. 5b). Though the relative harmonic amplitude is enhanced for small kh in all three cases (Jeffreys, Generalized, and unforced Stokes), notice that both pressure profiles grow slower than the unforced Stokes wave. That is, the pressure forcing appears to counteract shoaling-induced a 2 /(a 2 1 k) enhancement to some extent. Figure 5b also highlights the importance of restricting kh 1. As kh decreases, a 2 becomes large compared to a 1 and the perturbation expansion could become disordered.
Both the harmonic phase and the relative harmonic amplitude determine the wave shape. We consider their combined influence by plotting the surface profile under the action of the Jeffreys pressure profile. Figure 6a shows how the surface profile η versus phase θ varies with P k/(ρ w g) = 0, 0.1, and 1.0 for wind blowing to the right. The P k/(ρ w g) = 0 profile has skewness (A 56) S = 0.3 and asymmetry (A 57) A = 0 as expected for a kh = 1 Stokes wave. The P k/(ρ w g) = 0.1 deviates only slightly from the unforced profile. However, the P k/(ρ w g) = 1 profile shows a noticeable horizontal asymmetry, with both skewness S = 0.05 and asymmetry A = 0.3 that are fundamentally different from a Stokes wave. This follows from fig. 4a: P k/(ρ w g) = 0.1 generates a very small β J , while β J is significantly larger for P k/(ρ w g) = 1. We can also see that increasing depth kh decreases the influence of wind on asymmetry ( fig. 6b). The kh = ∞ profile (S = A = 0.1) is less asymmetric than the kh = 1 profile, in agreement with fig. 5. Note that both panels of fig. 6 show the wave tilted towards the wind direction, a result of the positive β J for Jeffreys type forcing ( fig. 3a).

Phase Speed and Growth Rate
In addition to influencing wave shape, the pressure forcing terms also affect the phase speed, as predicted by Jeffreys (1925) and Miles (1957). We normalize the phase speed c = Re{ω}/k by the unforced, linear phase speed c 0 = g tanh(kh)/k. The complete fractional phase speed change ∆c/c 0 is given in (A 61). If we consider small waves ε ≪ 1, then (A 61) simplifies considerably: (4.13) with ψ P = ± 1 2 π for the Jeffreys profile. If we instead assume the forcing is weak, O(P k/(ρ w g)) ≪ 1, we find For these limiting cases, we find that all three surface pressure profiles generate the same change to the phase speed; this is unsurprising since, at leading order, all three pressure profiles are equivalent (if ψ P = ± 1 2 π). The a 2 1 term is the amplitude-dispersion due to nonlinearity described by Stokes (1880).
Notice that, for both the Jeffreys and Generalized Miles profiles, the HP β and growth rate are related for small waves (ε ≪ 1) and weak wind (O(P k/(ρ w g)) ≪ 1) as (4.17) The connection with wave asymmetry (related to β) suggests a deeper link between wave growth and wave shape. This is potentially analogous to shoaling, weakly nonlinear waves that both grow and becomes asymmetric.

Discussion
We have shown that the wind, via a surface pressure profile expressed as a timeindependent η convolution, affects wave shape, in addition to the previously known changes in phase speed and growth rate. The different surface profiles (Jeffreys, Miles, and Generalized Miles) produce qualitatively different results. The Miles profile has no effect on wave shape, and the Generalized Miles (for ψ P = 3π/4) HP has the opposite sign as the Jeffreys HP. Here, we will use results from large eddy simulations (section 5.1) to provide guidance on the choice of ψ P and P k/(ρ w g), allowing comparison with a laboratory experiment (section 5.2). We then discuss the surface pressure profiles in the context of laboratory observations and LES as well as potential future directions (section 5.3).  (fig. 7). The nondimensional surface perturbation pressure p ′ H varies over a range from ±20 with maximum shifted ≈ 3π/4 windward of the wave crest (red bar in fig. 7), yielding our choice of ψ P ≈ 3π/4. The Husain et al.

Comparison of Theory to Laboratory Wave Shape Observations
Here, we compare our predicted harmonic phase to the laboratory experiments in Leykin et al. (1995). We cannot compare to Feddersen & Veron (2005) as their kh 1.2, and the u * /c 0 to γ/f 0 relationship ( fig. 2) needed for determining P k/(ρ w g) is for deep water. In Leykin et al. (1995), laboratory wind-generated surface gravity waves with ε ≈ 0.15 and kh = 2.5 had a quasi-linear relationship between the biphase β at the peak frequency (the statistical analogue of our harmonic phase β) and the inverse wave age u * /c 0 ( fig. 8). For comparison, our pressure magnitude P k/(ρ w g) must be converted to an inverse wave age u * /c 0 (section 2.4), where we assume the deep-water relationship between u * /c 0 and γ/f 0 ( fig. 2) holds for kh = 2.5, which is parameterized (Banner & Song 2002) as ( fig. 2, solid line) Using (2.19), we can relate γ/f 0 to P k/(ρ w g) for deep-water to give P k ρ w g = 32.5 sin ψ P ρ a ρ w u * c 0 2 (5.4) allowing comparison between theory and laboratory observations. Using (5.4), the measured inverse wave ages u * /c 0 = 0.5 to 1.5 correspond to pressure magnitudes P k/(ρ w g) = 0.01 to 0.1. Assuming a Generalized Miles pressure profile with ψ P = 3π/4, the predicted and measured β are in qualitative agreement (compare red curve to symbols in fig. 8). We emphasize that the conversion between u * /c 0 and P k/(ρ w g) (5.4) is only approximate. If the conversion factor were a factor of 3 larger, the results would match reasonably well. We also note that the relatively high wind speeds (u * up to 1.7 m s −1 ) likely caused additional physical processes, such as whitecapping or microbreaking, to occur. Such dissipative processes are not considered in our theoretical treatment.

The Surface Pressure Profile
Most theoretical treatments of wind-induced wave growth utilize a linear theory with monochromatic waves (e.g. Miles 1957;Belcher & Hunt 1993;Young & Wolfe 2014). In this scenario, for the same ψ P , the pressure profiles considered are identical at leadingorder and one need not distinguish between, for instance, the Miles or Generalized Miles  Leykin et al. (1995) laboratory experiments. The black, dashed line is the Leykin et al. (1995) linear fit. Theoretical HP β (solid red) are given for the Generalized Miles pressure profile with ψ P = 3π/4, kh = 2.5, and ε = 0.15, and conversion of u * /c 0 to P k/(ρ w g) is given by (5.4) (cf. section 5.2).
profiles. However, when considering higher-order corrections to the higher harmonics, differences arise and care must be taken when choosing the pressure profile.
Direct measurements of the surface pressure profile are challenging and rare (Donelan et al. 2006). However, our theory can offer insight by comparing the profiles' differing effects on wave shape parameters to simulations and measurements of wind-forced waves. We have shown that the Miles pressure profile (2.12) generates no change to wave shape, in particular HP β = 0. This disagrees with the experimental measurements that find a nonzero β (Feddersen & Veron 2005;Leykin et al. 1995). Thus, the Miles profile appears to be an inappropriate pressure profile. Both Feddersen & Veron (2005) and Leykin et al. (1995) measure a harmonic phase β < 0 for co-aligned wind and waves. However, the Jeffreys profile gives a positive β while the Generalized Miles profile with ψ P ≈ 3π/4 gives a negative β (figs. 3a and 4a). Additionally, the Jeffreys requirement of ψ P = ± 1 2 π appears inconsistent with numerical simulations (Hara & Sullivan 2015;Husain et al. 2019). This suggests that the Generalized Miles case is the more appropriate surface pressure profile of those considered here.
Throughout the derivation, we have maintained a somewhat general surface pressure profile p(x, t), namely any time-independent convolution with η (cf. section 2.3). However, LES atmospheric simulations over a simple sinusoidal (i.e. no harmonics) topography (e.g. Hara & Sullivan 2015;Husain et al. 2019) result in surface pressure profiles that are not simply sinusoidal (i.e. they have skewness or asymmetry, cf. fig. 7), counter to the assumption thatp m ∝η m . Our small ε theory could be extended to allow pressures with Fourier representationsp m = kP mηm + k 2 nP m,nηnηm which could better represent the LES simulated surface pressure of Husain et al. (2019). Additional surface pressure complexity is likely generated if LES atmospheric simulations used a Stokes wave bottom profile instead of a single sinusoid. Allowing the wind, via surface pressure profiles, to affect wave shape, as done here, likely induces further changes back to the airflow and surface pressure profile. That is, the air and water phases are coupled. Although this investigation relied on prescribed surface pressures, it lays the groundwork for a weaklynonlinear coupled theory. Future work will attempt to couple the wind and waves directly, providing insight into the surface pressure profile and the related wave shape and growth.

Summary
Here, we derive a theory for the wind effect on the shape of surface gravity waves. The influence of the wind on ocean waves has been studied in great detail theoretically, numerically, and observationally with a principal focus on wave growth. A few laboratory and numerical experiments have shown that wind can influence wave shape, though no theory for it exists. Two key weakly-nonlinear wave shape parameters are the harmonic phase β, encoding the relative phase between the primary and harmonic (zero for unforced Stokes waves), and the relative harmonic amplitude a 2 /(a 2 1 k). These two parameters can also be converted to the more conventional skewness and asymmetry. Motivated by prior wind-wave generation theories, three surface pressure profiles (Jeffreys, Miles, and Generalized Miles) based on a convolution with the wave profile η are prescribed. A multiple scales perturbation analysis is performed for the small wave steepness ε := a 1 k. The deep to intermediate water theoretical solutions are derived for quasi-periodic progressive waves yielding the wind-induced changes to β and a 2 /(a 2 1 k), as well as higher-order corrections to the previously known growth and phase speed changes. These parameters are functions of the four nondimensional parameters: the wave steepness a 1 k, the depth kh, the pressure magnitude P k/(ρ w g), and the wind angle ψ P . By formal replacement of the pressure magnitude P with P → εP or P → ε 2 P , our derivation permits a variety of pressure magnitudes (i.e. wind speeds).
The Miles profile had no effect on wave shape in contrast to laboratory observations. The relative harmonic ratio a 2 /(a 2 1 k) displayed a strong dependence on the forcing type, enhanced for Jeffreys but suppressed for Generalized Miles. The harmonic phase β had more complicated behaviour, including a local minimum for the Generalized Miles case as a function of the pressure magnitude. Despite restricting our analysis to intermediate and deep water, we found decreasing kh enhances the wind effect on wave shape. This suggests pressure forcing could play a larger role in wave shape for shallow water waves. We also found direct relationships between growth rates and β for the pressure profiles considered. Atmospheric large eddy simulations constrained both the pressure magnitude P and wind angle ψ P . Using the constrained ψ P , our HP predictions were qualitatively consistent with laboratory observations. Only the Generalized Miles profile could reproduce the observed sign for β, suggesting that Generalized Miles surface pressure profiles best represent the actual wave surface pressure profile. Future studies will investigate the shallow water limit. Other avenues for future work include dynamically coupling the air to the wave field. Such an approach would obviate the need to impose a specified pressure profile, increasing the applicability of the theory.

A.1. O ε 3 Equations
In section 3, we derived the leading order contributions to the HP β and relative amplitude a 2 /(a 2 1 k). Now, we will extend this derivation to the next non-zero correction. This will reveal a weak amplitude-and time-dependence to these shape parameters. Furthermore, by finding β and a 2 /(a 2 1 k) accurate to O ε 2 , we can formally take P → εP yielding solutions with P k/(ρ w g) = O(ε), or P → ε 2 generating P k/(ρ w g) = O ε 2 results (appendix A.5). However, the expressions begin to become unwieldy; therefore, we will only sketch the derivation. The third-order equations give with c.c. denoting the complex conjugate. Here, KIN 3,1 , KIN 3,3 , DYN 3,1 , DYN 3,3 ∈ C are constants that do not depend on A 1 , x, t n , or z (these dependencies have been explicitly factored out) and are composed entirely of known quantities from previous orders. In general, KIN n,m and DYN n,m are the constants (depending on h, ψ P , andP m only) for the n-th order, m-th Fourier component (i.e., exp(imkx)) term from the kinematic or dynamic boundary condition, respectively. Once again, inserting our Fourier transforms (3.3) to (3.5), we find m = 1 Fourier Component: (A 4) m = 3 Fourier Component: Notice that we did not evaluate the ∂/∂t 0 derivative the ( ∂/∂t 0 − iω 0 ) of (A 7); we will discuss this momentarily.
Preventing secular terms requires that coefficients of exp(−iω 0 t 0 ) for m = 1 vanish. Thus, we require Here, we encounter an issue: given that A 1 (t 2 , t 3 , . . .) is explicitly not a function of t 0 , there is no (nontrivial) way to satisfy the t 0 dependence of this compatibility condition. This could be dealt with rigorously by allowing the fast timescale t 0 to modulate the slower timescales when we defined our multiple scales expansion. We encounter this issue because the growth on the fast timescale affects the period of the slower timescales. This could be dealt with formally if we had instead defined our multiple scale expansion with additional, fast-timescale dependences: with the primes to make our new timescales distinct from the originally defined ones. Then, we can choose the form of µ n to remove secular terms. This modified multiple scales approach is similar to the one specified in Pedersen (2006). Using this freedom to remove these problematic secularities, we would find that This method would eliminate the need to be careful about the ( ∂/∂t ′ 0 − iω 0 ) ∂A 1 /∂t ′ 2 term previously mentioned, and would eliminate the e 2 Im{ω0}t ′ 0 term we are attempting to deal with currently. Later, to re-express the solution in terms of t, a simple integration yields where we required that t ′ n = 0 at t = 0. Note: t ′ 0 is not a special case; treating n as a continuous variable and taking the limit n → 0 recovers t ′ 0 = t. Note that, since our previous solutions had no t 1 dependence, making this post hoc change to t 2 does not alter any of our previous conclusions. Furthermore, we will see that only the even timescales (t 2 , t 4 , etc.) need this treatment; since we are only considering timescales up to t 3 , we will only make this replacement for t 2 .

(A 19)
With the compatibility condition solved, the m = 1 equation reduces to the homogeneous equation; for simplicity, we will chooseη 3,1 = 0 . Substituting (A 14) and our solution forη 3,1 into the surface boundary conditions allows us to solve forφ 3,1 . Assuming a solution of the form The second harmonic (m = 3) equation is solved forη 3,3 as usual. Then, substituting this solution into the surface boundary conditions permits solving forφ 3,3 .
Thus, we have the solutions Finally, going to fourth-order, we have Here, KIN 4,0 , KIN 4,2 , KIN 4,4 , DYN 4,0 , DYN 4,2 , DYN 4,4 ∈ C are constants that do not depend on A 1 , x, t n , or z (these dependencies have been explicitly factored out) and are composed entirely of known quantities from previous orders. Inserting the Fourier transforms (3.3) to (3.5) gives m = 2 Fourier Component: Preventing secular terms requires that ∂ t3 A 1 = 0. These can be solved as usual forφ 4,m ; using the surface boundary conditions, the solutions forη 4,m can then be determined as well.
The only terms worth discussing are the zero-modes,φ 4,0 andη 4,0 . Whileη 4,0 has physical meaning (this is a component of the setup or setdown),φ 4,0 has a gauge freedom: we may add a constant term (in x and t 0 ), as well as a term proportional to t 0 , without affecting any observables. Using this freedom, we will choose these two free constants such that theη 4,0 → 0 andφ 4,0 → 0 as P → 0.
The solutions at this order are  Here,Ã 1 := A 1 P =0 is the additive "constant" we were permitted from the m = 0 equation; note:Ã 1 could still be a function of slower timescales t 1 , t ′ 2 , etc. As mentioned previously, a term, linear in t 0 , was included inφ 4,0 ; this was necessary to include thẽ A 1 term inη 4,0 . For reference, the full solution for η is with A 1 (t 2 ) given by (A 17).

A.3. Shape Parameters
Now, we can calculate the shape parameters when pressure enters at leading order. Recall that we are seeking two parameters-the HP β, and the relative harmonic amplitude, a 2 /(a 2 1 k) (with a 2 the amplitude of the complete first harmonic, and a 1 the amplitude of the complete primary).
The primary wave is simply with A 1 (t ′ 2 ) given by (A 17).

(A 51)
We can see that the O ε 2 correction grows as a function of the fast timescale, t 0 , as well as the slow timescale, t ′ 2 (through its A 1 (t 2 ) dependence). Likewise, the HP β is the complex angle (3.45) of (A 50): Im C 2,2 + ε 2 |A 1 | 2 e 2 Im{ω0}t0 C 4,2 Re C 2,2 + ε 2 |A 1 | 2 e 2 Im{ω0}t0 C 4,2 with β 0 given in (3.46) by with a asterisk representing the complex conjugate. Notice that β now begins to show a weak amplitude, |A 1 |, dependence. Asymmetry and skewness are more common shape parameters than β and a 2 /(a 2 1 k). Therefore, we derive the asymmetry and skewness of our solution. The skewness S and