There exist different theories representing the effects of surface gravity waves on oceanic flow fields. In the past, the author has conjectured that the vertically integrated, two-dimensional fluid equations of motion put forward by Longuet-Higgins and Stewart are correct and that theories that differ from their theory cannot be entirely correct; this paper explores these differences. Longuet-Higgins and Stewart deduced vertically integrated, two-dimensional equations featuring a wave radiation stress term in the fluid dynamic, momentum equation. More recently, the author has proposed vertically dependent, three-dimensional equations that have required correction but when vertically integrated, agreed with the earlier, two-dimensional equations. This paper derives both vertically independent and vertically dependent equations from the same base and, importantly, using the same expression for pressure in the belief that the paper will contribute to the understanding and clarification of this seemingly difficult topic in ocean dynamics. An error in the classical papers by Longuet-Higgins and Stewart has been detected. Although the final phase-averaged result was correct, the error has had consequences in the development of vertically dependent equations. The prognostic equations in this paper are for the Eulerian current plus Stokes drift; toward the end of the paper these equations are contrasted with prognostic equations for the Eulerian current alone.
On the ocean surface, the correlation of wind pressure fluctuations and wave slope forces wave generation but also is a forcing term in the phase-averaged momentum equation. Phase averages of wave motion yield residual particle displacements (Stokes drift) identical to current displacements. The difference between the instantaneous wave velocity minus its phase average leads to the so-called radiation stress, first identified by Longuet-Higgins and Stewart (1962, 1964, hereinafter L-HS), which modifies the momentum equation much like turbulence Reynolds stress but, in addition, includes a component due to pressure.
L-HS, Phillips (1977, hereinafter Phillips), and Smith (2006) derived vertically integrated, wave radiation stress terms as an addition to the fluid dynamic equations of motion rendering them suitable for two-dimensional oceanic flows with a wind-driven wave surface. These derivations start from the basic vertically integrated, fluid dynamic equations after which the total velocity is partitioned into current and wave velocities resulting in phase-averaged equations applicable to vertically independent ocean models. A characteristic of the aforementioned papers is that prognostic equations are obtained for the combined (Eulerian) current and Stokes drift.
Following L-HS, this paper presents a concise review of the vertically integrated wave circulation equations but then derives the corresponding, vertically dependent equations directly from their integrated counterparts. The purpose is to demonstrate the intimate connection between the two equations sets and, at the same time, provide a corroborative and somewhat simpler derivation of the vertically dependent wave–current interaction equations in Mellor (2003). An error in the L-HS formulation that, however, disappeared in their final equations led to a similar error in the formulation by Mellor (2011) but is corrected here and in Mellor (2013a).
The basic equations of fluid dynamics are cited in section 2, wherein velocities and pressure are partitioned into their mean and wave components. Section 3 and appendix A are devoted to an expression for pressure, an important element of the ensuing derivations, which distinguishes this paper from papers that present an alternate expression for pressure. The vertically integrated equations are cited in section 4, from which the phase-averaged integral equations of L-HS are derived in section 5. The same vertically integrated equations are the basis of the derivation of vertically dependent, three-dimensional equations derived in section 6. Some relationships between the derivations of sections 5 and 6 are discussed in section 7, wherein a subtle but important error in the L-HS derivation is revealed. For convenient reference, the complete phase-averaged, three-dimensional equations are summarized in section 8, and small and neglected errors related to surface wave slope and bottom topographical slope are discussed in section 9. In section 10 and appendix B, there is discussion of theories by McWilliams and Restrepo (1999), Newberger and Allen (2007a,b), and Ardhuin et al. (2008) that apply to vertically dependent ocean models but which diverge from the findings of the present paper and, when vertically integrated, diverge from the findings of L-HS and Phillips.
Incorrect rendering of the pressure term in Mellor (2003) was later corrected in a PDF file (available online at http://shoni2.princeton.edu/ftp/glm/Corrected2003.pdf) and in Mellor (2013a). As mentioned above, a similar error occurs in L-HS; nevertheless, their final radiation stress result is correct, as is the final result in Mellor (2003).
2. The basic equations
The basic equations of fluid motion are
where the Greek subscripts denote the horizontal coordinates; z is the vertical coordinate; t is time; the three-dimensional velocity is ; p is the kinematic pressure; and g is the gravity constant. Repeated subscripts denote summation, for example, . We have omitted the Coriolis and baroclinic terms and turbulence and pressure–slope vertical momentum transfer (Mellor 2013b) in order to simplify most of the ensuing discussion; these terms will be reinserted in section 8.
In the above, is the free-surface elevation, and is the depth.
In the following, dependent variables are partitioned into current variables (whose temporal and spatial scales are large compared to inverse wave frequency and wavenumber) and wave variables, such that
The wave variables are taken to be
where we define ; is the wavenumber vector and ; a is the amplitude; is the intrinsic frequency, is the phase speed, and ; and is based on phase-averaged velocity and is defined in Mellor (2003). The phase-averaged, free-surface elevation is , and h is the depth. Also, is the water column depth. Equations (5a) – (5d) are the well-known linear wave relations; they are valid in the full instantaneous range (Mellor 2011). The dispersion relation is and is used below where appropriate; also so that, for example, (5d) can be written . Errors associated with (5) are discussed in section 9.
We denote phase averages by an overbar; thus, . For simple linear terms, the overbar and the overcaret are identical, for example, ; the phase average of product terms are those that are determined in section 4 and beyond. (Note that , and so on.) From the phase-averaged definition and using as an example, so that current variables are further defined as those where are small as stated above.
Pressure plays an important role in the following derivations. In appendix A, it is found that
where is the atmospheric pressure acting on the free surface. After insertion into the momentum equation below, the last term in (6) disappears, but terms of order (ka)4 will be neglected and represent an error of the same order. However, the phase average of (6) is exactly , in agreement with L-HS. To lowest order in ka, the free-surface pressure . Specifically, . This result is obtained from (5c) and (5d) by expanding and to obtain and, similarly, for insertion into (6).
Equation (6) differs from that in Mellor (2003), which was corrected in a PDF file (available online at http://shoni2.princeton.edu/ftp/glm/Corrected2003.pdf) and in Mellor (2013a). It also differs, for example, from the treatment of pressure in McWilliams and Restrepo (1999) and Ardhuin et al. (2008).
4. The vertically integrated equations of motion
5. The derivations as in L-HS and Phillips
The derivation of the vertically integrated equations is reviewed here. L-HS, in their pursuit of the then novel idea of radiation stress terms, generally neglect the nonlinear advective terms as in (8). The book by Phillips is rich in its omnibus description of surface (and internal) waves. However, his derivation of phase-averaged momentum equation is spread around different pages and a derivation of an important element of the stress radiation term [the last term in (16)] is missing. This term and its vertically distributed counterpart has been a source of controversy in the literature, as discussed in sections 7 and 9.
a. The continuity equation
Now, and plus higher-order terms in ka. Therefore, the integral continuity equation may be written
is the integral Stokes drift. Wave energy is . Now, define so that
b. The momentum equation
The last two terms on the right side, using (6), can be evaluated as ; generally, so that both terms can be neglected. Furthermore, we obtain
which derives from the two portions of (6) and ; phase averages of the integrals of the other portions of (6) result in terms. LH-S obtained the same result as in (13b) but for the wrong reason, as discussed in section 7.
where is the group speed, and is a Kronecker delta function (=1 when or = 0 otherwise). Following Phillips, let , after which the advective term in (15) is plus the term = O(ka)4, which can be neglected. [Phillips initially includes the term within , but subsequently the term disappears. Note that the assumption of constancy of can be avoided as seen below by vertically integrating (30).]
6. The vertically dependent equations
The present strategy, as in Mellor (2003), is to map the wave variables from x, y, and z coordinates to x, y, and coordinates such that
where is the vertical displacement of material surfaces due to waves. Notice that the phase-averaged , and thus can be recognized as the “conventional sigma” coordinate (reserving σ for frequency). For is the material surface displacement. For . Now see that, to lowest order in ka,
plus terms higher order in ka. Similarly,
In the above, definitions are
a. The continuity equation
Thus, (7) can be written
is the vertically resolved Stokes drift. After vertical integration, one obtains (10b). [Note that the derivation of (22) is alternative to the more conventional Lagrangian process (e.g., Phillips, 31–32) but with the same result.] Therefore,
where is defined so that the left side of (24) is nil everywhere in the range and therefore (24) conforms to the conventional sigma coordinate continuity equation as used in numerical ocean models [e.g., the Princeton Ocean Model (POM) and the Regional Ocean Modeling System (ROMS)]. The variable will be recognized as the component of velocity normal to surfaces of constant sigma as viewed in a Cartesian coordinate system.2
b. The momentum equation
Taking the partial differential outside of the integral and after phase-averaging results in only one extra term, ∂/∂xα, but that term can be neglected since ≪ gD (see footnote 1). Now using (18) and (19c), we have
Note that after expansion does yield plus an additional term that, however, is of order (ka)4 and is neglected. The term under the phase-averaged bar in (25) can be evaluated as since all other terms such as are nil. Thus,
The partial derivatives have been taken outside of the integral of (26) since horizontal derivatives of the integral limits are nil.
The expression under the integral can be written , where again . The term is of order (ka)4 and can be neglected. Thus, we arrive at
where the wave radiation stress is
Now add to the left side of (28) so that the integral equation is unchanged and the integrand can be written
Furthermore, the addition of the third term on the left side is necessary3 to conform to established fluid dynamic theory.
7. Relationships between the integral and differential equations
One advantage of sigma coordinate equations is that vertical boundary conditions are “built in,” and they are easy to vertically integrate. Thus, reintegration of (24) and (30) readily yields (11) and (15).
In section 5, it was seen that some of the phase-averaged wave properties were obtained from the portion of integrals whose limits were from , whereas in section 6, wave properties were imbedded in the integrand of integrals that spanned the range or . However, the two strategies are related. Thus,
so that the determinations of the integral Stokes drift in (10) and (22) are equivalent. The mathematical equality between the surface-concentrated second term and the distributed integrand on the right side may be the basis of the occasional misconception that, physically, Stokes particle drift is concentrated at the surface.
At this point, it is appropriate to point out that the derivation of L-HS introduced a subtle error in their determination of . They heuristically assumed that, hydrostatically, so that , as in (13b). In hindsight, the hydrostatic assumption is certainly incorrect. Instead, the portion of (6), , yields a negative value of , whereas the term contributes a positive . Furthermore, the vertically distributed term in square brackets in (29)—now deemed to be correct—would not be obtained using . [Following L-HS, Mellor (2011) also used the hydrostatic pressure relation that led him to erroneously replace the square bracketed term in (29) as a delta function concentrated at the surface.]
8. Summary of the full, vertically dependent equations
To provide convenient reference, we repeat the continuity equation:
To the momentum equation, add Coriolis, baroclinic, and vertical momentum transfer terms to (30) such that
where f is the Coriolis parameter. The last term on the left side is the baroclinic term—familiar to users of sigma coordinate equations—where ; is the phase-averaged density, and is a reference density. The baroclinic term is derived from the hydrostatic relation for vertically variable density after transformation to sigma coordinates. On the right side, is the conventional turbulence momentum transfer, whereas is a pressure–slope transfer introduced by Mellor (2003); it is a term absent from other theories. All of these additional terms could have been added to (2) and tracked through sections 4 and 6; their derivations are devoid of the complexity associated with the advective and wave radiation terms.
and where convenient definitions are
Equation (33b) integrates so that . Now, if the two terms whose coefficients are are combined, the result is identical to the expression for in Mellor (2003) but from an expression different than (29).
The pressure–slope transfer term in (32), , arises from consideration of a small wind pressure component in addition to (6) and proportional to ; it therefore correlates with . Thus, is recognized as surface form drag (e.g., Donelan 1999) that is projected into the water column as . For a full description of pressure–slope transfer and its interaction with turbulence transfer in the water column, refer to Mellor (2013b).
After wave parameters have been calculated, Stokes drift can be obtained from (22) and then the current may also be obtained.
Equations for scalar quantities may be written as
where T is any scalar variable, and q is the vertical flux of that variable. If T is taken to be temperature, then q is vertical heat flux including penetrative solar radiation. Equation (35) is derived in a similar way as the continuity equation in section 6.
In the preceding developments, sigma coordinates have been a convenience. However, the end results (31), (32), and (35) can be transformed to Cartesian coordinates for which reference is made to appendix A of Mellor (2005; wherein the term should everywhere be expunged).
For completeness, the wave energy equation is added; it is, ab initio, a vertically integrated equation and is
where and are source and sink terms. It has been determined that (Mellor 2003), where , and r is a weighting function so that is appropriately weighted near the surface. Equation (36) may be complemented by an equation for (e.g., Komen et al. 1994; Mellor 2003). Both equations apply to a monochromatic wave train, and both equations are found in Phillips except that his equations are restricted to a vertically constant . Whereas can include wave breaking, the transfer of wave momentum to current momentum is automatic in (32) and requires no added empiricism.
When dealing with wave spectra, , as described in Komen et al. (1994) and in many references cited therein, terms allowing for changes in frequency and wave direction should be added to (36). Spectral averages of the F terms in (34) are needed to feed into (32) and (33). This does create a burden on computational resources in addition to the need to deal with five independent variables and the need to deal with frequency shifts due to wave–wave interaction. An approximation whereby the spectral shape is parameterized has been proposed in Mellor et al. (2008) to greatly reduce the computational burden.
9. Errors of order (ka)4
In all developments of phase-averaged equations, it is assumed that temporal and spatial scales of amplitudes are small, that is, . Similar conditions apply to k and . Coupling between waves, vertical current gradients (Kirby and Chen 1989), and density gradients are excluded by stipulating that and ≪ 1; is the Brunt–Väisälä frequency.
Throughout this paper’s development, wave slope ka has been stipulated to be small. The wave terms contained in (29) or (33) and the addition of the Stokes drift term to are terms of order (ka)2. Terms of order (ka)3 are identically zero since they are coupled with triplet products of sines or cosines whose phase averages are nil. Terms of order (ka)4 have been neglected throughout the paper. Relative to the retained terms, the neglected terms represent an error of order (ka)2 = O(10−2). Note that, ab initio, waves are represented by only the first term in the Stokes series, as in (5); including the second-order term would again introduce terms of neglected order.
An error is embodied in the linear wave solutions [(5a)–(5d)], which are based on the boundary condition . Bennis et al. (2011) apparently state that the linear wave solutions cannot be invoked as done in this paper no matter how small the slope may be. This, of course, is not correct. Rather one should pose the question: what is the error associated with a finite slope? The error has been investigated in Mellor (2013a) and found to be proportional to . Therefore, consistent with the slope error of order (ka)4, one should require that . For most oceanographic applications, this requirement should not be overly restrictive, even less so if kD > 3.
10. Prognostic equations for the Eulerian velocity
The developments in this paper and those of L-HS and Phillips predominantly deal with prognostic equations for or , whereas the papers by McWilliams and Restrepo (1999) [see also the more complicated paper by McWilliams et al. (2004), which is extended to include phase-averaged long-wave equations] and Ardhuin et al. (2008) obtain prognostic equations for or that introduce a “vortex force” term. Restricting to integral equations, Smith (2006) does derive an equation for ; he first obtains an equation for from (36) and then subtracts it from to obtain . The latter does contain a vortex force term plus other terms. The derivation of his integral equations are revisited in appendix B and extended to obtain a vertically dependent, prognostic equation for . There is further discussion of similarities and differences between the equations of McWilliams and Restrepo (1999), Newberger and Allen (2007b), and Ardhuin et al. (2008) and those in this paper and the vertically integrated equations of L-HS, Phillips, and Smith (2006).
A difference is the treatment of pressure. The McWilliams and Restrepo (1999) derivation begins with the curl of the momentum equation wherein the pressure gradient term disappears; then, after considerable manipulation, the resulting phase-averaged vorticity equation is “uncurled” and a gradient pressure term can be reintroduced, but their choice of a new, phase-averaged pressure differs from that derived in appendix A and cited in (6). In McWilliams and Restrepo (1999) and Ardhuin et al. (2008), the Stokes vertical component is nonzero so that, instead of (31), the divergences of and are separately nil requiring the existence of a nonzero vertical Stokes drift. As derived in L-HS and Phillips and in this paper, Stokes drift has horizontal components [(22)], but the derived vertical component is nil; thus, currents and Stokes drift are coupled in the continuity equation [(31)].
Ardhuin et al. (2008) include terms proportional to in their formulation. This could have been incorporated in the equations of this paper. For example, add the term to the terms in the parentheses of (21) and complete the ensuing analysis. It is found that the new term relative to is of order and, to be consistent, must happily be neglected. The same is true if were taken into account in (35).
The theoretical development in Newberger and Allen (2007a) is quite complicated, but the final equation set in appendix B of Newberger and Allen (2007b) is not. These equations are meant to apply to shallow water (kD ≪ 1). The momentum equation for Eulerian current is vertically dependent but is forced by vertically integrated terms like (16) but where a vortex force term replaces the first term on the right side of (16). Also, gradients of component Stokes drift are concentrated at the surface and are a surface boundary condition for the underlying continuity equation. Additional modeling for breaking wave rollers is added to the wave energy equation.
A relatively recent paper by Aiki and Greatbatch (2013) seems to have bearing on the approach of this paper and the 2003 paper, but the paper is quite complicated and I could not extract equations such as (31) and (32) for comparison. However, they correctly complained about my earlier characterization (born out of the now discovered error in L-HS) of a portion of [(33b)] as a surface-concentrated Delta function.
The pressure–slope term in (32) is not found in papers cited in this section.
Bennis et al. (2011) have created an inviscid barotropic, shallow-water test case whereby unidirectional waves progress into an open channel. The channel bottom topography is first flat, then increases by 33% of the total depth, is flat again, and then decreases to the inflow value. The vortex force in this case is nil, and so they obtain zero currents. Stokes drift is given by (22). In contrast, a similar calculation by Mellor (2011, appendix 4) produced the same Stokes drift plus a vertical current profile stronger than the Stokes drift. This was because of the vertical variation of the second term of the radiation stress in (33)—the first term being vertically constant in this flow case. For, say the bottom radiation stress is less than 0.2% of the surface radiation stress, and therefore the variations in bottom topography have negligible effect.
The paper was significantly improved because of reviewer comments.
Derivation of (6)
where . The two differentiations have been taken outside of the integrals by using (3a). The boundary layer or hydrostatic approximation is invoked so that, in the absence of waves, ; therefore, the remaining terms in (A1) are wave terms.
The second integral in (A1) has been evaluated in detail, but by inspection
The third term on the right side equals −, cancelling in the first term. Therefore,
where is given by (5d) or (19c). In (A4), since is O(ka)2, we have retained the term here, but the phase averages of and the product , which otherwise appear in (25), are nil. Therefore, terms of order (ka)4 have been neglected. Note that whereas the phase resolved (A4) is needed, the phase average of (A4) could have been obtained directly and exactly from (A1) since the phase averages of and are nil.
Derivation of a Prognostic Equation for
This is an effort to derive a prognostic equation for . It leans heavily on the precursor derivation by Smith (2006) that became a guide through rather tortuous algebraic manipulations. The difference here is that is not assumed to be vertically constant throughout the derivation. Absent this generalization, the integral result below (B5) is identical to Smith’s result.
Excluding the Coriolis and buoyancy terms, (32) may be written as
where . Using , divide the term
When converting from E to , the term is generated resulting in the second and third term on the right side of (B3). It is has been determined that (Mellor 2003), where , and r is a weighting function so that is appropriately weighted near the surface.
Following Smith, manipulation of the second term on the left side and the first term on the right side of (B4) yields
where, again following Smith, [and ]. Notice that is vertically constant, whereas is not. However, if one assumes that is vertically constant on the right side of (B5), then ; the first term on the right cancels, and (B5) is identical to Smith’s (2.28) and (2.29) after noting that , , and . He uses subscripts i and j instead of and , the last term in (B5) is denoted by , and .
The integrand of all terms can be written
Equation (B6) is a differential equation for the current . It is a sigma coordinate equation but can be transformed to Cartesian coordinates (see appendix A of Mellor 2005). Aside from removing Coriolis and baroclinic terms—which can easily be reinstated—the term from (36) should have been included, although there is some evidence that it is small. It was not present in Smith’s derivation because he employed the wave action equation instead of (36) in which the term is subsumed but only if (Mei 1983). It is also this writer’s contention that only the pressure–slope portion of the total vertical stress divergence should balance the wave generation term, that is, . Finally, it seems appropriate that instead of , the vertically dependent form obtained from the corresponding portion of (33) should be substituted for J.
It is now recognized that (B6) coincides with the equations of Newberger and Allen (2007b) if the left side of (B6) remains vertically dependent but velocities on the right side are assumed be vertically constant. Thus, the first term and the third term (for are nil, whereas the second and third term (for shallow water, combine to yield their vertically constant
The left side of (B6) also corresponds to the equations of McWilliams and Restrepo (1999; however, advective terms are excluded), McWilliams et al. (2004), or Ardhuin et al. (2008). [In the latter case, refer to Bennis et al. (2011), which offers simplifications to the 2008 paper.] The right side agrees with all of the above authors’ equations only if one approximates or if one assumes velocities are vertically constant in which case the first term on the right side of (B6) vanishes and the second term is a vortex force term. Also, Bennis et al. (2011) have a term instead of the third term on the right side of (B6). In the present paper and in L-HS and Phillips, , as discussed in section 10.