Stokes drift through corals

We investigate the all-penetrating drift velocities, due to surface wave motion in an effectively inviscid fluid that overlies a saturated porous bed of finite depth. Previous work in this area either neglects the large-scale flow between layers [Phillips (1991)] or only considers the drift above the porous layer [(Monismith (2007)]. We propose a model where flow is described by a velocity potential above the porous layer, and by Darcy's law in the porous bed, with derived matching conditions at the interface between the two layers. The damping effect of the porous bed requires a complex wavenumber k and both a vertical and horizontal Stokes drift of the fluid, unlike the solely horizontal drift first derived by Stokes Stokes (1847) in a pure fluid layer. Our work provides a physical model for coral reefs in shallow seas, where fluid drift both above and within the reef is vitally important for maintaining a healthy reef ecosystem [Koehl et al. (1997), Monismith (2007)]. We compare our model with measurements by Koehl \&Hadfield (2004) and also explain the vertical drift effects described in Koehl et al. (2007), who measured the exchange between a coral reef layer and the (relatively shallow) sea above.


INTRODUCTION
The importance of wave-driven flows and their role in mass transport within the world's oceans has long been recognised on a macroscopic scale, carrying waste and driftwood [7]. One of the key processes which is understood to facilitate this large-scale transport ('drifting') is Stokes drift, a net motion in the direction of wave propagation, first introduced by Stokes [3], arising from the difference between the Lagrangian and Eulerian velocities of fluid particles undergoing oscillatory motion.
The case of Stokes drift arising from surface gravity waves is well-studied, and a detailed exposition of its derivation can be found, for example, in Phillips [8]. For a fluid particle with initial position x 0 and a Lagrangian velocity field u L (x 0 , t), the particle's position at time t is and the Eulerian velocity u (x, t) is given, upon expansion of a Taylor series, by u (x, t) = u (x 0 , t) + t 0 u L (x 0 , s) ds · ∇u + . . . .
(2) Following the approach of Longuet-Higgins in [9], this first-order difference term can be time-averaged over a period of oscillation to give the Stokes drift velocity which we will henceforth denote u S . Letting . . . denote an average over an oscillation, and dropping the subscript from the Lagrangian velocity, we obtain In the case of linear water waves on an ocean surface with mean position z = 0 and total depth D [described by z = A exp {i (kx − ωt)}, where the x coordinate is horizontal], [8] shows that the Stokes drift effect is entirely horizontal, with magnitude which decreases exponentially with depth. However, in many oceanographic contexts, it is not realistic to treat the water as a layer of fixed depth D with an impenetrable boundary at z = −D -perhaps the simplest example of this is the case of shorelines where the sea overlies a saturated bed of sand, and is known to induce some flow within the sand (as discussed by Phillips [1]). One other such potential extension is to a coral reef underlying the ocean. Unlike the dense sand beds discussed by Phillips, flow between the two layers is of great importance as the reef layer is much more permeable on average. A typical coral reef also has a much deeper vertical extent than the porous sand layer typically modelled, before one reaches an effectively solid rock layer below.
The importance of mass transport into and throughout coral reefs is indisputable. Monismith [2] states the importance of flow throughout the reef in trapping nutrients and plankton, as well as mass transfer's role in preventing coral bleaching events [10]. Further studies have investigated the effects of such flows on larval accumulation [11] and underlined the importance of vertical flows as well as the purely-horizontal effects seen in the absence of a porous layer [6].
Fluid-mechanically, however, existing studies of the hydrodynamics of coral reefs tend to consider the porous layer as a boundary condition, exerting drag on the flow above (see, for example, Rosman & Hench [12]). Monismith [2] uses the model of Longuet-Higgins & Stewart [13] to treat the effect of propagating waves as a body force on the mean flow, and also cites a number of studies on wave-breaking on the faces of reefs. However, on smaller scales, there has not, to the present authors' knowledge, been any studies on Stokes drift within the coral reef itself, and how the damping effect of this porous layer affects the drift velocities both above and within, aside from the model discussed here and first mentioned in [14].
In this article, we will describe a two-layer model, coupling potential inviscid theory (as first treated in [3]) above a viscously-dominated porous layer, where the flow is governed by Darcy's law [1], to derive not only expressions for how incident waves are damped by the presence of a porous layer, but also analytical expressions for the Stokes drift velocities that result. We find that the difference between Eulerian and Lagrangian velocities gives rise not just to a horizontal drift effect, but also to the vertical drifts observed by Koehl et al. [6], and offer quantitative explanations for this behaviour in terms of the damping of the waves. Finally, we discuss the applicability of the model to complicated real-world coral reefs, and compare the predictions of our model with measurements made in the field.
We consider a layer of water of total depth D bounded below by an impenetrable floor at z = −D. This lower boundary directly underlies a saturated porous medium, of permeability Π (z), which occupies the space between z = −D and z = −d. Such a configuration, with surface waves z = η (x, t), is shown in figure 1. In the case d → D or κ → ∞, we recover the classical result derived by Stokes, an example of which is shown in equation (4), which will later serve to provide a check on results.

Flow above the porous layer
Above the porous layer, we follow [3] in letting the fluid velocity be described by a velocity potential φ, such that u = ∇φ. Incompressibility then imposes that ∇ 2 φ = 0, which is to be solved in z > −d subject to boundary conditions which we will determine.
We start by making the usual assumption that where the wavenumber k may be complex, to allow for damping effects, which we will later see are key to describing the novel internal vertical drifts driven by horizontal wave motion at the surface. At the free surface z = η (x, t), the dynamic boundary condition indicates which, making the assumptions of linear theory, can be linearised to give ∂φ/∂z = ∂η/∂t at z = 0. This condition implies that φ must have the same form of x-and t-dependence as η (x, t), such that where the real part is implicit (henceforth we will always make this assumption). Solving Laplace's equation for φ with this assumed form allows us to determineφ up to constants to be determined. Thence, and the linearised form of equation (6) shows that Furthermore, we can impose a kinematic boundary condition arising from the Bernoulli equation for unsteady potential flow (see, for example, Batchelor [15]), such that at z = η (x, t), for f an arbitrary function of time. Without loss of generality, this f (t) can be set to equal zero, absorbing constants, and so the linearised boundary condition becomes providing a second relation between the constants C 1 and C 2 , namely that Flow within the porous layer The modelling of flow within porous layers is based on the use of Darcy's law (Phillips [1]), relating the volumetric flux to pressure gradients within the fluid. For a porous medium with isotropic permeability Π (z), this flux q is given by where p is the pressure field within the porous medium.
In order to facilitate matching with the oscillatory pressure field above the porous medium, we postulate that this pressure field must take the form for constant C 3 and a functionp (z) to be determined. Because we consider incompressible flow, and the matrix of the porous layer is fixed, ∇ · q = 0, so ∇ 2 p = 0. Further imposing the condition of zero vertical flux across the impermeable boundary at z = −D, the pressure field must take the form (16b)

Matching the two layers
Having defined our model, we are left with four constants C 1 , C 2 , C 3 and C 4 to be determined. We combine equations (9) and (12) to find that which, as expected, gives the well-known dispersion relation ω 2 = gk tanh kD ( [16]) in the case of no porous layer (fixing C 2 = kD to ensure zero vertical velocity at the lower boundary).
Matching of the inviscid flow above the porous layer and the viscously-dominated Darcy flow within is a nontrivial task. The velocities above the porous layer have a different intrinsic meaning to the spatially-averaged volume fluxes within the porous layer, and it is not immediately clear whether considering the classical matching conditions of velocities and pressures is a valid approach. It is, however, clear that mass conservation requires vertical velocities to match at the interface between the two flow regimes, allowing us to derive the matching condition It is clear that we should not match tangential velocities at the interface, however. Beavers & Joseph [17] remark that there is a 'slip' discontinuity between layers, and the forces exerted by the reef on the flow could complicate matters, for example by imposing shear stresses (as is the case in many existing models such as Rosman & Hench [12]). The model in [17] also requires continuity of pressure at the interface, and the linearisied Navier-Stokes equations suggest, in z ≥ −d, Then, if we take the (constant) pressure above the surface waves to be zero, (20) Hence, matching this with the expression in (15), it is seen that C 3 = ρgd and Deriving the dispersion relation This condition now means that the solution can be fullydetermined. Starting from (17), which determines C 2 , we can also determine C 1 from (12), namely that (22) It is then straightforward to determine the value of C 4 using either (18) or (21). This gives (23) Perhaps of more interest, however, (18) and (21) can be divided to give where Π I = Π (−d) is the interfacial permeability. This is a dispersion relation linking the frequency ω of surface waves to their complex wavenumber k. It is important to note from the outset that only the value of permeability at the porous layer interface, Π I , is of any importance in this relation; this is an artefact of the fact that there is no stress matching condition and instead all matching conditions depend only on values at the boundary.

Special cases of the dispersion relation
Define the dimensionless constant J = ρωΠ (−d) /µ such that the dispersion relation of equation (24) can be written (25) Owing to the nature of this dispersion relation, it will usually need to be solved numerically for k given ω and the parameters of the model. To investigate the properties of the relation, define such that, for fixed ω, zeros of f correspond to suitable values of the wavenumber k. Making the assumption that waves propagate in the positive x-direction and are therefore damped in this direction (i.e. their amplitude decreases as x increases), we will only seek (the physically-relevant) solutions with Re (k) > 0 and Im (k) > 0.
An initial investigation shows that it is possible to find multiple solutions for k that satisfy this assumption, even with all of the other parameters fixed. Figure 2 shows that there exists one propagating solution with a small imaginary part and larger real part, and one almost evanescent solution with a very small real part and larger imaginary part. For our purposes here, we will consider the propagating solutions because we would expect the almost evanescent solutions to decay quickly as they are damped over the reef.

Limit of no porous layer
As mentioned above, one would expect that the classical dispersion relation ω 2 = gk tanh kd [16] would be recovered in the limit of no porous layer. Indeed, taking either d → D or Π (−d) = 0 (i.e. setting the porous layer interface to have zero permeability) sets the lefthand side of equation (25) to zero, such that tanh arctanh ω 2 gk − kd ⇒ ω 2 = gk tanh kd. (27)

Limit of small frequencies
In the case of low-frequency waves, we would expect k to be small, giving longer wavelengths. Using the addition formula for tanh, and tanh x ≈ x for small x, so  simplifying to (30)

Limit of high frequencies
Conversely, we would expect k → ∞ as ω → ∞ and therefore, again starting from equation (28), but because J is real, ω 2 = gk, or, alternatively, k = ω 2 /g. Note here that k is real in this limit, so the asymptotic expression for Im (k) is simply zero to leading order. Considering the next order, as detailed in the appendix, gives the asymptotic expansion where E 1 = e −2ω 2 d/g and E 2 = e −2ω 2 (D−d)/g .

STOKES DRIFT VELOCITIES
Following the approach detailed in [8], given the expressions for velocities both above and within the porous layer, it is possible to calculate the Stokes drift velocities in these layers. Provided that particles do not cross between these two layers, a complication which will be later discussed in more detail, this gives an indication of the wave-driven fluid velocities.

Above the porous layer
Given the expression of (8) for the velocity potential, it can be seen that taking real parts as discussed above. Also, again taking real parts, Using the result that Re ae iωt Re be iωt + c = 1 2 Re ab , where a, b and c are time-independent complex numbers, and equation (3), we can find the Stokes drift velocity above the porous layer where the arguments of cosh and sinh, both kz + C 2 , are omitted for brevity. Then, we can use the identities to see that Within the porous layer Analogously within the porous medium, recall that is an analogue for velocity u, so Then, Thus, we can write the two components u S = (u S , v S ) in the form and It is also of interest to consider the case where d → D and there is no porous layer. As was seen in the case where there is no porous layer, the dispersion relation becomes ω 2 = gk tanh kd and thus C 2 = kd. This means that Then, above the porous layer, equation (37) becomes which simplifies precisely to the form of (4) in this limit, as would be expected.

STOKES DRIFT IN CORAL REEFS
It becomes useful at this stage to introduce representative dimensional paramters -dependence on properties of the porous layer such as the permeability and its geometry cannot be encapsulated in a simple dimensionless parameter like J above. As discussed in the introduction, an understanding of the hydrodynamics of coral reefs is of great importance, and there is a wealth of data available on their structure and composition. Though this can vary widely, we choose the well-documented reefs in Kaneohe Bay, Hawai'i as a basis for this modelling, which are primarily comprised of Porites compressa coral [5].
On the assumption that the water's kinematic viscosity is ν ≈ 10 −6 m 2 s −1 and the flow within the reef has characteristic velocity scale U = 10 −2 ms −1 [5], the Reynolds number of the flow in the porous layer is where l is a characteristic length scale for the problem -the grain diameter of the reef. Provided this diameter is of the order 10 −3 m (i.e. millimetre-scale), we can expect Darcy's law to be a valid description of the flow therein (Bear [18] states that Darcy's law is a valid model for Re 10). The reefs on Kaneohe Bay are especially simple to model, given that Koehl & Hadfield [5] remark that waves break on the outside of the reef and then propagate steadily, and uniformly, over the water surface, satisfying our modelling assumptions. For the sake of simplicity, we start by modelling reefs where the permeability is a constant Π (z) = Π 0 . The parameter values which are used in this model are summarised in table I.

Parameter
Value Density of water, ρ   (2019). Although it is difficult to determine the value of D, measurements suggest that the flow of water is negligible more than 2 metres below the ocean surface.

Flow velocities
In this particular scenario, the numerically-calculated wavenumber is k = 0.56 + 0.28i. As an initial check of the validity of our approach, we can compare the velocity magnitudes both above and within the reef as calculated by our model [using equations (8) and (16)] to the peak values presented in table 1 of [5]. This comparison is shown in table II.
Depth Measured peak (ms −1 ) Predicted peak (ms −1 ) z = −0.48m (above) 0.12 0.11 z = −0.91m (within) 0.03 0.08 TABLE II: A comparison of the peak horizontal velocities predicted by our modelling and reported in [5], showing approximate agreement both above and within the reef layer.

Stokes drift velocities
Plotting the paths of fluid particles undergoing oscillatory motion in this model, as shown in figure 5, we can see not only a drift effect in the positive x-direction (i.e. in the direction of wave propagation), as would be expected from Stokes' theory, but also a new vertical drift effect. This drift can be seen to arise from the damping of the waves, as the amplitude of vertical motion decreases with horizontal distance.
As an individual fluid 'parcel' moves forwards, it also moves downwards, before the direction of its horizontal velocity changes and it moves backwards, behind its original position. At this point, the parcel is moving upwards, but the magnitude of this upwards velocity is greater than that of the downward section of motion, owing to the damped amplitude. Therefore, over an entire period of motion, the fluid parcel not only experiences a net drift in the horizontal direction, but also a net upwards vertical drift. Such a simplified model, however, breaks down close to the boundary between the reef and the fluid overlying it. An example of this behaviour is shown in the second plot of figure 5, indicating that the expressions of equations (37) and (41) are no longer valid, and a direct numerical approach is required. In the case where these equations are valid, however, figure 6 shows representative horizontal and vertical drift velocities, of the order of a few millimetres per second, largely in agreement with the mean (drifting) velocities measured by Koehl & Hadfield [5].

Varying reef depths
It is apparent that changing the depth of the porous layer underlying the water will affect the amount by which waves are damped.

Role of spatial variation in permeability
Webber & Huppert [14] discuss in more detail the effect of taking the permeability to be a function of depth Π (z) as opposed to a constant Π 0 . Such a consideration is important physically, with marine reefs decreasing in permeability with depth. The dispersion relation is local, however, only dependent on the value of the permeability at the interface, and therefore any non-uniform permeability structure only manifests itself in the flows within the reef.
The expressions for the Stokes drift velocity within the porous layer, equations (41), show that permeabil- and so we can only reasonably neglect this term in cases where the permeability changes are slight or the waves are propagating at especially large or small wavenumbers.
A second application of such spatial variation is seen when an algal turf layer overlies a coral reef -this can have dramatic effects on flow within the reef itself and the habitat therein (Koehl et al. [4]). In this case, a thin layer of algal material, typically of lower permeability, overlies a reef which we assume to have uniform permeability. Restricting our attention to coralline algae, where issues of compression of the layer can be neglected [19], we see that varying the permeability of this layer can have a dramatic effect on the transport, both horizontally and vertically, within the reef, as discussed in [14].

CONCLUSION
It has been shown that an analysis of surface gravity waves can be extended to a system where the fluid sits atop a saturated porous bed, incorporating the damping effect of this bed on the fluid motion. By matching a viscously-dominated Darcy flow law in this lower porous layer with inviscid potential flow above, a complete description of flow both above and within the bed can be derived, from which one can derive Stokes drift velocities in both layers.
In addition to the well-known horizontal drift in the direction of wave propagation first derived in [3], it is seen that the damping of the waves on the surface leads to a vertical drift in both layers. Though this model does not apply for fluid that crosses the boundary between the porous medium and the overlying water, in both regions the calculated drifts match measurements made by Koehl & Hadfield [5] to a good degree of accuracy. Further field measurements in different reef locations around the world are needed to both confirm the generality of this result and provide more accurate parameter values for modelling.
This different, and possibly simplified, model can be extended to consider permeable layers where the permeability varies with depth, as explored in more detail in Webber & Huppert [14], and could then be further extended to consider other features of 'real-world' coral reefs, including variable topographies and sloping sea floors, as well as incorporating existing work on the breaking of waves on reef lagoons and its effects on mass transport [2]. However, we would expect to still see this novel vertical drift effect as a major contributor to mass transport in reefs.
Joseph Webber is thankful to the Heilbronn Fund at Trinity College, Cambridge, for funding much of his research into this topic. Both authors also thank Professor Mimi Koehl for a stimulating seminar in Cambridge that raised the problem considered in this paper, her continued encouragement and her provision of field measurements.
Appendix: Large-ω limit of Im (k) It is seen from (31) that, to leading order, k is real with value ω 2 /g as ω → ∞, but it is also desirable to have a leading order approximation to Im (k) in this limit to understand the damping of the waves. Using the fact that tanh x ≈ 1 − 2e −2x as x → ∞, (28) becomes where E 1 = e −2kd ≈ e −2ω 2 d/g and E 2 = e −2k(D−d) ≈ e −2ω 2 (D−d)/g to leading order. We then postulate that k = ω 2 g (1 + K) with K 1, (A.47) and work to first order in K. This results in Solving for Im (K) gives (A.49) Finally, using the fact that E 1 E 2 E 1 , E 2 and E 2 1, it is found that as ω → ∞, shown to be a good fit for even relatively small ω in figure 3.