Wave-induced Lagrangian drift in a porous seabed

The mean drift in a porous seabed caused by long surface waves in the overlying fluid is investigated theoretically. We use a Lagrangian formulation for the fluid and the porous bed. For the wave field we assume inviscid flow, and in the seabed, we apply Darcy’s law. Throughout the analysis, we assume that the long-wave approximation is valid. Since the pressure gradient is nonlinear in the Lagrangian formulation, the balance of forces in the porous bed now contains nonlinear terms that yield the mean horizontal Stokes drift. In addition, if the waves are spatially damped due to interaction with the underlying bed, there must be a nonlinear balance in the fluid layer between the mean surface gradient and the gradient of the radiation stress. This causes, through continuity of pressure, an additional force in the porous layer. The corresponding drift is larger than the Stokes drift if the depth of the porous bed is more than twice that of the fluid layer. The interaction between the fluid layer and the seabed can also cause the waves to become temporally attenuated. Again, through nonlinearity, this leads to a horizontal Stokes drift in the porous layer, but now damped in time. In the long-wave approximation only the horizontal component of the permeability in the porous medium appears, so our analysis is valid for a medium that has different permeabilities in the horizontal and vertical directions. It is suggested that the drift results may have an application to the transport of microplastics in the porous oceanic seabed. Surface waves in an upper fluid layer induce Lagrangian mean drift in a porous bottom layer. For spatially damped waves, the mean Eulerian current in the porous bed may be larger than the Stokes drift. Analysis relevant for transport of microplastics in the porous oceanic bottom layer. Surface waves in an upper fluid layer induce Lagrangian mean drift in a porous bottom layer. For spatially damped waves, the mean Eulerian current in the porous bed may be larger than the Stokes drift. Analysis relevant for transport of microplastics in the porous oceanic bottom layer.


Introduction
In a pioneering paper Reid and Kajiura [18] investigate the linear interaction between surface wave motion in a fluid layer of finite depth and the motion in an underlying saturated porous medium of infinite depth. The fluid is inviscid, and the flow in the porous medium is governed by Darcy's law. It is found that the exchange of fluid between the porous layer and the overlying fluid causes the waves to become spatially attenuated.
The present paper takes a new look at this problem, assuming that the porous layer has a finite depth; see also [20]. The main emphasis of the present investigation will be on the mean particle drift in the porous medium. For this purpose, we apply a Lagrangian description of fluid motion in the fluid layer and the porous bed. The novelty here is that the pressure gradient in Darcy's law becomes nonlinear in the Lagrangian description, yielding a direct wave-induced drift in the porous medium. In particular, for spatially damped waves, we consider the effect of radiation stress [12] at the surface on the mean motion. The Lagrangian mean surface gradient must be balanced by the gradient of the radiation stress in this case; see e.g. Phillips [16] for a discussion of this problem in Eulerian terms. Then, the Lagrangian mean drift in the fluid layer is given by the Stokes drift [19]. In the longwave approximation, the horizontal mean pressure gradient is continuous into the porous bed. Accordingly, the gradient of the radiation stress component in long waves induces a mean horizontal drift in the porous layer, in addition to the small Stokes drift.
We also point out that the exchange of mass between the upper fluid and the porous medium may cause temporal attenuation of the waves. This follows directly from the dispersion relation (52) in the "Appendix". Again, through nonlinearity, this leads to a small horizontal Stokes drift in the porous layer, but now damped in time.
The present problem has some interesting parallels to processes in the natural environment. For example, Luthar et al. [13] have found mean currents induced by surface waves in a seagrass meadow. This is of course a more complicated problem, since the fibrous porous matrix is not fixed but is moving in the oscillatory current. A closer similarity to the present problem is the investigation by Webber and Huppert [20,21], who calculate the wave-induced Stokes drift in corals reefs.
Another important environmental aspect of the present analysis is the application to the transport of microplastics in the porous oceanic bottom layer. Microplastics are defined as plastic fragments less than 5 mm in length. It originates from surface plastic debris and finally ends up at the ocean bottom. In natural pore systems, like coral reefs, the permeability typically is K = 5 × 10 −7 m 2 [20,21]. This corresponds to pore dimensions of order centimeters, and we realize that microplastics at the bottom may penetrate the bottom layer due to the oscillatory surface wave field. When inside the porous bottom layer, these plastic fragments may spread horizontally by the wave-induced Lagrangian mean current.
Apart from reefs and beaches, shallow banks in the open ocean are locations where surface waves may induce oscillatory flow at the ocean bottom. For example, on Spitsbergenbanken in the Barents Sea, the minimum water depths are 26-53 m [2], so the sea bottom is well in reach of long ocean swells. Spitsbergenbanken has a porous bottom layer of depth 5-20 m of loose carbonate material consisting of a mixture of shells, stones and gravel [24]. The permeability is typically K = 2.62 × 10 −7 m 2 [14], corresponding to an average grain diameter of about 2 cm. A considerable accumulation of floating debris (mostly plastic) is observed at the nearby Tromsøflaket [3]; see Fig. 1.
This makes Tromsøflaket a potential source of microplastics in the western Barents Sea. The proximity to the shallow Spitsbergenbanken, and the direction of the ocean currents, could lead to the accumulation of microplastics in the porous bottom layer of Spitsbergenbanken. The wave-induced Lagrangian transport would then, depending on the direction of the swell, extend laterally the contaminated region. We discuss this problem in some detail at the end of Sect. 7.

Mathematical formulation
We study surface gravity waves in a horizontal fluid layer of undisturbed depth H 1 overlying a fluid-saturated porous layer of thickness H 2 . The fluid is incompressible and homogeneous in both layers and has density . In the fluid, we neglect the effect of friction, while in the porous layer the motion is governed by Darcy's law [1]. We take that the X-axis is horizontal and situated at the interface between the porous layer and the fluid layer, while the Z-axis is vertical, and positive upwards. The pressure is P , and upper and lower layer variables are denoted by subscripts 1 and 2, respectively. The free surface of the upper layer is given by Here we assume that the pressure vanishes, i.e. P 1 = 0 . The impermeable bottom of the porous layer is given by Z 2 = −H 2 . This configuration is sketched in Fig. 2.
We assume that the surface wave in the upper layer has a wavelength that is much larger than the thickness of both layers, i.e. ≫ H 1 , H 2 . Then we can assume that the pressure is hydrostatic in the vertical. Hence, in both layers where g is the acceleration due to gravity. Utilizing that the pressure is continuous at the interface between the upper and lower layer, where Z 1 = Z 2 = (X, t) , we find by integration in the vertical that In the horizontal direction we have for the fluid layer where subscripts denote partial differentiation. In the porous medium, using Darcy's law [1], the momentum equation has the form where is the kinematic viscosity of the fluid in the pores, and K H is the horizontal component of the permeability of the porous medium. We note right away that the vertical component of the permeability K V does not enter the problem because of the shallow-water assumption.
For the mathematical analysis, we apply Lagrangian coordinates. Let a fluid particle (a, c) initially have coordinates X 0 , Z 0 . Its position (X, Z) at later times will then be a function of a, c and time t . Velocity components and accelerations are given by X t , Z t and X tt , Z tt , respectively. It is trivial to express Eqs. (4) and (5) in terms of the independent Lagrangian spatial coordinates (a, c) ; see [9]. These equations now become We note that because of the shallow-water assumption, the horizontal displacements X 1 , X 2 are not functions of the Lagrangian independent vertical variable c . Furthermore, it is interesting to notice that the governing Eq. (7) in the porous medium becomes nonlinear in the Lagrangian formulation. Even though (a, c) is not the initial particle position, it is convenient to write the displacements of the form (e.g., Pierson [17]) where the deviations ( x 1,2 , z 1,2 ) can be expanded in series after a small parameter proportional to the wave amplitude [23]. Equations (6) and (7) can thus be written The conservation of mass (here volume) leads to [9], From Eqs. (8), (11)

becomes in each layer
In the Lagrangian formulation, the boundary of the upper layer is given by c = H 1 . This is a freely moving material surface, along which the pressure must be constant. At the horizontal bottom boundary of the porous bed, given by c = −H 2 , the vertical velocity must vanish. Particles originally situated at the horizontal separation line between the fluid layer and the porous bed are given by c = 0 . After the perturbation by the surface wave, these particles will still be found at c = 0 , irrespective of whether they are located in the fluid or in the porous layer. The complicating factor is that the Darcy velocity in the porous medium is defined as an average, which makes the problem of mass and momentum transfer far from trivial; see the comprehensive discussion by Lācis et al. [8]. We here follow the approach in [7,10], and take that normal component of the velocity and the normal stress (here the pressure) are continuous at the interface, which characterizes a material surface.
Although Eqs. (9) and (10), together with the continuity Eq. (12), are nonlinear, it is seen that the hydrostatic approximation simplifies the problem considerably. In addition, as noted before, only the permeability in the horizontal direction enters the problem. This is an advantage since it appears that the permeability in nature is often anisotropic.

The Lagrangian mean drift
The linear problem is trivial, and the details have been deferred to the "Appendix". Here we consider the nonlinear problem to second order for the Lagrangian mean quantities. The mean is defined by averaging over the wave cycle, and this process is denoted by an overbar. From Eq. (9), we obtain to second order in the fluid layer where real parts are used for the linear variables and the Lagrangian mean drift in the fluid is defined as u (L) 1 = x 1t . Longuet-Higgins [11] derived the Stokes drift u (S) 1 to second order in Eulerian variables. In Lagrangian terms it can be written where we in the last expression have utilized that the flow is irrotational, i.e. x 1c =z 1a . It should be pointed out that for irrotational flow, Eames and McIntyre [5] have given an exact expression for the Stokes drift (not limited to second order in wave steepness).
In the Lagrangian formulation the Stokes drift is generally obtained from the governing equations by allowing for a small viscosity (whatever small, as long as it is non-zero). This singular behavior is discussed in the review [23]. In the present inviscid case, the vorticity vanishes in the fluid layer. By expressing this fact in Lagrangian terms, the Stokes drift can be calculated for irrotational waves [4,22].
We define the mean quantity Q by In the shallow-water approximation, the last term on the right-hand side is negligible, and we notice that Q appears in Eq. (13). Generally, by simple manipulations We note from Eq. (14) that the first term in Eq. (16) is the time derivative of the Stokes drift. The second term is the spatial gradient of the kinetic energy per unit depth and density for long gravity waves. Since we have equal partition of kinetic and potential energy in nonrotating waves, we can write for the total mechanical energy E, We finally utilize that wave energy E and wave mean momentum M are connected by (Phillips [16]). Here U (S) 1 is the vertically integrated Stokes drift (the Stokes flux) in the fluid layer, and C 0 = ∕k = gH 1 1∕2 is the phase speed (and the group velocity) of long waves. Since Eqs. (18) and (19) also are valid in the long-wave approximation, we obtain by insertion into Eq. (13) that This equation for the Lagrangian drift in high-frequency long gravity waves is also found in [23], (Eq. 114), when the effects of viscosity and the Earth's rotation are disregarded.
In the porous layer, we obtain from Eq. (10) to second order where u (L) 2 = x 2t is the horizontal Lagrangian mean drift in the porous bed. We observe from Eq. (14) that the first term on the right-hand side of Eq. (21) is just the Stokes drift for long waves in the porous bed, and we can write We thus see that the surface mean gradient in the fluid layer yields an additional drift contribution in the porous bed. It will be demonstrated in Sect. 5 that the surface gradient is proportional to the tensor component of radiation stress per unit density in shallow-water waves when the mean flow is steady.

Drift in temporally damped waves
When = 0 in Eq. (42), the drift problem becomes particularly simple. We can then assume that the mean variables do not depend on the horizontal coordinate a . Hence, from Eq. (20) we find where we have inserted from the "Appendix" for the real part of the variables to obtain u (S) 1 from Eq. (14). We note that the horizontal Lagrangian mean flow in this case is just the Stokes drift.
In the porous layer Eq. (22) reduces to where the small parameter R is defined by Eq. (55). Again, the Lagrangian mean flow in the porous bed equals the Stokes drift. We notice that the Stokes drift here is a factor R 2 smaller than the Stokes drift in the fluid layer.
We also remark that since x ta = 0 in both layers, we obtain from the continuity Eq. (12) that z tc = 0 . With a rigid bottom at c = −H 2 , and continuity of mean velocities at c = 0 , this means that z t = 0 in both layers. In the next section where we consider spatially damped waves, we shall see that this change and we have vertical mean drift velocities in the porous bed as well as in the fluid layer.

Drift in spatially damped waves
For spatially damped waves, i.e., = 0 in Eq. (43), we have that a ≠ 0 . A steady balance in Eq. (20) then requires Following Longuet-Higgins and Stewart [12], we can write the tensor component S (22) of radiation stress per unit density in shallow-water waves along the 1-axis as where we have used Eq. (19). Hence, from Eq. (25) we find This means that the mean pressure gradient at the surface due to the presence of spatially damped waves must be balanced by the radiation stress gradient per unit depth in order to achieve steady horizontal mean motion. In fact, this balance must hold also in the case when the wave amplitude varies slowly in time.
Again, the Lagrangian horizontal mean flow in the fluid layer is just the Stokes drift, exactly as it was for temporally damped waves; see Eq. (23).

Now the Stokes flux becomes
Hence, from Eq. (25) we obtain for the mean surface elevation that where A in this case is a slowly varying function of time, as will be demonstrated in Sect. 6. Using the spatial damping rate given by Eq. (56), we can write Eq. (22) in the porous bed, when we insert from Eq. (25), We remark that the effect of the radiation stress on the horizontal Lagrangian drift in the porous bed (last term in the brackets) will dominate when H 2 > 2H 1 .

Vertical mean drift for spatially damped waves
Continuity in the porous layer yields x 2ta = −z 2tc from Eq. (12). By integrating in the vertical, and using that z 2t c = −H 2 = 0 , we find where w (L) 2 = z 2t is the vertical Lagrangian mean drift in the porous bed. We note that this is very small. From Eq. (12) we find in the fluid layer to second order

3
From Eq. (32) we note that the first term on the right-hand side of Eq. (36) is of O R 3 compared to the second term, which is of O(R) from Eq. (28). Hence, the first term can be neglected. We may then write Inserting from Eqs. (28) and (30) into Eq. (37), we find Accordingly, in the case of spatially damped waves the mean free surface has the form Here A 0 is the initial wave amplitude at a = 0 , where the surface wave enters our domain. We thus see that the free surface in this case has a small upward motion in addition to the small horizontal attenuation.

Discussion and concluding remarks
It is worth noting that if we require a steady position of the free surface, i.e. This means that the volume fluxes for the drift in each layer would be equal in magnitude and oppositely directed. However, this is highly unrealistic. The mean flow in a porous medium must be much smaller than in the fluid above, since the Reynolds number based on the pore diameter must be of order 1 or less for Darcy's law to be valid [1,15]. In fact, the vertical drift velocity at the interface between the porous bed and the fluid must depend directly on the pore diameter through the permeability; see Eq. (32). When the pore diameter approaches zero, and hence K H → 0 , the interface becomes impermeable, which requires w (L) 1 (c = 0) = w (L) 2 (c = 0) = 0 . This shows that the free surface position cannot be stationary when we consider drift due to spatially damped waves in this problem.
Assuming that marine microplastic particles are non-inertial and passive, our results in Sect. 5 may be applied to Spitsbergenbanken, which has a lateral dimension of about 100 km; see Fig. 1. We take H 1 = 30 m for the ocean layer, and H 2 = 10 m , K H = 2.62 × 10 −7 m 2 and = 10 −6 m 2 s −1 for the porous bed. Furthermore, we consider a long swell of period 16 s (wavelength 400 m) and typical amplitude 1 m, arriving at the bank from the Norwegian Sea. From (28) and (31), this swell will for a period of 24 h cause a net horizontal drift of particles above the bed of about 1200 m. In the porous bed, the corresponding drift distance will be about 15 m. From the spatial damping rate (56) such waves will typically attenuate over a distance of 4 km, so swell due to repeated storm events in the Norwegian Sea will cause accumulation of microplastics at the western part of the bank.
To the authors' knowledge, the application of a Lagrangian description to flow in porous media is not common. It yields directly the mean particle drift in oscillatory wave motion. Of particular interest is the application of the long-wave approximation, which is relevant for a considerable range of physical and geophysical applications. As pointed out before, a consequence of this approximation is that the vertical component of the permeability of the porous medium does not enter the problem. This simplifies the analysis of the drift in the porous bed considerably. In addition, our results are valid even if the medium is anisotropic.
From the continuity equation in the porous bed: Integrating in the vertical, and applying that z 2 c = −H 2 = 0 , we find Finally, at the interface we must have that z 1 (c = 0) =z 2 (c = 0) (see the discussion at the end of Sect. 2), which yields the complex dispersion relation:

Temporally damped waves
We first take that = 0 . Then, to lowest order from the real and imaginary part of Eq. (52), respectively: and where we have defined the small dimensionless parameter R by According to Reid and Kajiura [18], R is the fundamental parameter of the porous problem.

Spatially damped waves
Now we take that = 0 in Eq. (52). Again, the frequency to lowest order becomes as before; see Eq. (53), while the spatial damping rates becomes It is easily seen that for the next approximation we have Since long waves are nondispersive, i.e. C g = C 0 = gH 1 1∕2 , we readily find from Eqs.   (57) 2 = gH 1 k 2 + O R 2 .
as first shown by Gaster [6]. Finally, real values of Eqs. (45), (47), (49) and (51) are used when we calculate the forcing terms for the nonlinear mean solutions.