Boundary layer formulations in orthogonal curvilinear coordinates for flow over wind-generated surface waves

The development of the governing equations for fluid flow in a surface-following coordinate system is essential to investigate the fluid flow near an interface deformed by propagating waves. In this paper, the governing equations of fluid flow, including conservation of mass, momentum and energy balance, are derived in an orthogonal curvilinear coordinate system relevant to surface water waves. All equations are further decomposed to extract mean, wave-induced and turbulent components. The complete transformed equations include explicit extra geometric terms. For example, turbulent stress and production terms include the effects of coordinate curvature on the structure of fluid flow. Furthermore, the governing equations of motion were further simplified by considering the flow over periodic quasi-linear surface waves wherein the wavelength of the disturbance is large compared to the wave amplitude. The quasi-linear analysis is employed to express the boundary layer equations in the orthogonal wave-following curvilinear coordinates with the corresponding decomposed equations for the mean, wave and turbulent fields. Finally, the vorticity equations are also derived in the orthogonal curvilinear coordinates in order to express the corresponding velocity–vorticity formulations. The equations developed in this paper proved to be useful in the analysis and interpretation of experimental data of fluid flow over wind-generated surface waves. Experimental results are presented in a companion paper.

wave and turbulent equations were derived in a hybrid, wave-following coordinate system similar to the one used by Sullivan et al. (2014) and Sullivan et al. (2018). The triple decomposition formulations in wave-following coordinates define the expression for the wave-induced stress as a sum of the wave and pressure stresses compared to the traditional definition (e.g. Makin & Kudryavtsev 1999;Hara & Belcher 2004) in terms of the wave-coherent velocity components. Finally, there are some recent attempts to express the viscous tangential stresses at the air-water interface in an orthogonal wave-following coordinate system (see Buckley & Veron 2017;Iafrati, De Vita & Verzicco 2019), but we note that these estimates are valid only in the linear wave limit.
In the study of turbulent flows over propagating surface waves, the necessity of employing a wave-following coordinate system has been recognized (e.g. Hsu et al. 1981;Sullivan et al. 2000;Hara & Sullivan 2015), and, therefore, as mentioned above, the algebraic mapping of rectangular coordinates into wave-following coordinates has commonly been employed in theoretical and numerical studies, and just recently, in laboratory measurements. However, the further step of fully transforming the governing equations with the dependent flow variables into those coordinate systems is rarely taken, and to the best of the authors' knowledge, the decomposed mean, wave and turbulent equations required for the complete treatment of wave-turbulence interaction are not yet fully derived. This is in part due to the mathematical difficulties associated with the curvilinear coordinates and the fact that experimental studies have not been able to estimate the additional geometric terms appearing in the transformed equations. Therefore, in the current study we present the fully transformed governing equations of fluid flow including continuity, momentum and kinetic energy equations in the wave-following orthogonal curvilinear coordinate system along with the triple-decomposed form of those equations. The orthogonal general coordinate system is physically intuitive and appropriate for the study of the turbulent flow over surface waves since it can provide an alternative framework in which the surface-parallel continuity, momentum and energy budget equations can be thoroughly investigated, which leads to a better physical interpretation of many quantities in the governing equations. The complete transformation of the governing equations also directly allows us to account for the streamline curvature. Finally, we simplify the equations for the mean, wave and turbulent fields using classical boundary layer approximations. The fluid flow governing equations in orthogonal curvilinear coordinates are presented in § 2, and the corresponding triple-decomposed equations are expressed in § 3. We narrow our focus to the weakly nonlinear waves in § 4 and perform a boundary layer scaling. The triple-decomposed equations for the boundary layer are then derived in § 4.1. We finally offer a brief conclusion in § 5. In appendix A, we provide the velocity-vorticity formulations in orthogonal curvilinear coordinates. Analyses of experimental data using the framework developed in this paper are presented in a companion paper (Yousefi, Veron & Buckley 2020).

Fluid flow governing equations
The governing equations of fluid motion, including continuity, momentum and kinetic energy equations in a turbulent flow are first derived in an orthogonal curvilinear coordinate system. These equations are then decomposed into mean, wave-induced and turbulence components by employing the triple decomposition technique. To express the equations of motion, the flow is assumed to have a constant density, a constant kinematic viscosity and velocity components u = (u 1 , u 2 , u 3 ) in x 1 , x 2 and x 3 directions, respectively. The continuity, momentum and energy equations for an incompressible fluid can be written in an invariant vector form as ∇ · u = 0, (2.1) ρ ∂u ∂t + ρ(u · ∇)u = −∇p + ∇ · τ , (2.2) where u is the velocity vector, ρ is the density and p is the pressure which includes the static gravity term gz. Also τ = 2µS is the viscous stress tensor, µ is the dynamic viscosity, 2S = (∇u + (∇u) T ) is the strain rate tensor and e = (u · u)/2 is the kinetic energy. The velocity gradient tensor is noted as ∇u and (∇u) T denotes its transpose.
2.1. Orthogonal curvilinear coordinates Fixed rectangular coordinate systems are not able to capture the effects of streamline curvature and, particularly in the study of wind waves, no spatially or temporally averaged information can be obtained beneath the highest wave crest (e.g. Sullivan et al. 2000;Hara & Sullivan 2015). In order to investigate the wave-induced motions near the interface and/or below the wave crest, it is thus necessary to employ a coordinate system that closely follows the wave shapes. For example, Hara & Sullivan (2015) and Yang & Shen (2017) recently employed a coordinate system that follows the vertical displacement due to the surface waves. However, we find that the strictly orthogonal curvilinear coordinate system is practical and physically intuitive. In this section, we derive the continuity, momentum and energy equations in an orthogonal curvilinear coordinate system. Let x i = (x 1 , x 2 , x 3 ) represent the Cartesian coordinate system and ξ i = (ξ 1 , ξ 2 , ξ 3 ) represent a set of arbitrary orthogonal curvilinear coordinates. It is well understood that the distinction between covariant and contravariant components in general non-orthogonal curvilinear coordinates vanishes in orthogonal curvilinear coordinate systems. The orthogonal coordinate basis can then be defined as follows (Vinokur 1974;Redzic 2001;Shikhmurzaev & Sisoev 2017): where e k are the corresponding Cartesian orthonormal basis. The base vectors in orthogonal curvilinear coordinates are not necessarily unit vectors. The orthonormal base vectors in orthogonal curvilinear coordinates can be obtained by where the quantities h i are the so-called scale factors of the orthogonal curvilinear coordinate system given by K. Yousefi and F. Veron In order to express (2.1) to (2.3) in orthogonal curvilinear coordinates, we will first spell out the differential vector operators including gradient, divergence, curl and Laplacian in orthogonal coordinates. These operators are fairly standard and can be found in the literature (e.g. Aris 1962;Batchelor 1967;Anderson, Tannehill & Pletcher 1984;Redzic 2001;Nikitin 2006), but we choose to include them here for completeness. Furthermore, in order to substantially simplify the notation in the remainder of this paper, we propose a notation that builds on and modifies the standard Einstein summation convention. Indeed, throughout this work the index notation is such that no summing is carried whenever the indices are enclosed within parentheses. Accordingly, the indices within parentheses only take the value of dummy or free indexes. All other aspects of the summation convention remain unchanged; any index that is repeated twice in any term of expression is called a dummy or repeated index to be summed over the range of its values and any index not repeated is called a free index taking any value in its range. For instance, the suffix (i) in the following expression only takes the value of the free index i, i.e. u i u j κ (i)j = u 1 u 1 κ 11 + u 1 u 2 κ 12 + u 1 u 3 κ 13 , and the suffix ( j) in the following expression takes the value of the dummy index j, i.e. u j u j κ ( j)i = u 1 u 1 κ 1i + u 2 u 2 κ 2i + u 3 u 3 κ 3i . With this notation, the standard vector and tensor operators can then be expressed in the orthogonal curvilinear coordinate systems in a way that we find elegant and conveniently compact. If ϕ is an arbitrary scalar, A is an arbitrary vector and T is an arbitrary tensor field, the expressions for the gradient, divergence, curl and Laplacian operators in the orthogonal curvilinear coordinates become In the above equations, a i , b i and T ij are respectively the components of the vector A, vector B and tensor T defined based on the orthonormal basis, and κ ij are the components of the curvature matrix which account for the curvature of the coordinate system and can be defined as (2.13) Conservation equations for flow over surface waves 888 A11-7 Thus, all terms in the continuity, momentum and energy equations, (2.1)-(2.3), can now be expressed in orthogonal curvilinear coordinates. The vector operators presented in (2.7)-(2.12) are consistent with those given in the literature, see e.g. Aris (1962), Vinokur (1974), Anderson et al. (1984), Redzic (2001) and Nikitin (2006), but expressed in a compact and practical form.
2.2. Continuity, momentum and energy equations In order to express the fluid governing equations, consider a set of arbitrary orthogonal curvilinear coordinates ξ i = (ξ 1 , ξ 2 , ξ 3 ) with corresponding velocity components U = (U 1 , U 2 , U 3 ). The velocity components in orthogonal curvilinear coordinates are projections of the velocity vector into the coordinate axes and thus related to the Cartesian velocities by the coordinate transformation. The components of the velocity vector u in both Cartesian and orthogonal curvilinear coordinate systems are schematically illustrated in figure 1. Using the vector operators introduced in (2.7)-(2.12), the continuity, momentum and energy equations can be then written as FIGURE 2. Illustration of the viscous stress field along with the curvature parameter in the ξ 1 -ξ 3 plane of an arbitrary orthogonal curvilinear coordinate system. Here, r 31 = 1/κ 31 and r 13 = 1/κ 13 are the radius of curvature along the constant ξ 1 -and ξ 3 -coordinate, respectively, with respect to the ξ 1 -ξ 3 plane. The point where two tangent vectors ξ 1 and ξ 1 + dξ 1 intersect is called the centre of curvature and denoted by the point O.
Here, it is noted that the governing equations (2.14)-(2.16) are obtained using a time-independent transformation. The continuity, momentum and energy equations given here for the orthogonal coordinates are quite similar to the conventional Cartesian coordinate system, except for the additional curvature terms that, in fact, account for the curvature of the coordinate system and produced due to the spatial dependence of the base vectors. The Cartesian coordinate equations are simply recovered by noting that h i = 1 for rectangular coordinates. The reader may directly verify that the expanded form of the governing equations (2.14)-(2.16) are identical to the continuity and momentum equations given in the literature for an orthogonal curvilinear coordinate system (e.g. Raithby, Galpin & Van Doormaal 1986;Blumberg & Herring 1987;Nikitin 2006;Shen et al. 2015).
Again, the viscous stress is where the components of the strain rate tensor can now be expressed in an orthogonal coordinate system as In order to provide a more in-depth explanation of the acceleration and stress terms in (2.15), which do not traditionally appear in the Cartesian equations, we confine the following discussion to 2-D flows in which we neglect the curvature in the lateral direction. The stress components are schematically illustrated along with the curvature parameters in figure 2 in the ξ 1 -ξ 3 plane of an orthogonal, wave-following coordinate system. The radius of curvature, defined as the reciprocal of the curvature, along the constant ξ 1 -coordinate with respect to the ξ 1 -ξ 3 plane is r 31 = 1/κ 31 , whereas the radius of curvature is r 13 = 1/κ 13 along a line of constant ξ 3 . Let us rewrite the ξ 1 -momentum equation for the 2-D steady case described in figure 2 as The contributions of the curvature, or equivalently the radius of curvature, to the fluid-particle acceleration components appear as additional terms. For example, (U 1 U 3 κ 13 − U 3 U 3 κ 31 ) are due to the curvature of the coordinates in the ξ 3 and ξ 1 directions, respectively. As noted by Raithby et al. (1986), the additional component τ 13 κ 13 in the viscous stress term arises because τ 13 has a net component in the ξ 1 direction, as shown in figure 2(a). In a similar fashion, the term τ 33 κ 31 arises due to the fact that τ 33 stresses are seen to have a net component in the negative ξ 1 direction (see figure 2b). The interpretation of the momentum equations in the other directions is quite similar to the ξ 1 -momentum equation and it is not presented here to keep brevity.
In orthogonal curvilinear coordinates, both the radius of curvature and the centre of curvature are functions of the arc length and, consequently, vary from point to point. Moreover, at spatial locations where the ξ 1 -coordinate has inflection points, the curvature is zero and the radius of curvature is infinite, i.e. κ 31 = 0 and r 31 = ∞.

Application to surface waves
In the section above, we reviewed the conservation equations for a fluid in motion using an arbitrary orthogonal curvilinear coordinate system. Our interest, however, is to examine the mean, wave-induced and turbulent flow on both sides of a wavy interface. Specifically, we are interested in the airflow above surface gravity waves but do not consider waves with large curvatures such as capillary waves. Experimental results of airflow measurements over wind waves will be presented in the accompanying paper (Yousefi et al. 2020). We note here that in the derivations above, the curvilinear system does not vary in time. Strictly speaking, this simply means that the wave shape is assumed to remain unchanged as the waves propagate. For monochromatic waves, the flow is thus examined after an initial Galilean transformation in which the Cartesian system is moving with the wave phase speed so that the wave shape becomes steady. In the case of wind waves with a spectrum of wave modes, the peak wave speed can be used for the Galilean transformation (Sullivan et al. 2018). In this case, the wave shape is quasi-steady. This restricts the analysis of the flow to time scales that are faster than that over which the wave shape evolves substantially. We expect this to be the case for the airflow over wind waves when the wind speed is larger than the wave phase speed (i.e. wind forced waves). Indeed, based on the wind-wave spectrum (Elfouhaily et al. 1997;Mueller & Veron 2009), we estimate that, in the laboratory and for winds ranging from 5 to 25 m s −1 , the wave shape remains correlated at 80 % for at least 0.7 peak wave period. In the field, after 0.7 peak wave period, the correlation reduces to 70 %. Therefore, over times corresponding to a fraction of the peak wave period, the coordinate transformation outlined above (with Galilean transformation using the peak wave speed) holds for wind waves with multiple modes. We note, however, that the water-side flow is likely to evolve on time scales that are comparable or longer than the time scales over which the wave shape changes. 888 A11-10

K. Yousefi and F. Veron
The issue of separating waves and turbulence has long been a recurrent challenge in the study of surface waves (e.g. Hussain & Reynolds 1970;Lumley & Terray 1983;Thais & Magnaudet 1995). Here, in order to extract the organized wave-coherent fluctuations in the flow field from the background turbulence, we decompose instantaneous variables into a phase-averaged component, f (ξ , t), and a turbulent fluctuation component, f (ξ , t), as (Hussain & Reynolds 1970) (3.1) The phase-averaged quantity is defined as where N is the number of realizations and λ is the surface wave wavelength. The phase average is the average of the values of f at a particular phase of the wave. The phase average can be further decomposed into the sum of a mean,f (ξ ), and a wave- The mean is defined as the ensemble average over all possible phases. This separation leads to the following so-called triple decomposition of an instantaneous quantity: Here, the wave-induced motion has a zero mean but is phase-coherent with the surface waves and thus not considered turbulent per se. The general properties of the ensemble and phase averages can be found in reports by Hussain & Reynolds (1970, 1972 and . The equations of motion including continuity, momentum and energy for an organized wave in a turbulent shear flow are then derived using this triple decomposition approach.

Continuity and momentum equations
The starting point to derive the decomposed equations is to substitute the decomposed field quantities into the governing equations, and then averaging; the ensemble averaging is applied first and then the phase averaging. For an incompressible fluid, the mean, wave-induced and turbulent continuity equations in an orthogonal curvilinear coordinate system can therefore be expressed as Using continuity, the instantaneous momentum equation (2.15) can be written as Substituting the decomposed velocity and pressure fields into the momentum equation (3.6) and then applying ensemble averaging yields the momentum equation for the mean flow, where the ratio of the vertical scale to the principal radius of curvature of the surface is assumed to be small. Also, the mean material derivative defined as and wherē is the mean viscous stress,S ij is the mean strain rate tensor, −U i U j is the turbulent stress, and −Ũ iŨj is the wave-induced stress (e.g. Hsu et al. 1981;Buckley & Veron 2016). The wave-induced stress term evidently makes the mean momentum equation (3.7) different from the conventional Reynolds-averaged equations for the turbulent flows in orthogonal curvilinear coordinate systems (see, for example, Nash & Patel 1972;Richmond et al. 1986;Chen, Patel & Ju 1990). The wave-induced stress is not only significant in the exchange of momentum and energy between the wind and waves particularly in the wave generation process (e.g. Hsu et al. 1981;, but also is a substantial portion of the total stress close to the surface (e.g. Janssen 1989; Makin, Kudryavtsev & Mastenbroek 1995;Makin & Kudryavtsev 2002). Equations (3.3) and (3.7) together describe the motion of the mean field. The wave-induced momentum equation can be obtained by applying the phaseaveraged operator to the instantaneous momentum equation ((3.6) after substituting the decomposed velocity and pressure terms therein), and substracting the mean momentum equation. The wave-induced momentum equation in an orthogonal curvilinear coordinate system can therefore be expressed as is the wave-induced viscous stress,S ij is the wave-induced strain rate tensor, r ij = U i U j − U i U j is the wave-induced turbulent stress (e.g. Hsu et al. 1981;Einaudi, Finnigan & Fua 1984) which represents the oscillation of the turbulent stress due to the presence of surface waves, andR ij = U i U j −Ũ iŨj is the fluctuating part of the wave stress (e.g. Einaudi et al. 1984;Rutgersson & Sullivan 2005). Analogous to wave-induced turbulent stress,R ij can be interpreted as the nonlinear wave contribution to the total fluctuation stress. Moreover, the term U i U j in the wave-induced wave stress term describes the momentum flux due to wave fluctuations. These terms are of considerable importance in coupling the wave and turbulence fields. Equations (3.4) and (3.10) together describe the wave-induced motion.
In addition to the mean and wave-induced momentum equations, the momentum equation for the background turbulence can also be derived by subtracting (3.7) and (3.10) from the decomposed, instantaneous momentum equation. Thus, the momentum equation for the turbulence in an orthogonal general coordinate system is where is the turbulent viscous stress and S ij is the turbulent strain rate tensor. The turbulence of the fluid motion can then be fully expressed by (3.5) and (3.12).
3.2. Mean, wave-induced and turbulent kinetic energy equations In this section we expand the kinetic energy budget equations and specifically look at mean, wave-coherent and turbulent kinetic energy. Following , Finnigan & Einaudi (1981) and Einaudi & Finnigan (1993), an equation for the mean kinetic energy budget can be obtained by multiplying the mean momentum equation (3.7) by U i and then successively phase-and ensemble-averaging: 14) whereē = U i U i /2 is the mean kinetic energy per unit mass. The left-hand side of this equation describes the rate of change of the mean kinetic energy and the righthand side represents different mechanisms that precipitate such changes. The first four terms on the right-hand side of (3.14) are in flux divergence form and, consequently, describe the spatial transport or redistribution of the mean kinetic energy by the mean pressure, wave perturbations, turbulent stresses and viscous stresses, respectively. The fifth term 2νS ijSij represents the viscous dissipation of the mean kinetic energy. The last two terms in (3.14) are analogous to the well-known shear production term and represent the exchange of energy between the mean flow and the wave-coherent and turbulent fields, respectively. The mean turbulent stress is likely to be positive over the ocean surface waves (e.g. Borue, Orszag & Staroselsky 1995;Shen et al. 2003;Buckley & Veron 2016), while the mean wave stress is positive below the critical layer and negative above (e.g. Townsend 1972;Hsu et al. 1981;Sullivan et al. 2000;Yang & Shen 2010). It is also observed that the mean shear is positive over propagating surface waves (e.g. Hara & Sullivan 2015;Husain et al. 2019).
The balance of the kinetic energy for the wave-induced motion can be similarly derived by multiplying the wave-induced momentum equation (3.10) with the wave-induced velocity fields and then successively applying the phase-and ensembleaveraged operators, whereẽ = U i U i /2 is the wave kinetic energy per unit mass andē =Ũ iŨi /2 is the ensemble-averaged wave kinetic energy. As previously noted, the left-hand side of (3.15) represents the rate of change of the mean kinetic energy, and the right-hand side represents the transport, production and dissipation mechanisms producing such changes. The first four terms on the right-hand side of (3.15) express the transport of the wave kinetic energy by pressure, wave-induced turbulent stresses, wave-induced wave stresses and viscous stresses or molecular diffusion. The fifth term, similar to its counterpart in (3.14), represents the viscous dissipation rate due to the wave-induced motion. The sixth term is the production term due to the periodic wave which represents the exchanges between mean shear and wave motion. This term also appears in the mean kinetic energy budget equation (3.14) but with the sign reversed. The last two terms in (3.15) are the rate of energy transfer between the wave-induced flow and in turn wave-induced wave stress and wave-induced turbulent stress. These terms are involved in the exchange of kinetic energy between the wave and turbulent fields (see Cheung & Street 1988;Einaudi & Finnigan 1993;Rutgersson & Sullivan 2005). In order to obtain further insights, it should be noted that the third and seventh terms on the right-hand side of the wave kinetic energy equation can be combined to give 1 h which explicitly describes the transport of the wave kinetic energy by the wave fluctuating part of the total fluctuation stress (Einaudi et al. 1984). It is a redistribution term (Rutgersson & Sullivan 2005) often neglected in previous studies.

K. Yousefi and F. Veron
Finally, an equation for the balance of the turbulent kinetic energy can be obtained via multiplying the momentum equation for the background turbulence by the turbulent velocity and averaging: where e = U i U i /2 is the turbulent kinetic energy per unit mass. Similarly to the mean and wave kinetic energy budgets, the three first terms on the right-hand side of (3.17) represent the transport of the turbulent kinetic energy within the flow. They are the pressure transport, the turbulence transport and viscous diffusion terms, respectively. The fourth term represents the viscous dissipation due to the turbulent motion. The fifth term is the shear production term which describes exchanges between the mean shear and turbulence. This term appears as an energy source term in the mean kinetic energy equation (3.14), but with the opposite sign. The sixth term is the turbulencewave interaction term (e.g. Rutgersson & Sullivan 2005; Davis & Monismith 2011) describing the turbulent energy production by the waves through the action of the wave-induced turbulent stresses Cheung & Street 1988). This term is also found in the wave kinetic energy equation but with the opposite sign. It is also noted that U i U j S ij =r ijSij . The last term is the advection of turbulent kinetic energy by waves (e.g. Rutgersson (3.18) The last two terms on the right-hand side of (3.17) describe the interactions between wave and turbulent fields and only appear in a triple-decomposed turbulent kinetic energy equation. The U i U jSij ,Ũ iŨjSij andr ijSij terms are present in (3.14), (3.15) and (3.17), and clearly denote the interaction among the mean, wave and turbulent fields. Finally, in obtaining (3.14)-(3.18), we note that the doubly contracted product of a symmetric tensor with another tensor is equal to the doubly contracted product of the first tensor with the symmetric part of the second tensor (Kundu & Cohen 2002) (because the doubly contracted product of any symmetric tensor with an asymmetric tensor is zero (Boyer & Fabrie 2012)). Therefore, for example, in equation (3.14): τ turb :S =τ turb : ∇U andτ wave :S =τ wave : ∇U, whereτ turb is the turbulent stress tensor andτ wave is the wave stress tensor.

Boundary layer scaling
The governing equations of motion, developed in the preceding sections, can be considerably simplified assuming boundary-layer type flows in which the vertical length scale of the motion is small compared to the horizontal length scale. Thus, we are further considering the flow over surface waves wherein the wavelength of the disturbance is large compared to the wave amplitude, i.e. ak 1. A non-dimensionalization of the problem described by (2.14)-(2.16) is then performed to obtain the reduced form of the equations through identifying the small parameters, after which the resulting equations are made specific for surface waves. To this end, a unified scaling was developed by allowing the vertical scales to be different from the horizontal scales. Accordingly, the following set of dimensionless variables, denoted by stars, are introduced: Here a is the wave amplitude, k is the wavenumber, ak is the wave slope, and U 10 is the wind speed measured at a 10 m height. It can be noted here that time is nondimensionalized using the peak frequency as where c is the phase velocity and U 10 /c is the inverse of wave age. By substituting the dimensionless variables into the governing equations of an orthogonal curvilinear coordinate system, the resulting non-dimensional equations can be expressed, for the continuity, as and, for the momentum equations, as 888 A11-16 K. Yousefi and F. Veron (4.8) where ε = ak is the wave slope and Re = U 10 /νk is the wave Reynolds number based on the wind speed at a 10 m height. The wave Reynolds number is related to the wave roughness Reynolds numbers, first introduced by Zhao & Toba (2001), through Re = Re H 4C 1/2 d ak in which C d is the air-sea drag coefficient and Re H = u * H/ν is the wave roughness Reynolds number, where u * is the friction velocity and H = 4a is proportional to the significant wave height. A non-dimensionalized equation for the kinetic energy budget can be similarly derived. It is not presented here for the sake of brevity as it leads to a long and tedious derivation. Choosing a small vertical length scale compared to the horizontal length scale (ak 1) implies that the range of validity of the boundary layer equations is limited to kξ 3 ≈ O(ak), i.e. approximately within a wave height of the surface. Wave-coherent motions are known to penetrate the airflow up to heights on the order of the wavelength rather than the wave height; the coordinate transformations, just like wave-induced motion, do indeed show an exp (−kξ 3 ) dependency. Instead, the boundary layer scaling offered here is intended to be utilized to evaluate, for example, the scale of high order terms such as Reynolds and wave stresses compared to viscous stresses. In fact, although the viscous stress in the airflow over wind-generated waves has been the topic of many studies over the past decades (e.g. Longuet-Higgins 1969;Hsu & Hsu 1983;Banner 1990;Banner & Peirson 1998;Veron, Saxena & Misra 2007;Buckley & Veron 2016), it has just recently been shown that the contribution of the viscous stress to the total air-water momentum flux is not negligible, at least in low wind speeds (Grare et al. 2013;Buckley & Veron 2017). Therefore, in order to retain the viscous effects in the boundary layer equations, the largest viscous term is required to be of the same order of magnitude as the inertia terms, i.e. 1/Re must be of the order of magnitude of ε 2 .
In order to pursue the order of magnitude analysis, it is also necessary to determine the order of magnitude of scale factors and curvature parameters in (4.5)-(4.8) and the kinetic energy budget equation. To this end, the specifics of the coordinate transformation need to be prescribed.
Following Benjamin (1959), we adopt an orthogonal wave-following coordinate system in the frame of reference moving with the wave, i.e.
where (ξ 1 , ξ 2 , ξ 3 ) are wave-following coordinates, (x 1 , x 2 , x 3 ) are rectangular coordinates, k 1 is the wavenumber in the streamwise direction, k 2 is the wavenumber in the lateral direction, and k = |k| is the wavenumber. The actual wave profile is, to the first order in ak, given by ξ 3 = 0. This type of coordinate transformation is thoroughly used in the literature (e.g. Belcher, Newley & Hunt 1993;Sullivan et al. 2000;Náraigh, Spelt & Zaki 2011;Tseluiko & Kalliadasis 2011) as it permits a linear analysis for small-slope waves. Considering ak 1 ∼ ak 2 ∼ ak are of the same order, say ε, and retaining first-order terms, the coordinate differential can be expressed as where ϕ = (k 1 ξ 1 + k 2 ξ 2 ). Consequently, the scale factors, introduced in (2.6), are of the order of h i ≈ 1 + O(ε). (4.11) The order of magnitude of dimensionless curvature parameters (2.13) can be further estimated to where ϕ * = (ξ * 1 + ξ * 2 ). Thus, κ * i3 ∼ O(ε 2 ) and all other dimensionless curvature parameters are of order O(ε). The order of magnitude analysis can now be thoroughly established for the terms involving the scale factor and curvature parameter. The details of estimating the order of scale factors and their derivatives can be found in Tseluiko &Kalliadasis (2011) andNáraigh et al. (2011).
We note here that the transformation outlined above in (4.9) is valid for monochromatic waves. In the case of wind waves, however, multi-modal transformations are generally utilized. For example, this was done by Sullivan et al. (2014Sullivan et al. ( , 2018 for numerical simulations and by Buckley & Veron (2016) for experimental studies. Extending the transformation (4.9) to a summation of multiple Fourier modes of amplitude a n and wavenumber k n leads to h i ≈ 1 + O( a n k n ). Because of the specifics of the spectral shape for wind waves (i.e. k −5/2 to k −3 above the wave peak in the equilibrium and saturation ranges), h i ≈ 1 + O(a p k p ), where a p and k p are respectively the peak amplitude and wavenumber.
We are now in a position to complete an order of magnitude analysis for the governing equations. Neglecting all terms with the order of magnitude of O(ε 2 ) and higher yields the governing equations for the flow over wind-generated surface waves with modest slopes. Therefore, the governing equations in terms of dimensional variables are reduced to the following forms: for the continuity equation, for the momentum equations, (4.14) (4.17) In which, e = 1 2 U 1 U 1 + 1 2 U 2 U 2 is the kinetic energy. We note here that (4.14)-(4.17) reduce to the conventional boundary layer equations in the rectangular coordinates at order O(1). These equations can be compared with the viscous boundary-layer equations over a solid curved surface derived by, for example, Cebeci, Kaups & Moser (1976), Degani & Walker (1993) and Cebeci & Cousteix (2005).

Triple decomposition
The boundary layer equations given above can also be decomposed into mean, wave-induced and turbulent components by using the triple decomposition technique introduced in § 3. However, before proceeding with the triple decomposition of continuity, momentum and energy equations, we need to introduce a new set of dimensionless variables since employing the triple decomposition technique leads to the double correlation of wave and turbulent velocities including the turbulent and Conservation equations for flow over surface waves 888 A11-19 wave stresses. Experiments (e.g. Buckley & Veron 2017 indicate that, across the boundary layer, the mean velocity component in the vertical direction can be assumed to be smaller than the streamwise and lateral mean velocity components, while the wave-induced and turbulent velocities are presumably all of the same order of magnitude and much smaller than the horizontal mean velocities. The triple-decomposed equations will be then non-dimensionalized, very much in the same manner as before, but using the following dimensionless variables: The continuity equations for the mean, wave and turbulent flow fields can be readily obtained by substituting the decomposed velocities into the continuity equation (4.13) and then applying the ensemble-and phase-average operators: Deriving the decomposed momentum and energy budget equations for the boundary layer, however, requires more effort. Following substituting the decomposed fields into the momentum equations (4.14) and (4.15) and applying the ensemble and phase averaging operators, it is necessary to perform the order of magnitude analysis using the dimensionless variables introduced in (4.18)-(4.23) due to the new double-velocity correlation expressions. The details, however, are eliminated to keep the brevity. Neglecting all terms with the order of magnitude of O(ε 2 ) and higher renders the ξ 3 equations trivial. The mean momentum equations in the ξ 1 and ξ 2 directions can be expressed as (4.28) Applying the phase-averaged operator to the decomposed wave-induced momentum equations (4.14) and (4.15), subtracting the mean momentum equations, and neglecting terms of order O(ε 2 ), the wave momentum equations can be written as (4.30) In these equationsr ij = U i U j − U i U j andR ij = U i U j −Ũ iŨj . Finally, we are deriving the momentum equation for the background turbulence in the boundary layer by subtracting the mean and wave momentum equations from the decomposed momentum equations and neglecting terms of order O(ε 2 ). Consequently, the momentum equations for the turbulent field can be written in the ξ 1 and ξ 2 directions as (4.32) In addition, the mean viscous stress, which appears in (4.27), can be shown to reduce to which to the leading order is An equation for the mean, wave and turbulent kinetic energy budgets for the flow in the boundary layer can also be obtained by multiplying the mean, wave and turbulent momentum equations by the mean, wave-induced and turbulent velocities, respectively, and then successively applying the phase-and ensemble-averaging operators, is the total fluctuation stress tensor, and are the mean, wave-induced and turbulent kinetic energy, respectively. Moreover, D/Dt is the mean material derivative defined in (3.8). The reader should note that deriving the triple-decomposed equations, particularly energy budget equations, in an orthogonal curvilinear coordinate system describing the fluid motion in the boundary layer is an exercise that, while straightforward, is quite tedious. The details of the derivation are omitted here for the sake of brevity.

Validation of scale analysis
The equations for the mean, wave-induced and turbulent flows in the boundary layer adjacent to surface waves with modest slopes (equations (4.27)-(4.32)) were derived by assuming that terms O(ε 2 ) and higher are negligible. In this section we assess some of these assumptions using the experimental laboratory data of Buckley & Veron (2016 who measured 2-D velocity fields in the airflow above moving surface waves for different wind-wave conditions with wind speeds ranging from U 10 = 0.89 m s −1 to U 10 = 16.59 m s −1 . To this end, the measured velocity fields were projected and transformed into the curvilinear coordinate system, then averaged in the ξ 1 -direction. A detailed analysis of this experimental data is presented in the companion paper (Yousefi et al. 2020). For the mean momentum equation in the ξ 1 -direction (4.27), for example, the complete dimensionless equation (in two dimensions) is ∂h 2 ∂ξ * and thus excluded from the averaging process. The streamwise derivative of the horizontal wave-induced variance and, especially, the additional wave-induced stresses are almost zero everywhere except near the surface where the curvature effects are strong. They fall to zero from their maxima at a distance of nearly k p ζ = 0.075 and k p ζ = 0.05 for term I and term III, respectively.
Similarly, the streamwise divergence of the horizontal turbulent variance (term IV in (4.37)), the vertical divergence of the turbulent stress (term V in (4.37)) and the additional turbulent stresses (term VI in (4.37)) are, respectively, of the (leading) order of (4.46) The laboratory measurements of these terms are also consistent with the order of magnitude analysis; the divergence of the horizontal turbulent variance in the streamwise direction and the additional turbulent stresses are respectively one and two orders of magnitude smaller compared to the vertical divergence of the turbulent stress (not shown here to maintain brevity). Evidently, the wave and turbulent stresses that appear because of curvature in the coordinate system, are asymptotically zero and hence negligible for the flow over surface waves with a small slope.

Conclusions
In the current work, the dynamical governing equations of three-dimensional fluid motion, i.e. the continuity, momentum and kinetic energy equations, have been transformed into the orthogonal curvilinear coordinate system in such a way that flow direction and coordinate directions coincide. These equations are then separated into the mean, wave-induced and turbulent components by employing the triple decomposition technique. The complete transformation of governing equations involves curvature parameters or, equivalently, the local radius of curvature and their higher order derivatives. The transformed equations therefore involve explicit extra geometric terms, for example, the additional production, advection and diffusion terms representing the effects of streamline curvature on the structure of fluid flow including the acceleration of the mean flow. These formulations are valid down to the surface and naturally incorporate the curvature effects. Furthermore, the precise expressions for the mean, wave-induced and turbulent viscous tresses are explicitly spelled out.
We also simplified considerably the continuity, momentum and kinetic energy equations for the boundary-layer type flows based on the assumption that the vertical length scale of the motion is small compared to the horizontal length scale. Considering the flow over periodic surface waves wherein the wavelength of the disturbance is large compared to the wave amplitude, the order of magnitude analysis is then performed to derive the boundary layer equations in the orthogonal wave-following curvilinear coordinates. The boundary layer equations are also decomposed into the mean, wave-induced and turbulent components. Implementing this triple decomposition leads to the appearance of double velocity correlations of wave and turbulent velocities that requires being treated with the assumption that wave-induced and turbulent velocity components are all of the same order.