Gerstner waves in the presence of mean currents and rotation

We present a Lagrangian analysis of nonlinear surface waves propagating zonally on a zonal current in the presence of the Earth’s rotation that shows the existence of two modes of wave motion. The first, ‘fast’ mode, one with wavelengths commonly found for wind waves and swell in the ocean, represents the wave–current interaction counterpart of the rotationally modified Gerstner waves found first by Pollard (J. Geophys. Res., vol. 75, 1970, pp. 5895–5898) that quite closely resemble Stokes waves. The second, slower, mode has a period nearly equal to the inertial period and has a small vertical scale such that very long, e.g. $O(10^{4}~\text{km})$ wavelength, waves have velocities etc. that decay exponentially from the free surface over a scale of $O(10~\text{m})$ that is proportional to the strength of the mean current. In both cases, the particle trajectories are closed in a frame of reference moving with the mean current, with particle motions in the second mode describing inertial circles. Given that the linear analysis of the governing Eulerian equations only captures the fast mode, the slow mode is a fundamentally nonlinear phenomenon in which very small free surface deflections are manifestations of an energetic current.

In 1802 the Austrian physicist and engineer Franz Josef von Gerstner published a remarkable solution to the water-wave problem, remarkable because his solution, obtained through consideration of the Lagrangian equations of motion, is an exact solution to the Euler equation with a free boundary -see Gerstner (1802); the solution was rediscovered by Froude (1862) and Rankine (1863). In contrast, Stokes' solution published 45 years later only approximately satisfies the free surface condition; see Stokes (1847). As argued by Lamb (1932), because the Gerstner flow is rotational, it cannot be generated by conservative forces and so has largely been neglected in the literature. However, there exist observations that do suggest the possibility of Gerstner waves: Monismith et al. (2007) report laboratory experiments for which † Email address for correspondence: adrian.constantin@univie.ac.at 512 A. Constantin and S. G. Monismith mean (wave-averaged) Eulerian flows exactly cancelled the Stokes drift, behaviour that is the lowest-order observable difference between Gerstner and higher-order (nonlinear) Stokes waves. Besides their own experiments, Monismith et al. (2007) showed that other experiments, notably those of Swan (1990), also displayed this behaviour, although other work in the same flume reported in Swan & Sleath (1990) showed mean velocity profiles that agreed better with higher-order Stokes wave theory. Weber (2011) developed a surface boundary layer theory including surface tension for both Stokes and Gerstner waves, and compared the results to laboratory experiments of Law (1999). Objective examination of his results suggests that both theories are of nearly equal accuracy.
The necessity of mean Eulerian flows that cancel the Stokes drift has also been suggested by Ursell (1950) and Hasselmann (1970), both of whom examined wave propagation on a rotating Earth. In particular, Hasselmann showed that at steady state, there would be a wave-averaged force directed along wave crests (that is, perpendicular to the wave propagation direction) due to distortion of the planetary vorticity by the waves. The so-called Hasselmann force (see e.g. Lentz et al. 2008) represents the strong rotation limit of the vortex force identified by Craik & Leibovich (1976) that arises from interaction of the wave velocity field with the mean vorticity field, in this case associated with the vertical vorticity (which is nearly equal to the planetary vorticity, f) and the horizontal, depth-variable Stokes drift. Pollard (1970) derived an extension of Gerstner's theory applicable to rotating flows, showing that for gravity waves, the effect was a very slight (1 part in 10 4 ) correction in the dispersion relation, and a very slight cross-wave tilt to wave orbital motions. Pollard noted however the important aspect of this solution: it satisfies Hasselmann's mean flow constraint; i.e. as with other Gerstner-type wave motions (e.g. Constantin 2001Constantin , 2012Hsu 2015; Henry 2016a), particle trajectories are closed. Moreover, like non-rotating Gerstner waves, the difference between Stokes waves and Pollard's wave is only apparent in a change in the mean flow at second order in wave slope, a difference that is difficult to discern observationally although there are suggestions for such a difference in the literature, notably simultaneous ocean observations of Eulerian and Lagrangian motions reported by Smith (2006). Xu & Bowen (1994) used an approximate, linear Eulerian analysis to show that because motions parallel to wave crests are in phase with vertical velocities, there is a weak wave stress. In the absence of viscosity, this wave stress produces the same steady Eulerian mean flow that Pollard found for Coriolis-affected Gerstner waves. However, the observations of Smith (2006) point to a second weakness of Gerstner wave theory: no one has developed a theory for Gerstner wave groups or has proved that such groups can exist. Indeed, the generation of the Stokes drift cancellation via the Hasselmann force presumably requires a time that is comparable to the inertial period (cf. the discussion in Pollard (1970) or Xu & Bowen (1994)), and thus is not likely to be applicable to wave group dynamics. Nonetheless, while there are no known solutions to initial value problems that produce Gerstner waves, neither has it been proven that such solutions do not exist. At the same time, we are unaware of any explanation based on Stokes wave theory that can describe Smith's observations that mean Eulerian flows developed that cancelled the Stokes drift associated with passing wave groups.
An important aspect of Pollard's solution is that it does not include the effects of mean currents. In the present paper we extend Pollard's solution to include a depthinvariant mean current. The dispersion relation for the wave motions in this case is defined by a sixth-order polynomial that we find has 4 real roots. The two larger of these roots describe the waves that Pollard found, i.e. Gerstner waves as slightly altered by rotation, whereas the two small roots define a new type of wave that decays strongly with depth from the free surface such that very long wavelength motions still satisfy the requirements of the Lagrangian deep water-wave analysis, i.e. that the motions are not influenced by the presence of the bottom of the ocean.

The governing equations
We investigate surface water waves propagating zonally in a relatively narrow ocean strip less than a few degrees of latitude wide and so we regard the Coriolis parameters as constant. Here φ denotes the latitude and Ω = 7.29 × 10 −5 rad s −1 is the rate of rotation of the Earth. For example, the values f =f = 10 −4 s −1 are appropriate to 45 • latitude in the Northern hemisphere; see Gill (1982). In terms of the Cartesian coordinate system with the zonal coordinate x pointing east, the meridional coordinate y pointing north and the vertical coordinate z pointing up, the governing equations in the f -plane approximation we solve are the Euler equations coupled with the equation of mass conservation and with the condition of incompressibility Here t is time, (u, v, w) is the fluid velocity, P is the pressure, g = 9.8 ms −2 is the (constant) gravitational acceleration at the Earth's surface and ρ is the water's density. The appropriate boundary conditions for deep water waves are the dynamic boundary condition P = P atm at the free surface z = η(x, y, t), where P atm is the constant pressure of the atmosphere at the surface of the ocean, and w = η t + uη x + vη y on z = η(x, y, t), (2.6) together with the requirement that the wave motion is insignificant at great depths. For our purposes it is convenient to use the Lagrangian framework. We will show that an explicit solution to the governing equations (2.2)-(2.6), with a depth-invariant mean current, is provided by specifying, at time t, the positions of the fluid particles in terms of the labelling variables (q, r, s) and the real parameters a, b, c, d, k, m. Here c 0 is the current velocity, k = 2π/L is the wavenumber corresponding to the wavelength L, c is the wave speed, and b, d and c are suitably chosen in terms of k > 0, m > 0 and a > 0; we shall see, cf. (2.30) below, that b > 0 while the sign of d varies (d > 0 in the Southern hemisphere, d < 0 in the Northern hemisphere and d = 0 on the Equator). Anticipating the relation b 2 = a 2 + d 2 , cf. (2.28) below, we note that (2.7) represents a wave with crests parallel to the y-axis, propagating in a West-East direction (see figure 1), with the particles moving in trochoidal orbits -curves described by a fixed point on a circle of radius b e ms , with a centre at (q − c 0 t, r, s), rolling at constant speed c 0 eastward along the x-axis, in a plane that is at an angle arctan(d/a) to the vertical (see figure 2). Note that a circle in three dimensions is parametrized by six numbers (two for the orientation of the unit vector ν which is normal to the plane of the circle, one for the radius δ of the circle and three for the circle's centre C): where µ is a unit vector orthogonal to ν and θ is the angular position on the circle; we recover (2.7) with (2.9) Since a > 0, b > 0 and the sign of d changes across the Equator, an examination of the components of the vector µ shows that only the equatorial waves (for which d = 0) are symmetric with respect to the local vertical in the meridional direction, all other waves being tilted towards the pole in their hemisphere with an angle of inclination from the vertical that increases with latitude. Let us also point out that the above considerations show that the labelling variables do not represent the initial position of the particle they define. In (2.7), the labelling variable q runs over the entire line, while r covers an interval (−r 0 , r 0 ) for some r 0 , and at every fixed value of r we allow s s 0 (r) for some value s 0 (r) < 0 that corresponds to the free surface, to be determined later on: at the fixed latitude r, setting s = s 0 (r) in (2.7) provides us with a parametrization of the surface wave at time t as q varies. Beneath the surface waves, the constraint m > 0 forces the amplitude of the vertical oscillations of a particle to decay exponentially with increasing depth. With c 0 interpreted as the mean flow velocity, the wave-current interaction specified by (2.7) therefore presents a negligible wave component at large depths -in such regions the flow is merely a pure horizontal current. Assuming constant density ρ 0 throughout the fluid, let us find the relations between the various parameters in (2.7) required for it to represent an exact solution. Set The Jacobian of the map relating the particle positions to the Lagrangian labelling variables is obtained by computing the determinant D of the matrix FIGURE 1. While the equatorial waves are symmetric with respect to the local vertical, all other waves are slightly tilted towards the pole in their hemisphere. There is no change in the y-direction but in some regimes the wave has pronounced flat troughs and sharp crests: for φ = 60 • S, note the trochoidal shape of the wave with c 0 = 0.3 m s −1 , k = 6.28 × 10 −2 m −1 , a = 10 m, seen in (a), versus the sinusoidal shape for the smaller amplitude (and much longer) wave with c 0 = 0.3 m s −1 , k = 6.28 × 10 −2 m −1 , a = 2 cm, seen in (b). FIGURE 2. Wave-induced particle motions on the water surface for a = 2 cm, c 0 = 0.3 m s −1 , k = 6.28 × 10 −7 m −1 and φ = 60 • S: (a) trajectory as seen in a fixed coordinate system; and (b) trajectory as seen in a coordinate system moving with the mean flow, where the closed orbit that typifies Gerstner waves is readily apparent.

A. Constantin and S. G. Monismith
The flow is volume-preserving if and only if this determinant is time independent. Since D = (am − bk) e ξ cos θ + 1 − abmk e 2ξ , this yields the relation am − bk = 0. (2.12) Furthermore, in addition to (2.12), we also have to impose the constraint 1 − abmk e 2ξ = 0 throughout the flow, (2.13) in order to ensure, by means of the inverse function theorem, that the labelling (2.7) represents a valid local diffeomorphic change of coordinates. Due to (2.12) and the fact that we allow all values ξ s 0 (r) < 0, (2.13) is equivalent to a 2 m 2 e 2ms 0 (r) < 1. (2.14) A glance at (2.7) confirms that ae ms 0 (r) is the surface wave amplitude at latitude r, so that (2.14) expresses the fact that 1/m is an upper bound for the wave amplitude.
Let us now write the Euler equation (2.2) in the form where D/Dt stands for the material derivative. From (2.7) we can compute the velocity and acceleration of a particle as respectively. We can therefore write (2.15) as Gerstner waves in the presence of mean currents and rotation

517
The change of variables The equality of the mixed partial derivatives of P with respect to the labelling variables (q, r, s) leads to relations fb + dk c = 0 (2.25) and mbkc 2 + mf dc = ak 2 c 2 , (2.26) if we take (2.12) into account. With the constraints (2.25)-(2.26) in place, note that, for every constant P 0 , the gradient of the expression with respect to the labelling variables is given by the right-hand side of (2.22)-(2.24). Therefore the free surface is identified at the meridional distance r from the reference latitude by specifying a value s 0 (r) of the label s for which the expression in (2.27) equals to P atm , thus validating (2.5) and (2.6). This is possible if and only if the righthand side of the expression (2.27) is time independent, in which case the value of s 0 (r) is found by solving the implicit functional equation P = P atm for a fixed value of r. The time independence of the expression (2.27) holds if (2.31) The remaining relation between the parameters, (2.29), takes on the form (2.32) Squaring (2.32), we can take advantage of (2.31) to obtain the dispersion relation (2.33) The non-zero roots of (2.33) are theoretically admissible wave propagation speeds if c 0 = 0. However, if c 0 = 0, then the physically relevant roots of (2.33) are those that differ from c = ±f /k; thus, the latter are extraneous roots, ruled out by (2.31).
For equatorial waves we have f = 0 andf = 2Ω, so that (2.33) reduces to In the absence of an underlying current (for c 0 = 0) we recover a result of Hsu (2015).
On the other hand, at mid-latitudes we have f = 0 and it is convenient to write the dispersion relation (2.33) as in terms of the non-dimensional variables We can consider the extent to which surface wind waves might be affected by rotation by noting that their wavelengths typically range from 16 m to 400 m, so that √ k/g ranges from 0.2 to 0.04 s m −1 , while √ gk ranges from 2 s −1 to 0.4 s −1 . Note that X = O(1) corresponds to a wind-wave speed c in the range from 18 km h −1 to 80 km h −1 . Given that f is of the order of 10 −4 s −1 , we see that ε = O(10 −4 ), and so, as concluded by Pollard (1970), the effect of rotation is small. Typical currents do not exceed 1 m s −1 , so that X 0 = O(10 −1 ).
The non-dimensional dispersion relation at mid-latitudes, (2.36), can be written out explicitly as a polynomial equation of degree six P(X) = 0, where In linear wave propagation theory one can take advantage of the superposition principle by seeking a Fourier decomposition of an oscillatory wave motion in eigenmodes e ik(q−ct) , with the following interpretation for the disturbance from a state of equilibrium, at a fixed wavenumber k > 0:  (i) a real root c of the dispersion relation generates two oscillatory eigenmodes of fixed amplitude, cos[k(q − ct)] and sin[k(q − ct)], which differ only by a phase; (ii) a complex root c = α + iβ, with α ∈ R and β ∈ R \ {0}, is indicative of evanescent or amplifying modes (according to whether β > 0 or β < 0), e −kβt cos[k(q − αt)] and e −kβt sin[k(q − αt)], propagating at the speed α with an amplitude that decays or is amplified in time by the exponential factor k|β|.
In contrast to this, (2.33) is the dispersion relation for the nonlinear wave-current interaction specified by the flow (2.7), and imaginary complex roots c ∈ C \ R have to be discarded.

Solution of the dispersion relation and a description of the resulting wave motions
The exact roots of the polynomial P(X) that defines the dispersion relation at midlatitudes cannot be found analytically. An example plot of P with parameters chosen to match the Antarctic circumpolar current (ACC), i.e. for 65 • S f ≈ −1.3 × 10 −4 s −1 and a mean current c 0 = 0.3 m s −1 (see Constantin & Johnson (2016) for field data), is shown in figure 3 where two roots near ±1 (figure 3a) can be seen and two roots near ±ε (figure 3b) can also be seen. One can prove that the polynomial P has exactly four real roots, two positive and two negative (see the discussion in the Appendix).
Examination of figure 3 suggests that accurate analytical approximations of the roots to P(X) can be found by a perturbation expansion in ε: (3.1) 520

A. Constantin and S. G. Monismith
After substitution into P(X) and collecting different powers of ε, we find that at O(1) Note that the roots to the sixth-order polynomial that are near ±1 correspond to Gerstner waves that are very slightly modified by rotation, of the type found by Pollard (1970). Since only real roots are relevant, in the rest of what follows we focus on the O(ε) roots. Assuming x 0 = 0, then 1 − x 2 1 = 0 at O(ε), which derives from the O(ε 2 ) expansion, so that (3.4) Consideration of the O(ε 3 ) terms shows that x 2 = 0 and finally, from the O(ε 4 ) terms, we find that X 2 0 ∓ 2x 3 = 0, that is, Thus an approximate solution to (3.1) shows two possible small roots that are approximately (3.7) For example for a 500 m long wave (k ≈ 1.3 × 10 −2 m −1 ) on the ACC, c ≈ 1 cm s −1 since f ≈ −1.25 × 10 −4 s −1 . On the other hand, for a 10 4 km long wave (k ≈ 6.2 × 10 −7 m −1 ), c ≈ 200 m s −1 , comparable to typical open ocean barotropic long wave speeds. Note that the period T for these waves is quite close to the inertial period since i.e. these slow waves are inertial waves, with inertial period T i . The correction to X at O(ε 3 ) is important in light of (2.31), which can be written as This gives a wavenumber that is much larger than k since g/( fc 0 ) 1. On the other hand, this also means that waves that are quite long with respect to the depth can produce surface-trapped motions that satisfy the fundamental assumption that the wave motions are not affected by the ocean's bottom. For example, a wave with a wavelength of 10 4 km (k ≈ 6.28 × 10 −7 m −1 ) on a mean flow of 0.3 m s −1 produces a vertical decay scale of approximately 6 m. We note that using satellite altimetry, White & Peterson (1996) show evidence of waves in the ACC with wavelengths of approximately 11 × 10 3 km, and so it seems possible that waves can exist in the ACC that are sufficiently long to produce physically plausible decay scales according to our model. Moreover, consideration of the vertical scale also highlights the importance of the mean current in that if c 0 = 0, then (2.33) and (3.9) with c ≈ √ g/k yield (3.11) and so the corresponding decay scale is much smaller than in the presence of a mean current.
We now compute the amplitude d of the lateral motions for the O(ε) roots: in view of (2.30), (3.7) and (3.9). Thus, while the vertical scale of these waves is small, they should have large sideways motions. Moreover, by (2.30), (3.7) and (3.12), i.e. particles travel in large inertial circles of radius b (or d) with very slight vertical motions. For the maximum possible amplitude wave (a = 1/m), the limiting lateral excursion would be d = ± g mfc 0 ≈ 1 k , (3.14) by (2.30), (3.7), (3.10), whereas the limiting longitudinal excursion is also 1/k: due to (2.30), (3.7) and (3.14). However, given that these excursions must take place over an inertial period, implying huge velocities, the surface amplitude a must be quite small relative to the (theoretically) maximum possible amplitude. Indeed, (2.12), (2.16) and (2.28) yield with amc ≈ ag/c 0 due to (3.7) and (3.10), so that it is easily seen that the maximum particle velocity is v max = c 0 + ga c 0 . (3.17) Thus, choosing v max = 1 m s −1 and assuming that c 0 = 0.3 m s −1 , the corresponding amplitude of the free surface deflection would be a ≈ 2 cm. For this velocity, the inertial circle radius would be approximately 5 km, in view of (3.13)-(3.15). For the sake of comparison the small roots to P(X) were computed numerically using the Matlab TM function 'fzero' which uses a combination of bisection, secant and 522 A. Constantin and S. G. Monismith inverse quadratic interpolation methods; see Forsythe, Malcolm & Moler (1976). A summary comparison of the perturbation expansion results for c etc. with numerically determined values is given in figure 4. These results show that the perturbation expansion is remarkably accurate, something that is hardly surprising given that ε is generally quite small.
We conclude our discussion of the wave motion induced by (2.7) by pointing out that the flow velocity is independent of the y variable. Indeed, inverting the matrix in (2.11) yields Note that (2.19) in conjunction with (2.30) ensure P y = ρ 0 fc 0 , (3.21) so that in the absence of an underlying current the pressure is also independent of y.

Linear analysis
The linearized momentum and continuity equations for a wave-current interaction with an underlying uniform current propagating eastwards at the speed c 0 , in a homogeneous layer of water (whose density we normalize to unit value), are Assuming that the flow velocity is independent of y, the exact boundary conditions on the surface z = η(x, y, t) are  Linearization of (4.2) and transfer to z = 0 gives The linearization of (4.3) is more complicated since we do not assume that the flow is irrotational. Constancy of the pressure of the free surface implies that the gradient of the pressure along the free surface defined by the tangent vector, τ , is zero: To leading order, the dynamic boundary condition (4.3) is A. Constantin and S. G. Monismith For a flow velocity that is independent of y, the last equation in (4.1) permits us to introduce a streamfunction Ψ (x, z, t), i.e.
The xand z-momentum equations in (4.1) can be combined to produce (4.9) whereas differentiation of the x-momentum equation with respect to y yields P xy = 0, allowing the y-momentum equation to be rewritten as follows, in a way that does not (4.10) We now seek travelling-wave solutions with m > 0 and k > 0, where the surface wave is defined as (4.12) Substitution of these forms into (4.10) produces the relation where λ = c + c 0 (4.14) is the wave phase speed relative to the underlying current. Likewise, (4.9) yields Combining (4.13) and (4.15) results in (4.16) noting that m is still unknown, as is c. The linearized kinematic free surface condition (4.4) connects Ψ and η: Finally, using (4.13) and (4.17), the linearized dynamic boundary condition (4.7) becomes mk 2 λ 2 +f k 2 λ − f 2 m − gk 2 = 0, (4.18) which gives Note that choosing f =f = 0 and c 0 = 0, the classical deep water gravity-wave dispersion wave is recovered from (4.18), with c 2 = g/k, since (4.16) would give c = 0 unless m = k. On the other hand, multiplying (4.18) by m and using (4.16), we findf (4.20) which yields by (4.16), after squaring both sides, the quartic polynomial equation allowing us to solve for m > 0 in terms of k. If we make the above equation nondimensional using the equation for the dimensionless vertical decay scale is The quartic polynomial P in (4.23) has two positive roots between 1 + ε 2 /2 ± ε 3 , since Note that, by Viète's formulas, the product of the four roots of P is positive, while the sum of all possible subproducts, taken three at a time, is negative. This ensures that there are no further positive roots. Since only values ξ > 0 are physically relevant, we deduce that ξ = 1 + 1 2 ε 2 + O(ε 3 ). (4.25) Substitution of (4.25) into (4.16) shows that (4.26) Therefore the linear analysis only identifies the 'fast' wave that is a standard deep water gravity wave modified very slightly by the Earth's rotation.

Conclusions
The Lagrangian analysis we present shows that including rotation and mean currents yields Gerstner wave solutions like that of Pollard (1970), excepting that we find that there are two kinds of waves possible. The first, described by Pollard, is a slightly modified Gerstner wave, that, like other Gerstner waves, does not have net wave transport and thus satisfies the steady state condition of no net transport on a rotating Earth given by Ursell (1950) and Hasselmann (1970). The second, much slower, type of wave motion requires a mean current in order that the vertical decay scale is finite; i.e. in the absence of currents, the second mode of motion would not be possible.
The second mode of motion, which we propose here to call an inertial Gerstner wave, produces waves that have extremely small vertical scales for gravity-wave wavelengths, while for very long wavelengths (thousands of km) it is one with waves with vertical decay scales of metres that can coexist with substantial currents even for relatively small amplitudes. At this time, we are unaware of any observations that might prove its existence, although given that it is at the inertial period it would be hard to detect using satellite methods (cf. White & Peterson 1996) since these tend to average out motions at time scales shorter than a few days at best. Likewise, given the ubiquity of inertial motions in the ocean, it might also be difficult to detect in data collected using fixed moorings (e.g. Gordon 1988). Nonetheless, given that the inertial Gerstner wave is an exact solution to the Euler equations and one that satisfies the exact free surface conditions on the deformed free surface, it must be seen to represent a possible form of large-scale oceanic wave motion, and thus, ultimately may be important to the dynamics of energetic currents like the ACC.
This appendix describes an analytic approach towards the identification of the location of the roots of the polynomial equation (2.38). Note that, due to the lack of significant surface waves at latitudes within 15 • from the poles, we have |F| 2 − √ 3 > 0.25, while |F| < 2.4 outside the tropical zone (at latitudes that exceed 23 • 26 16 ), with F > 0 in the Northern hemisphere, while F < 0 in the Southern hemisphere sincef > 0 but f changes sign across the Equator.