1.

## Introduction

In a pioneering paper (Reid and Kajiura, 1957), 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, is investigated. The fluid is inviscid, and the flow in the porous medium is governed by Darcy’s law (Bear, 1972). It is found that the exchange of fluid between the porous bed and the overlying fluid causes the waves to become spatially attenuated.

In the present paper we consider internal gravity waves in the fluid layer, which is situated above a permeable bottom layer of finite thickness. The waves become attenuated due to the interaction with the lower layer. In the nonlinear analysis we consider temporally damped internal waves (real wave number and complex frequency) as well as spatially damped waves (real frequency and complex wave number). The interesting physics in this problem is that although the fluid layer is inviscid, the interaction with the frictional permeable bed causes damping in the inviscid fluid through vertical motion at the common boundary; see e.g. Reid and Kajiura (1957) for surface waves. This phenomenon is of course well-known in rotating fluids where the horizontal divergence in a frictional Ekman layer may change the vorticity in the interior inviscid fluid by Ekman pumping, leading to damping in time; see for example the classic paper by Charney and Eliassen (1949). However, we do not consider the effect of the Earth’s rotation, so here the permeable layer acts a sink of wave energy.

The main emphasis of the present investigation is on the mean particle drift in the system. For this purpose we apply a Lagrangian description of the wave motion, and of the motion in the permeable layer. The advantage here is that the pressure gradient becomes nonlinear in the Lagrangian formulation, yielding directly the wave-forcing terms in the equation for the mean drift in the permeable bed.

Internal waves are generated by phenomena like flow over topography, convective storms, imbalance of large-scale circulation, gravity currents, river plumes (Sutherland, 2010). Existing literature on various aspects of internal gravity waves is quite extensive. For a classic account in the ocean; see Phillips (1976), while Sutherland (2010) covers both the atmosphere and the ocean. For a review of the wave-mean flow interaction in internal waves using generalized Lagrangian-mean theory, the reader is referred to Bühler (2014).

A pioneering paper on the Stokes drift in internal coastal Kelvin waves with constant Brunt-Väisälä frequency $N$ is that of Wunsch (1973). Later Al-Zanaidi and Dore (1976) calculated the Stokes drift in internal waves in the case of a thin thermocline with constant $N.$ The Stokes drift results were generalized to arbitrary stable stratification by Weber et al. (2014), with application to equatorial internal Kelvin waves.

In Lagrangian form nonlinear internal waves with zero mean drift was first studied by Sanderson (1985), while in Weber (2019a) it was pointed out that internal waves with no mean drift (Gerstner-type waves) and waves with a mean forward drift (Stokes-type waves) were both possible solutions in the Lagrangian formulation.

The novel approach in this paper is that the internal wave layer is bounded below by a permeable bed, which is modelled as a macroscopic porous medium of finite thickness. As pointed out by Webber and Huppert (2020), surface wave motion over a porous (coral) bottom leads to a slow Stokes drift in the coral. In a similar way we show that the presence of internal waves above a permeable bed leads to a slow Lagrangian mean drift in the interstitial fluid. The present problem has some interesting parallels to other processes in the natural environment, such as wave-induced flow in aquatic vegetation. For example, Luhar et al. (2010) 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 rigid, but is moving in the oscillatory current.

2.

## Mathematical formulation

We study internal gravity waves in a horizontal fluid layer of undisturbed depth ${H}_{1}$ overlying a fluid-saturated permeable bed of thickness ${H}_{2};$ see Fig. 1,

Fig. 1.

A diagram showing the configuration, with internal waves above a permeable bottom layer.

The fluid in the upper layer is incompressible and stably stratified, while the interstitial fluid in the permeable bed has constant density. We assume inviscid motion in the fluid layer, while in the bottom layer the flow is dominated by friction. The effect of the Earth’s rotation is neglected in this problem (we comment on that in Section 6).

For the mathematical analysis we apply Lagrangian coordinates. Let a fluid particle $\left(a,c\right)$ initially have coordinates $\left({X}_{0},{Z}_{0}\right).$ Its position $\left(X,Z\right)$ at later times will then be a function of $a,c$ and time $t.$ Here the $X$-axis is horizontal, and situated at the interface between the permeable layer and the fluid layer, while the $Z$-axis is vertical, and positive upwards. Velocity components and accelerations are given by $\left({X}_{t},{Z}_{t}\right)$ and $\left({X}_{tt},{Z}_{tt}\right),$ respectively, where subscripts denote partial differentiation. The pressure is $P,$ and the density is $\rho .$ Upper and lower layer variables are referred to by subscripts 1 and 2, respectively. Furthermore, in the Lagrangian formulation the impermeable bottom boundary of the lower layer is given by $c=-{H}_{2},$ the permeable interface between the layers by $c=0,$ and the upper boundary of the fluid layer by $c={H}_{1}.$

We assume that the internal wave in the fluid layer has a wavelength $\lambda$ that is much larger than the thickness of both layers, i.e. $\lambda \gg {H}_{1},{H}_{2}.$ Then we can assume that the pressure is hydrostatic in the vertical. In the fluid layer we thus have for the momentum balance

((2.1))
${X}_{1tt}=-{P}_{1X}/{\rho }_{1},$
((2.2))
${0=-P}_{1Z}/{\rho }_{1}-g,$
where $g$ is the acceleration due to gravity.

The transformation of $X,Z$ from independent variables in the Eulerian description to dependent Lagrangian variables $X\left(a,c,t\right),Z\left(a,c,t\right)$ is trivial; see Lamb (1932). The conservation of mass can be written

((2.3))
${\rho }_{1}J\left({X}_{1},{Z}_{1}\right)={\rho }_{10}J\left({X}_{10},{Z}_{10}\right),$
where subscript 10 denotes initial values in the upper layer, and $J\left(F,G\right)\equiv {F}_{a}{G}_{c}-{F}_{c}{G}_{a}$ is the two-dimensional Jacobian. We assume that the fluid is incompressible, that is
((2.4))
${\rho }_{1}\left(a,c,t\right)={\rho }_{10}\left(a,c,0\right),$
so (2.3) reduces to
((2.5))
$J\left({X}_{1},{Z}_{1}\right)=J\left({X}_{10},{Z}_{10}\right).$

Utilizing (2.4), (2.1) and (2.2) now become in Lagrangian terms

((2.6))
${X}_{1tt}{X}_{1a}+g{Z}_{1a}=-{P}_{1a}/{\rho }_{10},$
((2.7))
${X}_{1tt}{X}_{1c}+g{Z}_{1c}=-{P}_{1c}/{\rho }_{10}.$

In this problem we assume that ${\rho }_{10}\left(c=0\right)={\rho }_{2},$ where ${\rho }_{2}$ is a constant. Using the bottom layer density as a reference density, we apply the Boussinesq approximation. Introducing the Brunt-Väisälä frequency $N$ by

((2.8))
${N}^{2}=-\left(g/{\rho }_{2}\right)\left(d{\rho }_{10}/dc\right),$
the momentum balance in the fluid layer becomes
((2.9))
${X}_{1tt}{X}_{1a}+g\left({\rho }_{10}/{\rho }_{2}\right){Z}_{1a}=-{P}_{1a}/{\rho }_{2},$
((2.10))
${X}_{1tt}{X}_{1c}+g{\left({\rho }_{10}/{\rho }_{2}\right)Z}_{1c}=-{P}_{1c}/{\rho }_{2}.$

We assume that $N$ is constant in this problem. Hence, in (2.9) and (2.10)

((2.11))
${g\left(\rho }_{10}/{\rho }_{2}\right)=g-c{N}^{2}.$

A treatise in Lagrangian coordinates of internal waves with constant $N,$ was first done by Sanderson (1985), without assuming the long-wave approximation.

For laminar slow fluid motion in a porous medium with narrow pores (the microscopic scale), like sand or sandstone, Darcy’s law applies (Bear, 1972). For a macroscopic permeable medium with a much coarser solid matrix, we take the fluid motion in the voids is turbulent. Actually, the modelling of turbulent flows in porous media is an open-end field; see the comprehensive review by Wood et al. (2020), and references therein. In the present study we apply a rather simple flow model, assuming that the Reynolds-averaged Navier-Stokes (RANS) equations are valid in the voids. By averaging the RANS equations over a representative elementary volume of the porous medium domain, a macroscale equation for the momentum balance is obtained, with an effective permeability $K$ and an overall eddy viscosity ${\nu }_{2};$ see for example Masuoka and Takatsu (1996) for details on how to define the eddy viscosity.

In the present paper we consider a thin horizontal macroscopic porous layer, and assume that the balance of forces in the vertical direction is hydrostatic. The motion in the permeable layer is forced by long internal waves with typical frequency $\omega$ and wavelength $\lambda .$ Then we find that if the parameter $\omega K/{\nu }_{2}$ is small, the local acceleration in the macroscopic momentum equation can be neglected compared to the macroscopic Darcy friction term. Similarly, if the typical horizontal velocity in the macroscopic medium is $U,$ the convective acceleration can be neglected if $UK/\left(\lambda {\nu }_{2}\right)$ is small. Finally, the Laplacian eddy viscosity term is negligible if $K/{\lambda }^{2}$ is small. In summary, our application of Darcy’s law in a macroscopic porous medium is valid when

((2.12))
$\begin{array}{c}\omega K/{\nu }_{2}\ll 1,\\ UK/\left(\lambda {\nu }_{2}\right)\ll 1,\\ K/{\lambda }^{2}\ll 1.\end{array}}$

In the permeable layer we have for the momentum balance, when (2.12) is fulfilled,

((2.13))
$0=-{P}_{2X}/{\rho }_{2}-\left({\nu }_{2}/K\right){X}_{2t},$
((2.14))
${0=-P}_{2Z}/{\rho }_{2}-g,$
where (2.13) is the 0-equation by Masuoka and Takatsu (1996). Since ${\rho }_{2}$ is constant, we get from (2.14) that
((2.15))
${P}_{2}/{\rho }_{2}=-g{Z}_{2}+{Q}_{2},$
where ${Q}_{2}$ is independent of ${Z}_{2}.$ Hence, by transforming into Lagrangian coordinates, (2.13) and (2.14) become
((2.16))
$0=-{Q}_{2a}-\left({\nu }_{2}/K\right){X}_{2t}{X}_{2a},$
((2.17))
$0=-{Q}_{2c}-\left({\nu }_{2}/K\right){X}_{2t}{X}_{2c}.$

Here, since ${Q}_{2}$ is independent of ${Z}_{2},$ and thereby ${Q}_{2c}=0,$ the last equation implies that ${X}_{2}={X}_{2}\left(a,t\right).$

Even though $\left(a,c\right)$ is not the initial particle position, it is convenient to write the displacements and pressure as (Pierson, 1962),

((2.18))
$\begin{array}{c}{X}_{1}=a+{x}_{1}\left(a,c,t\right)\\ {X}_{2}=a+{x}_{2}\left(a,t\right)\\ {\mathrm{}Z}_{1,2}=c+{z}_{1,2}\left(a,c,t\right)\\ {P}_{1}={P}_{10}\left(c\right)+{p}_{1}\left(a,c,t\right)\\ {P}_{2}=-{\rho }_{2}gc+{p}_{2}\left(a,c,t\right)\end{array}},$
where $d{P}_{10}/dc=-{\rho }_{10}g.$ Applying (2.11), we find, when we put an insignificant constant equal to zero
((2.19))
${P}_{10}=-{\rho }_{2}g\left[c-{N}^{2}{c}^{2}/\left(2g\right)\right].$

The deviations (${x}_{1,2},{z}_{1,2},{p}_{1,2}\right)$ in (2.18) can subsequently be expanded in series after a small parameter proportional to the wave steepness; see for example Weber (2019b). Using (2.18), we can write (2.9) and (2.10) in the fluid as

((2.20))
${x}_{1tt}{\left(1+x}_{1a}\right)+g\left({\rho }_{10}/{\rho }_{2}\right){z}_{1a}=-{p}_{1a}/{\rho }_{2},$
((2.21))
${x}_{1tt}{x}_{1c}+g{\left({\rho }_{10}/{\rho }_{2}\right)\left(1+z}_{1c}\right)=-{p}_{1c}/{\rho }_{2}.$

In the permeable layer, the momentum equation becomes from (2.16)

((2.22))
$\left({\nu }_{2}/K\right){x}_{2t}\left(1+{x}_{2a}\right)=-{Q}_{2a}.$

The conservation of mass (2.5) is the same in both layers. Hence

((2.23))
${x}_{1,2ta}+{z}_{1,2tc}+{\left({x}_{1,2a}{z}_{1,2c}\right)}_{t}-{\left({x}_{1,2c}{z}_{1,2a}\right)}_{t}=0.$

Although our governing equations (2.20)-(2.23) are nonlinear, it is seen that the hydrostatic approximation simplifies the problem considerably.

3.

## Linear analysis

We linearize our governing equations (2.20)-(2.23), and denote the linearized variables by a tilde. From the curl of the linearized versions of (2.20) and (2.21), we obtain

((3.1))
${\stackrel{\sim }{x}}_{1ttc}-{N}^{2}{\stackrel{\sim }{z}}_{1a}=0.$

Applying (2.23), from which ${\stackrel{\sim }{x}}_{1a}=-{\stackrel{\sim }{z}}_{1c},$ we find

((3.2))
${\stackrel{\sim }{z}}_{1t\mathit{tcc}}+{N}^{2}{\stackrel{\sim }{z}}_{1aa}=0.$

We assume that the variables separate, and write

((3.3))
${\stackrel{\sim }{z}}_{1}=\phi \left(c\right){e}^{i\theta },$
where $\theta =\kappa a+i\delta t$ is the phase function. Here we have defined
((3.4))
$\kappa =k+i\alpha ,$
((3.5))
$\delta =\beta +i\omega .$

In (3.4)-(3.5) $k$ is the real wave number and $\omega$ the real wave frequency, while $\alpha$ and $\beta$ are the real small spatial and temporal attenuation coefficients, respectively. Inserting into (3.2), we find

((3.6))
${\phi }^{\text{″}}+{\gamma }^{2}\phi =0,$
where primes denote derivative with respect to $c.$ In (3.6) we have defined
((3.7))
$\gamma =i\kappa N/\delta .$

To eliminate the small barotropic response in the fluid layer, we use the rigid-lid approximation for the upper boundary; see Gill and Clarke (1974). Since then $\phi \left(c={H}_{1}\right)=0,$ the solution to (3.6) can be written

((3.8))
$\phi =A\text{sin}\gamma \left(c-{H}_{1}\right),$
where $A$ is the internal wave amplitude. The linear displacements thus become
((3.9))
${\stackrel{\sim }{x}}_{1}=\left(i\gamma A/\kappa \right)\text{cos}\gamma \left(c-{H}_{1}\right){e}^{i\theta },$
((3.10))
${\stackrel{\sim }{z}}_{1}=A\text{sin}\gamma \left(c-{H}_{1}\right){e}^{i\theta }.$

In the permeable bed the linearized version of (2.22) is

((3.11))
${\stackrel{\sim }{x}}_{2t}=-\left(K/{\nu }_{2}\right){\stackrel{\sim }{Q}}_{2a}.$

Introducing ${\stackrel{\sim }{Q}}_{2}={q}_{2} \text{exp} i\theta ,$ and applying that the lower boundary of the bottom layer is impermeable, that is ${\stackrel{\sim }{z}}_{2}\left(c=-{H}_{2}\right)=0,$ we find by using the (3.11), together with the continuity equation (2.23), that

((3.12))
${\stackrel{\sim }{x}}_{2}=\left[i\kappa K{q}_{2}/\left({\nu }_{2}\delta \right)\right]{e}^{i\theta },$
((3.13))
${\stackrel{\sim }{z}}_{2}=\left[{\kappa }^{2}K{q}_{2}/\left({\nu }_{2}\delta \right)\right]\left(c+{H}_{2}\right){e}^{i\theta }.$

3.1.

### Conditions at the interface and the dispersion relation

In the present problem we consider an inviscid stratified fluid layer overlaying a permeable medium saturated with the same fluid, but of constant density. The surface of the upper layer at $c={H}_{1},$ and the impermeable bottom of the porous layer at $c=-{H}_{2},$ are material surfaces. On the other hand, the horizontal separation between the fluid layer and the porous bed at $c=0$ is an immaterial surface. Here particles from the upper layer may penetrate into the lower layer, and vice versa. 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. (2020). We here follow the approach by Levy and Sanchez-Palencia (1975), and take that normal component of the velocity and the normal stress (here the pressure) are continuous at the interface; see also Kandel and Pascal (2013). However, these conditions are the same as if the interface was a material surface. Hence, we can use our Lagrangian formulation where particles characterized by $c\ge 0$ belong to the fluid layer, and particles with $c\le 0$ are confined to the porous bed. Accordingly, we take that ${\stackrel{\sim }{z}}_{1}\left(c=0\right)={\stackrel{\sim }{z}}_{2}\left(c=0\right).$ We then obtain from (3.10) and (3.13)

((3.14))
$A=-{\kappa }^{2}{H}_{2}K{q}_{2}/\left({\nu }_{2}\delta \text{sin}\gamma {H}_{1}\right),$
which relates the internal wave amplitude to the amplitude ${q}_{2}$ of the fluctuating pressure gradient in the permeable layer; see (3.11).

With a hydrostatic pressure distribution in both layers, the pressure must be continuous at the interface, or ${P}_{1}\left(c=0\right)={P}_{2}\left(c=0\right).$ Utilizing (2.18), this means that ${p}_{1}\left(c=0\right)={p}_{2}\left(c=0\right).$ But since ${\rho }_{10}\left(c=0\right)={\rho }_{2},$ we realize that

((3.15))
${\stackrel{\sim }{p}}_{1}+{\rho }_{2}g{\stackrel{\sim }{z}}_{1}={\stackrel{\sim }{p}}_{2}+{\rho }_{2}g{\stackrel{\sim }{z}}_{2},\mathrm{}c=0.$

Hence, using (3.15), we find from (2.20) and (2.22)

((3.16))
${\stackrel{\sim }{x}}_{1tt}\left(c=0\right)=-{\stackrel{\sim }{Q}}_{2a}=-{i\kappa q}_{2}{e}^{i\theta }.$

Inserting from (3.9) and (3.14) into (3.16), we finally obtain the complex dispersion relation for this problem

((3.17))
$i\kappa {H}_{2}KN/{\nu }_{2}=\text{tan}\left(i\kappa {H}_{1}N/\delta \right).$

3.2.

### Spatially damped waves

For spatially damped waves we take that $\beta =0$ in (3.5), and introduce the small dimensionless parameter $M$ by

((3.18))
$M=k{H}_{2}KN/{\nu }_{2}\ll 1$

Then, from (3.17)

((3.19))
$\left(i-\alpha /k\right)M=\text{tan}\left(D\left(1+i\alpha /k\right)\right),$
where we have defined
((3.20))
$D=k{H}_{1}N/\omega .$

In our analysis we assume that the waves are slowly damped over one wavelength, that is

((3.21))
$\alpha /k\ll 1.$

Utilizing (3.18) and (3.21), a simple series expansion yields to lowest order in $M$ and $\alpha /k$

((3.22))
$iM=\text{tan}D+i\alpha D/{\left(k\text{cos}}^{2}D\right).$

The real part of this equation requires that

((3.23))
$\text{tan} D=0,$
which is fulfilled when $D=n\pi ,n=1,2,3,\dots$ In this case, with given frequency, the permissible wave numbers obtained from (3.20) are
((3.24))
${k}_{n}=n\pi \omega /\left(N{H}_{1}\right).$

From the imaginary part of (3.22) we obtain

((3.25))
$\alpha D/{\left(k\text{cos}}^{2}D\right)=M.$

Hence, using that $D=n\pi ,$ we find

((3.26))
${\alpha }_{n}/{k}_{n}=M/\left(n\pi \right)={\left(H}_{1}/{H}_{2}\right)R,$
where $R,$ as mention in the abstract, is a small parameter given by
((3.27))
$R=\omega K/{\nu }_{2}.$

For surface waves over an infinitely deep porous layer, Reid and Kajiura (1957) state that $R$ (with molecular $\nu$ and a different notation) is the fundamental parameter of the porous problem. We note that for a macroscopic permeable bed of finite depth, this is also the case for internal wave forcing, but now modified by the depth ratio ${H}_{1}/{H}_{2}.$

3.3.

### Temporally damped waves

For temporally damped waves we take that $\alpha =0$ in (3.4). Assuming $\beta /\omega \ll 1,$ we then obtain from (3.17)

((3.28))
$iM=\text{tan}D+i\beta D/{\left(\omega \text{cos}}^{2}D\right).$

Again, the real part is satisfied for $D=n\pi ,n=1,2,3,.$ In this case for given $k,$ we obtain for the permissible frequencies

((3.29))
${\omega }_{n}=k{H}_{1}N/\left(n\pi \right).$

The imaginary part of (3.28) reduces to

((3.30))
$\beta D/{\left(\omega \text{cos}}^{2}D\right)=M.$

Hence,

((3.31))
${\beta }_{n}/{\omega }_{n}=M/\left(n\pi \right)={\left(H}_{1}/{H}_{2}\right){R}_{n},$
where now ${R}_{n}={\omega }_{n}K/{\nu }_{2}.$

Long internal waves are non-dispersive, so the group velocity is ${C}_{g}=\omega /k.$ Thus, we note from (3.26) and (3.31) that for each wave mode

((3.32))
$\beta =\left(\omega /k\right)\alpha ={C}_{g}\alpha ,$
as first shown by Gaster (1962).

4.

## Lagrangian mean drift induced by temporally damped waves

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 over-bar. In the mean drift calculations we use real values of the linear wave variables.

Traditionally, the Stokes drift (Stokes, 1847) to second order is obtained from the formula in Eulerian variables by Longuet-Higgins (1953). In Lagrangian terms Longuet-Higgins’ expression can be stated as

((4.1))
${u}^{\left(S\right)}=\stackrel{-}{{\stackrel{\sim }{x}\stackrel{\sim }{x}}_{ta}}+\stackrel{-}{\stackrel{\sim }{z}{\stackrel{\sim }{x}}_{tc}}.$

As pointed out by Longuet-Higgins (1953), the Lagrangian mean drift ${u}^{\left(L\right)}$ can be written as

((4.2))
${u}^{\left(L\right)}={u}^{\left(S\right)}+{u}^{\left(E\right)},$
where ${u}^{\left(E\right)}$ is the viscosity-dependent Eulerian mean velocity. We demonstrate in Section 5 that for spatially damped waves, ${u}^{\left(E\right)}$ becomes important in the permeable bed.

4.1.

### Drift in the fluid layer

We first focus on temporally damped waves, which is the simplest case. Then we can assume that the mean variables do not depend on the horizontal coordinate $a.$ Accordingly, from (2.20),

((4.3))
${u}_{1t}^{\left(L\right)}=-\stackrel{-}{{\stackrel{\sim }{x}}_{1tt}{\stackrel{\sim }{x}}_{1a}},$
where ${u}_{1}^{\left(L\right)}={\overline{x}}_{1t}.$ We note that if we consider permanent waves (constant amplitude in time and space), the phase differences on the right-hand side of (4.3) would result in ${u}_{1t}^{\left(L\right)}=0,$ or ${u}_{1}^{\left(L\right)}=F\left(c\right),$ which is correct, but inconclusive. For irrotational inviscid wave motion, one may utilize the fact that the mean vorticity vanishes identically. This yields the correct expression for ${u}_{1}^{\left(L\right)},$ which in this case is equal to the Stokes drift; see Clamond (2007) for surface waves.

However, internal waves possess vorticity, which makes the determination of ${u}_{1}^{\left(L\right)}$ from a direct Lagrangian approach more complicated. In the present problem the amplitudes of the wave variables on the right-hand side of (4.3) are exponentially damped in time. From a Lagrangian point of view, this small damping is related to a weak viscosity ${\nu }_{1}$ in the fluid. So to obtain the Lagrangian mean drift in this case, a small non-zero viscosity has to be included into the problem; see the discussion in Weber (2019b) for surface waves.

Hence, we introduce the effect of a small viscosity in the stratified case, and consider temporally damped waves. Again, since the mean variables do no depend on the horizontal coordinate $a,$ we can write the horizontal mean momentum (2.20) when viscosity is included as

((4.4))
${\overline{x}}_{1tt}+\stackrel{-}{{\stackrel{\sim }{x}}_{1tt}{\stackrel{\sim }{x}}_{1a}}-{\nu }_{1}\left(\mathrm{}\stackrel{-}{{\stackrel{\sim }{x}}_{1tcc}{\stackrel{\sim }{x}}_{1a}}+\overline{{\nabla }^{2}{x}_{1t}}\right)=0.$

The Laplacian operator ${\nabla }^{2}$ is nonlinear in the Lagrangian formulation. Utilizing the results by Pierson (1962), and assuming long waves, that is $|{\partial }^{2}/{\partial c}^{2}|\gg |{\partial }^{2}/{\partial a}^{2}|,$ (4.4) becomes to second order

((4.5))
${{\nu }_{1}u}_{1cc}^{\left(L\right)}-{u}_{1t}^{\left(L\right)}=\overline{{\stackrel{˜}{x}}_{1tt}{\stackrel{˜}{x}}_{1a}}-{\nu }_{1}\left[4\overline{{\stackrel{˜}{x}}_{1tcc}{\stackrel{˜}{x}}_{1a}}+3\overline{{\stackrel{˜}{x}}_{1tc}{\stackrel{˜}{x}}_{1ac}}\right],$
where we have utilized that ${\stackrel{\sim }{x}}_{1a}=-{\stackrel{\sim }{z}}_{1c}$ from continuity. The forcing-terms on the right-hand side of (4.5) are proportional to $\mathrm{exp}\left(-2{\beta }_{n}t\right).$ From the linear part of the viscous wave problem, it is easily shown that
((4.6))
${\beta }_{n}={\nu }_{1}{l}_{n}^{2}/2,$
were ${l}_{n}=n\pi /{H}_{1}$ is the vertical wave number (Weber, 2019a). In the present case the wave damping (3.31) is not determined by the viscosity in the fluid layer, but by the vertical motion at the interface. For consistency, the fluid viscosity in (4.5) must then be modelled from (4.6) as ${\nu }_{1}=2{\beta }_{n}/{l}_{n}^{2}.$

Inserting from the real part of (3.9), we find from (4.5) to lowest order in ${\beta }_{n}/\omega ,$

((4.7))
${\nu }_{1}{u}_{1ncc}^{\left(L\right)}-{u}_{1nt}^{\left(L\right)}=-\left(3/2\right){\nu }_{1}{C}_{n}{l}_{n}^{4}{A}_{n}^{2}{e}^{-2{\beta }_{n}t}\text{cos}2{l}_{n}\left(c-{H}_{1}\right),$
where ${C}_{n}=\omega /{k}_{n}=N{H}_{1}/\left(n\pi \right)$ is the phase speed of mode $n.$ The solution to (4.7) contains the Stokes drift as well as the Eulerian mean flow due to viscosity; see (4.2). Remark, however, by using the linear solution (3.9) as forcing term, we have neglected the Eulerian mean currents resulting from viscous effects at the boundaries. These are in effect treated as free-slip boundaries in our analysis. Utilizing that ${\nu }_{1}=2{\beta }_{n}/{l}_{n}^{2},$ direct integration of (4.7) yields the Stokes drift for mode $n,$
((4.8))
${u}_{1n}^{\left(S\right)}=\frac{1}{2}\left(n\pi N/{H}_{1}\right){A}_{n}^{2}{e}^{-2{\beta }_{n}t}\text{cos}2{l}_{n}\left(c-{H}_{1}\right).$

We note that the Stokes drift modes vary trigonometrically in the vertical (Wunsch, 1973) in such a way that the integral over the layer (the Stokes flux) vanishes identically. The last result is valid for an arbitrary stable stratification, as shown by Weber et al. (2014). We also remark that each Stokes drift mode is positive at the upper $\left(c={H}_{1}$) and lower $\left(c=0$) boundary.

We should point out here that the solution (4.8) for the Stokes drift is valid however small ${\beta }_{n}$ is, as long as it is nonzero. In the present problem ${\beta }_{n}\to 0$ when $K\to 0,$ that is when the lower permeable layer becomes solid. Finally, it is a simple exercise to show that the application of (4.1) yields the Stokes drift (4.8).

4.2.

### Drift in the permeable bottom layer

For temporally damped waves, we find to second order when averaging (2.22)

((4.9))
${u}_{2}^{\left(L\right)}=-\stackrel{-}{{\mathrm{}{\stackrel{\sim }{x}}_{2t}\stackrel{\sim }{x}}_{2a}}={-\left(\stackrel{-}{{\mathrm{}{\stackrel{\sim }{x}}_{2t}\stackrel{\sim }{x}}_{2}}\right)}_{a}+\stackrel{-}{{\mathrm{}{\stackrel{\sim }{x}}_{2}\stackrel{\sim }{x}}_{2ta}},$
where ${{u}_{2}^{\left(L\right)}=\overline{x}}_{2t}$ is the horizontal Lagrangian mean drift in the permeable bed. We observe right away that the first term on the right-hand side is zero for temporally damped waves. Furthermore, since ${\stackrel{\sim }{x}}_{2}={\stackrel{\sim }{x}}_{2}\left(a,t\right);$ see (3.12), we note from (4.1) that (4.9) is just the Stokes drift in the bottom layer,
((4.10))
${u}_{2}^{\left(L\right)}={u}_{2}^{\left(S\right)}=\stackrel{-}{{\mathrm{}{\stackrel{\sim }{x}}_{2}\stackrel{\sim }{x}}_{2ta}}.$

To calculate the forcing in (4.10), we must first relate the linear pressure gradient ${q}_{2n}$ in the bottom layer to the amplitude ${A}_{n}$ for mode $n$ in the case of temporally damped waves. Inserting into the general expression (3.14) from (3.29) and (3.31) we find that

((4.11))
${q}_{2n}={\left(-1\right)}^{n}{NC}_{n}{A}_{n}.$

Hence, from the real part of (3.12), using (4.11), the Stokes drift (4.10) for mode $n$ becomes

((4.12))
${u}_{2n}^{\left(S\right)}=\frac{1}{2}\left(n\pi N/{H}_{1}\right){R}_{n}^{2}{A}_{n}^{2}{e}^{-2{\beta }_{n}t},$
where the small parameter ${R}_{n}$ is defined by (3.31). We notice that the Stokes drift in the permeable bed is a factor ${R}_{n}^{2}$ smaller than the Stokes drift in the fluid layer at $c=0.$

We also remark that since ${\overline{x}}_{ta}=0$ in both layers, we obtain from the continuity equation (2.23) that ${\overline{z}}_{tc}=0.$ With an impermeable boundary at $c={-H}_{2},$ plus continuity of mean velocities at $c=0,$ this means that ${\overline{z}}_{t}=0$ in both layers. We demonstrate in Section 5 that this changes when we consider spatially damped waves.

5.

## Lagrangian mean drift induced by spatially damped waves

In the case of spatial damping, we have nonzero horizontal mean gradients in the system. Then, from (2.20)

((5.1))
${\overline{x}}_{1tt}=\mathrm{}{\omega }^{2}\stackrel{-}{{\mathrm{}{\stackrel{\sim }{x}}_{1}\stackrel{\sim }{x}}_{1a}}{-\left[{\overline{p}}_{1}/{\rho }_{2}+g\left({\rho }_{10}/{\rho }_{2}\right){\overline{z}}_{1}\right]}_{a},$
where we have used that for spatially damped waves ${\stackrel{\sim }{x}}_{1tt}=-{\omega }^{2}{\stackrel{\sim }{x}}_{1}.$ In fact, the first term on the right-hand side may be developed further, such that
((5.2))
${\overline{x}}_{1tt}=\mathrm{}-\alpha {\omega }^{2}\stackrel{-}{{\stackrel{\sim }{x}}_{1}^{2}}{-\left[{\overline{p}}_{1}/{\rho }_{2}+g\left({\rho }_{10}/{\rho }_{2}\right){\overline{z}}_{1}\right]}_{a}.$

We realize that in order to avoid infinite drift velocities as time increases, the right-hand side must vanish identically, that is

((5.3))
${\left[{\overline{p}}_{1}/{\rho }_{2}+g\left({\rho }_{10}/{\rho }_{2}\right){\overline{z}}_{1}\right]}_{a}=-\alpha {\omega }^{2}\stackrel{-}{{\stackrel{\sim }{x}}_{1}^{2}}.$

Just as for temporally damped waves, this means from (5.1) that the Lagrangian mean drift reduces to the Stokes drift in the upper layer, but now slowly damped in space,

((5.4))
${u}_{1n}^{\left(L\right)}={u}_{1n}^{\left(S\right)}=\frac{1}{2}\left(n\pi N/{H}_{1}\right){A}_{n}^{2}{e}^{-2{\alpha }_{n}a}\text{cos}2{l}_{n}\left(c-{H}_{1}\right).$

For spatial damping, we now obtain from (2.15) in the permeable bed

((5.5))
${\left[{\overline{p}}_{2}/{\rho }_{2}+g{\overline{z}}_{2}\right]}_{a}={\overline{Q}}_{2a},$
where the mean pressure gradient ${\overline{Q}}_{2a}$ is a forcing term for the mean drift; see (2.22). At the interface $c=0$ we must have ${\overline{p}}_{1}={\overline{p}}_{2}$ and ${\overline{z}}_{1}={\overline{z}}_{2}.$ Hence, from (5.3) and (5.5), using that ${\rho }_{10}\left(c=0\right)={\rho }_{2},$ we obtain
((5.6))
${\overline{Q}}_{2a}=-\alpha {\omega }^{2}\stackrel{-}{{\stackrel{\sim }{x}}_{1}^{2}},\mathrm{}c=0.$

Inserting to lowest order from (3.9), we find for mode $n$

((5.7))
${\overline{Q}}_{2na}={\frac{1}{2}\alpha }_{n}{N}^{2}{A}_{n}^{2}{e}^{-2{\alpha }_{n}a}$

Hence, from (2.22),

((5.8))
${u}_{2}^{\left(L\right)}=-\stackrel{-}{{\stackrel{\sim }{x}}_{2t}{\stackrel{\sim }{x}}_{2a}}\mathrm{}{-\left(K/{\nu }_{2}\right){\overline{Q}}_{2a}=u}_{2}^{\left(S\right)}-\left(K/{\nu }_{2}\right){\overline{Q}}_{2a}.$

With reference to (4.2), we note that the mean pressure gradient ${\overline{Q}}_{2a}$ induces an Eulerian mean current in the permeable bed that comes in addition to the Stokes drift. Inserting from (3.12) and (5.7), we obtain for mode $n$

((5.9))
${u}_{2n}^{\left(L\right)}=\frac{1}{2}n\pi N\left[1/{H}_{1}+1/{H}_{2}\right]{R}^{2}{A}_{n}^{2}{e}^{-2{\alpha }_{n}a}.$

We notice that if ${\mathrm{H}}_{2}<{H}_{1},$ the last term in the parenthesis (from the Eulerian mean flow) will dominate the Lagrangian drift in the permeable bed. In the suggested application in Section 6 of this theory to the slow drift of particles in the tropical rainforest, modelled as a macroscopic porous medium, we obviously have that ${H}_{2}\ll {H}_{1},$ where ${H}_{1}$ typically is the height of the troposphere. Hence, the Eulerian mean flow will dominate in this problem.

We also remark that since now ${\overline{x}}_{1,2ta}\ne 0,$ we have a slow vertical mean drift in the system. The continuity equation (2.23) in the upper layer reduces to

((5.10))
${\overline{z}}_{1tc}=-{\overline{x}}_{1ta}.$

Integrating in the vertical, and applying that the total Stokes flux in internal waves is zero (Weber et al., 2014), we obtain from (5.10) that

((5.11))
${\overline{z}}_{1t}\left({H}_{1}\right)={\overline{z}}_{1t}\left(0\right).$

By integrating ${\overline{z}}_{2tc}=-{\overline{x}}_{2ta}$ in the lower layer, and utilizing that the bottom boundary is impermeable, we find

((5.12))
${\overline{z}}_{2t}=2\alpha {\left(c+H}_{2}\right){u}_{2}^{\left(L\right)}.$

From (3.26) and (5.9), we note that this vertical mean drift is of $O\left({R}^{3}\right).$ Since we from mass conservation must have ${\overline{z}}_{1t}\left(0\right)={\overline{z}}_{2t}\left(0\right),$ we realize from (5.11) that

((5.13))
${\overline{z}}_{1t}\left({H}_{1}\right)=2\alpha {H}_{2}{u}_{2}^{\left(L\right)}.$
Hence, also the upper boundary moves vertically with a velocity of $O\left({R}^{3}\right)$ (it should be pointed out that this does not contradict the rigid-lid assumption, which was applied to the linear problem in Section 3 to eliminate the barotropic wave response). Anyway, the very small vertical mean drift discussed here is probably negligible for all practical purposes.

6.

## Discussion and final remarks

To the author’s knowledge the application of a Lagrangian description to flows in porous media is not common. This formulation has the advantage that it yields directly the mean particle drift in the medium. Of particular interest is the use of the long-wave approximation when the layer thickness is small compared to the scale of the external forcing. This assumption simplifies the drift problem considerably. Furthermore, it is relevant to geophysical applications; see the discussion below. It should be mentioned, however, that any application to nature faces the problem of estimating the macroscopic permeability $K$ and the overall eddy viscosity ${\nu }_{2}$ in the voids. But this is not peculiar to flows in porous media. Similar estimates must be done when we apply sub-grid quantities like eddy viscosity or eddy diffusivity in large-scale atmosphere or ocean modelling. When it comes to turbulent flows in porous media, Wood et al. (2020) contain references to studies suggesting various ways (measurements and/or numerical simulations) of solving the interstitial turbulent fluid problem, but this is outside the scope of the present paper.

The wave-induced drift in a permeable bottom layer has several interesting geophysical applications. As mentioned in the Introduction, surface gravity waves in the ocean may cause vertical and horizontal drift in coral reefs (Webber and Huppert 2020), which is important for maintaining a healthy reef system. Likewise, surface waves may induce a mean transport of microplastics in the porous bottom layer of shallow banks, where the bottom consists of a mixture of shells, stones and gravel (Weber and Ghaffari 2021). In the atmosphere, where the gravity waves are of the internal type, a comparable case would be the (potentially) wave-induced drift in large boulder fields (felsenmeere) found in many places of the earth. The analysis in this problem could also pertain to internal waves in the atmosphere interacting with a dense forest canopy, which by far appears to be the most interesting application.

For oceanic applications we need not worry about the additional presence of slowly varying ocean currents, since the gravity wave motion is dominating, and the mean currents are second order quantities. In the atmosphere, it is opposite. Here fluctuating wind fields dominate near the ground, and the gravity wave motion is secondary.

For the boreal forest, there exists an extensive literature on wind-driven flows in forest canopies; see the reviews by Raupach and Thom (1981), Finnigan (2000), Belcher et al. (2012), and the many references therein. The main aim in many of these studies is the modelling of the turbulent boundary layer near the ground in various complex terrains. One particularly important purpose is to quantify the $\mathrm{C}{\mathrm{O}}_{2}$ exchange between forests and the atmosphere.

In normal circumstances the action of the wind will subdue the weak mean drift induced by the wave motion in the forest canopy, as described in this paper. Exceptional cases with no or little winds occur, but they are sporadic. However, one possible exception of a more permanent character could be the tropical rainforest, which experiences very light winds. The average wind speed above the canopy is 2-3 $\mathrm{m}\mathrm{}{\mathrm{s}}^{-1},$ and often less than half of this value. Wind speeds are even smaller when measured below the tree tops. Because the rainforest consists of large stands of broad-leaved trees, any light breeze passing above the dense canopy is obstructed in the understory. For example, wind speeds near the forest floor in a Colombian jungle were typically between 1 and 5 percent of the speed recorded above the canopy (Baynton et al., 1965). Due to its vast size, the Amazonian tropical rain forest is of particular interest. Here Acevedo et al. (2004) observed median nocturnal winds of 0.75 $\mathrm{m}\mathrm{}{\mathrm{s}}^{-1}$ at the 5.7 m level, while on many nights the average wind speed fell considerably below the median value. At another measurement site in the Amazonas wind velocities were well below 0.5 $\mathrm{m}\mathrm{}{\mathrm{s}}^{-1}$ at 25 m height, both day and night (Fitzjarrald et al. 1990). In addition, the wind direction appears to be disorganized (do Couto-Santos and Luizão, 2010).

The primary production term for the turbulent kinetic energy within the mixing-layer region of a forest canopy is proportional to the vertical wind shear, while the dominant balance of forces on an air parcel is between advective deceleration, the pressure gradient, and the canopy drag (Belcher et al. 2003, 2012). With this in mind, the light winds and the weak horizontal and vertical wind shear in the tropical forest may justify the use of the lowest order momentum balance between the pressure gradient and a linear drag (Masuoka and Takatsu 1996) in modelling the rain forest as a macroscopic porous layer. Hence, the wave-induced drift in a permeable bottom layer explored in this paper could provide a model for the vertical fluctuations and the slow net transport of small substances in the tropical rainforest, like bio aerosols emitted from the ecosystem and smoke particles from deforestation fires.

With the application to the equatorial region in mind, it is a reasonable approximation to neglect the Coriolis force in the momentum equation for the upper layer in Section 2. Alternatively, we could use the beta-plane approximation. Then the waves become internal equatorial Kelvin waves, propagating eastwards and trapped within the baroclinic equatorial Rossby radius. Near the equator the result for the Stokes drift will be approximately the same as that derived in Section 4 and Section 5, if we direct the $X$ axis eastward along the equator; see the analysis by Weber et al. (2014).

Although the present theory appears to be sound from a fluid dynamic point of view, there are unfortunately no field observations to corroborate it. Generally, wave-induced drift is notoriously hard to observe, and adequate experiments may be particularly difficult, both methodically and practically, to perform in the dense rainforest. A future possible way to approach this problem theoretically and experimentally is by the wave tank set up of Luhar et al. (2010), but with suitably arranged vertical cylinders in the bottom layer. This could provide interesting results with relevance to wave-induced drift in rigid fibrous materials, like the dense rainforest.