Modelling and simulation of waves in three-layer porous media

The propagation of gravity waves in an emerged three-layer porous medium is considered in this paper. Based on the assumption that the flow can be described by Darcy’s Law, an asymptotic theory is developed for small-amplitude long waves. This leads to a weakly nonlinear Boussinesqtype diffusion equation for the wave height, with coefficients dependent on the conductivities and depths of each layer. In the limit of equal conductivities of all layers, the equation reduces to the single-layer result recorded in the literature. The model equations are numerically integrated in the case of an incident monochromatic wave hitting the layers. The results exhibit dissipation and also a downstream net height rise at infinity. Wave transmission coefficient in three-layer porous media with conductivity of mangrove is discussed. Numerically, propagation of an initial solitary wave through a porous medium shows the emergence of wave reflection and transmission that both evolve as permanent waves. Additionally we examine the impact of a solitary gravity wave on a porous medium breakwater.


Introduction
As reported in Dahdouh-Guebas (2006) and Kathiresan and Rajendran (2006), mangrove forests along coastal areas provide shelter against storms and floods, as well as tsunami.Mangrove belts can absorb wave energy and reduce wave amplitude and velocity.The amount of wave amplitude reduction in mangroves depends on factors such as water depth, coastal profile, and also mangrove species and mangrove thickness.Due to the complexity of the problems, research on mangrove protection mostly relies on field measurement data from certain areas.For instance, Teh et al. (2009) conducted a numerical and experimental study of tsunami mitigation by mangroves (Avicennia officinalis) on the coast of Penang, Malaysia.Liu et al. (2003) studied flow resistance of mangrove trees using a depth-averaged hydrody-namic numerical model.The model and the vegetative parameter were tested using on-site mangrove data (Kandelia sp.) along Keelung River, Taiwan.Masda et al. (1997) reported on the experimental study of mangroves (Kandelia sp.) as a coastal protection in the Tong King delta, Vietnam.
Motivated by the protective function of mangrove in mitigating high waves, here we develop a study of gravity waves in a three-layer porous medium.The three layers here are leaves, trunks, and roots layers of mangrove forest.Wave interaction with non-uniform mangrove has been considered in the following literature.Vo-Luong and Massel (2008) study wave energy dissipation in a three-layer mangrove above arbitrary depth.They use velocity potential formulations with dissipation to develop a predictive model for random waves in a non-uniform forest in water of changing depth.Massel et al. (1999) study wave attenuation in mangrove forests consisting of roots and trunks.
Generally, studying wave propagation through a porous structure has many practical applications in coastal and ocean engineering.As a breakwater, porous structures reduce the amplitude of waves propagating into the structure.Studies in this area are quite extensive; here we mention only a few of them.Sollitt and Cross (1972), Madsen (1974), and Sulisz (1985) developed a linear wave model for progressive waves in a finite thickness of breakwater, and studied wave reflection and transmission far away from the porous structures.Experimental and numerical studies of wave interaction with a finite thickness of porous structure are discussed for instance in Fernando et.al. (2008), van Gent (1995), Lynett et.al. (2000), Lin and Karunarathna (2007), Liu et al. (1999), andScarlatos andSingh (1987).These authors use the modified shallow water fluid flow equations to take into account frictional dissipation.Other authors, Chwang andChan (1998), andLiu andWen (1997), use porous media flow equations as described by Darcy's Law to describe flow through the structures.Since for thick groves, viscous effects dominate inertial effects, for our study here we adopt Darcy's assumption, and we develop an asymptotic model for gravity waves in a rigid isotropic three-layer porous structure.
The organization of this paper is as follows.The governing equations are presented in Sect. 2. In Sect. 3 asymptotic methods are used to derive a simpler set of equations; the wave height is shown to satisfy a Boussinesq equation.It is a diffusion equation with higher-order nonlinear and dispersion terms.In Sect. 4 a predictor-corrector finite-difference scheme is used to simulate the movement of a monochromatic wave through the medium.Simulation shows amplitude reduction of waves in three-layer porous media.Numerical simulations clearly display the nonlinear effect: the rising of water levels at the right ends of rather long porous media.In Sect. 5 we use typical hydraulic conductivity of mangrove to calculate the wave transmission coefficient for a porous belt with certain length.In Sect.6 we show that the newly derived three-layer Boussinesq model can be combined with the Boussinesq equation for free water to simulate a solitary wave passing through a porous medium, in which the predictor-corrector scheme follows through.Conclusions are given in Sect.7.

Governing equations for flow in three-layer porous media
In this section, we formulate the governing equations for free-surface flow in a three-layer porous medium over an infinitely long horizontal axis.Let η( x, t) denote the free surface elevation measured from the undisturbed water level; see Fig. 1.The porous medium consists of three layers: an upper, middle, and lower layer with thicknesses h1 , h2 − h1 , and h3 − h2 , respectively, denoted as We assume that each layer ¯ i for i = 1, 2, 3 is a rigid isotropic porous medium containing an ideal fluid of depth h3 ( x).
Inside the porous medium we employ Darcy's law: with K i hydraulic conductivity for each layer i = 1, 2, 3, with dimension m/sec.In each of the three layers, ¯ i mass conservation requires Along the free surface, the pressure P is constant and can be taken as zero; we then obtain a boundary condition We assume that each layer Ωi , for i = 1,2,3 is a rigid isotropic porous medium containing an ideal Inside the porous medium we employ Darcy's law with K i is hydraulic conductivity, for each layer i = 1,2,3, with dimension m/sec.In each of the 75 three layers Ωi mass conservation requires Along the two interfaces z = − hi the following boundary conditions hold: Equation ( 7) comes from continuity of the pressure P , while Eq. ( 8) comes from the assumption that fluid influx and outflux across the interface are the same.The kinematic boundary condition along the surface is Along the impermeable bottom z = − h3 ( x), the normal flux vanishes: We summarize: the governing equations for ideal fluid flow in the three-layer porous media are Eqs. (5-10).

Asymptotic expansion method leading to Boussinesq equation
Assuming small amplitude long wavelength, we look for asymptotic solutions ¯ i for i = 1, 2, 3 of the governing Eqs.(6-10).Let L be typical length of the porous medium, h3 the fluid depth, and a typical free surface displacement.
Introduce the following non-dimensional variables, written without bars as follows: Note that conductivity K 1 is scaled in the non-dimensional time variable t, and also that h 3 ≡ h3 / h3 = 1.Further, we denote two small parameters Under the same assumption, Liu and Wen (1997) implement the asymptotic expansion method for the case of a single layer.In this section we extend this approach for the case of three-layer porous media.Compared to this approach, our approach is more direct.Rewriting the governing Eqs.(6-10) in non-dimensional variables, equations for the first layer are The equation for the second layer is The equations for the third layer are The matching conditions along the two interfaces z = −h i are Further, we restrict to Boussinesq approximation, in which , and look for a solution in the form of the series We implement the asymptotic expansion method in the normalized governing Eqs. ( 13) and ( 20), by solving the corresponding equations layer by layer, starting from the first O(1) terms, to higher order O(µ 2 ) terms, and further to O(µ 4 ) terms.Hence, an explicit approximate expression of i1 for i = 1, 2, 3 can then be obtained; see Eqs. (A5), (A7), and (A10) in Appendix A, respectively.Substituting the approximate potentials into the interface condition (20) along z = −h 2 yields the following equation of Boussinesq type: with coefficients A and B dependent on the conductivity and thickness of each layer, or explicitly with The detail derivation of Eq. ( 21) is given in Appendix A. The Boussinesq equation above, written in normalized variables, is a model for gravity waves in an emerged three-layer porous medium; it is a diffusive equation with higher-order nonlinearity and dispersion.
The horizontal flux of fluid in the three-layer porous media is given by Since the potential i for each layer was already approximated, substituting ( i ) x , for i = 1, 2, 3 yields the following relation: with Hence, Eqs. ( 21) and ( 25) are Boussinesq models for gravity waves in three-layer porous media.It is interesting to note that the flux and surface displacement equations are only weakly coupled; after solving for the displacement using Eq. ( 21), the flux can be obtained using Eq. ( 25).This is different from the fully coupled equations that arise from the corresponding shallow water (free) flow problem.The reason for this can be explained as follows.Under Darcy's assumption, we get boundary condition ( 14) which gives us the first-order term 10 (x, t) = η(x, t).From the asymptotic expansion framework, the first-order terms of Eq. ( 19) result in 20 (x, t) = 10 (x, t) = η(x, t), and further 30 (x, t) = 20 (x, t) = η(x, t).Since the first-order term of all potential i , i = 1, 2, 3 is just η(x, t), the leading-order term in the continuity equation η t + Q x = 0 is just a diffusion equation with a diffusive term proportional to η xx .
Next, we take a zero test where the three layers reduce to one layer.For this purpose we take h 1 = h 2 = h 3 for which K 3 = K 2 = K 1 , then the dispersion coefficient and diffusion coefficients become A = 1 3 h 2 3 and B = h 3 respectively, and further C = 1 3 h 2 3 .Hence, Eqs. ( 21) and ( 25) reduce to respectively.The set of Eqs. ( 27) and ( 28) is equivalent to the Boussinesq equations for flow in one-layer porous media as derived in Liu and Wen (1997) using a slightly different asymptotic expansion.21) and ( 25) for the three-layer model with the corresponding results (27, 28) for the single layer, it can be seen that waves propagating into three layers behave like those propagating into a single layer.Note that the model developed in this paper is based on the assumption that the waves are long compared to depth.For waves of wavelength comparable with the depth (or shorter), the velocity profile would vary significantly with depth, and the single-and three-layer models would likely produce different results.
For small amplitude displacements η, Eq. ( 27) can be approximated by the linear diffusion equation with wave-like solutions of the form exp i(kx − ωt).Substituting this yields iω The relation above means that waves of frequency ω propagating into the system decay over a distance on the order of a wavelength 1/k.The penetration length is strongly dependent on wave frequency ω, which is approximately √ 2B/ω.Numerical solutions to Eq. ( 21) are obtained using the predictor-corrector MacCormack method.The method has second-order accuracy O( x 2 , t 2 ), with x and t the partition length of spatial and time variables, respectively.Since there is a combine derivative in term η xxt , we first rewrite Eq. ( 21) as The first step in the MacCormack method is the predictor Fig. 2. A periodic tidal wave in a porous medium as the result of ( 21) with a significant wave r The above procedure is then used to simulate wave evolution in a three layer porous tially, we have a still water level η(x,0) = 0, then we impose a monochromatic wave i left η(0,t) = sin(t) or η n 1 = sin(t n ).Along the right boundary, we implement a trans ary, simply by η n+1 N x = η n N x−1 .For computations we take ε = µ 2 = 0.3, and ∆x = ∆ take a three layer porous model of the same thickness, with conductivity K 1 = 0.5, K 205 K 3 = 0.5.The surface profile reduced in the porous media is depicted in Figure 2.This matrix is diagonally dominant for small values of µ 2 , hence it is invertible.Correction for n+1 j , denoted as n+1 j , is calculated using the average of old values at time t n and predicted values at time t n+1 , as follows: Having n+1 j for j = 1, • • • , N x, we can directly compute η n+1 j for the new time step, using η = D −1 .The above procedure is then used to simulate wave evolution in a three-layer porous medium.Initially, we have a still water level η(x, 0) = 0, then we impose a monochromatic wave influx from the left: η(0, t) = sin(t) or η n 1 = sin(t n ).Along the right boundary, we implement a transparent boundary, simply by η n+1 N x = η n N x−1 .For computations we take ε = µ 2 = 0.3, and x = t = 0.02.We take a threelayer porous model of the same thickness, with conductivity K 1 = 0.5, K 2 = 0.47, and K 3 = 0.5.The surface profile reduced in the porous media is depicted in Fig. 2.
From the same computation we plot waves entering the porous medium in four different scaled frequencies: ωt = 0, π/2, π, 3π/2; see Fig. 3. Apart from wave reduction, we clearly observe the increase of mean water level in the far right end of the porous medium.This nonlinear phenomenon is derived in Liu and Wen (1997) as the analytical asymptotic solution of the single-layer Boussinesq model ( 27), in which η → εR 31 4h 3 for x → ∞.This well-known but surprising result was observed by Nielsen (1990) in the field.1027 Fig. 2. A periodic tidal wave in a porous medium as the result of ( 21) with a significant wave reduction.The above procedure is then used to simulate wave evolution in a three layer porous medium.Initially, we have a still water level η(x,0) = 0, then we impose a monochromatic wave influx from the left η(0,t) = sin(t) or η n 1 = sin(t n ).Along the right boundary, we implement a transparent boundary, simply by η n+1 N x = η n N x−1 .For computations we take ε = µ 2 = 0.3, and ∆x = ∆t = 0.02.We take a three layer porous model of the same thickness, with conductivity K 1 = 0.5, K 2 = 0.47, and 205 K 3 = 0.5.The surface profile reduced in the porous media is depicted in Figure 2. 8 Fig. 3. Monochromatic waves enter the porous medium with four different scaled frequencies: ωt = 0, ωt = π/2, ωt = π, ωt = 3π/2, computed using ε = µ 2 = 0.3.All four cases show the rising of mean water level as x → ∞.

Discussion of possible application to a mangrove belt
Consider a mangrove belt as a three-layer porous medium with upper, middle and lower layers corresponding to leaves, trunks and roots of mangrove.We propose an approximate model for wave evolution in mangrove as Eq. ( 21).As a wave breaker, our main concern is the amount of wave energy absorbed by the breakwater.Since wave energy is proportional to the square of wave amplitude, we concern ourselves with the amount of amplitude reduction in a porous belt.Wave transmission coefficient K T is a quantity defined as a ratio between transmitted wave and incident wave amplitudes.If the surface profile in a mangrove belt as a solution of Eq. ( 21) is known, then the K T curve is just the normalized displacement of wave entering the porous medium with zero phase.An experimental study in Susilo and Ridd (2005) reported hydraulic conductivity for mangrove (Rhizophora sp.) as about 10 m day −1 or approximately 0.01 cm s −1 .Here, we solve Eq. ( 21) for three-layer mangrove of equal thickness, with conductivity K 1 = 0.01 cm s −1 , K 2 = 0.03 cm s −1 , K 3 = 0.02 cm s −1 , and ε = µ 2 = 0.3.We impose a monochromatic wave with period 10 s, or frequency ω = 2π/10 rad s −1 .The resulting wave transmission coefficient is given in Fig. 4. Note that the curve is plotted with respect to the physical coordinate x = Lx = ( h3 / √ ε)x.In physical dimension, a point P in Fig. 4 means the following.For water depth h3 and the incident wave amplitude a that satisfies a/ h3 = ε = 0.3, after passing a porous medium with length x = 102.2cm, the wave transmission coefficient is 0.8374, which corresponds to an amplitude reduction of ≈ 16 %.
For a mangrove forest with certain hydraulic conductivity, the wave transmission curve K T clearly depends on the choice of ε(= µ 2 ), and on wave frequency ω.Further, we perform calculations using mangrove conductivity as above, using ε = 0.2 − 0.6 for several wave frequencies, and the results are given in Table 1.For mild wind waves with typical frequency 2π/6 rad s −1 , a porous belt length 100 cm gives For a mangrove forest with certain hydraulic conductivity, the wave transmission curve depends on the choice of ε(= µ 2 ), and on wave frequency ω.Further, we perform calcul 235 mangrove conductivity as above, using ε = 0.2 − 0.6 for several wave frequencies and t given in Table 1.For mild wind waves with typical frequency 2π/6 rad/sec, a porous bel cm gives wave transmission coefficient 0.74 to 0.82.We note that the recommendation holds for Boussinesq waves, i.e. small amplitude long waves, and the water depth is as constant.Moreover, reflection from the beach has been neglected.In coastal engineering context it is important to determine how much of the energy of a 'soliton' is absorbed by a breakwater, here the three layer porous media.With this examine the impact of a soliton propagating across the ocean (depth h 3 ) on a finite le medium of the above type.The soliton is assumed to impact normally on the medium.O 245 interest is the dependence of the transmission coefficient K T on the length L of the por with certain conductivities.The impacting wave is assumed to be long and of small a that the standard Boussinesq approximating equations Fig. 4. Wave transmission coefficient curve from numerical computation using = µ 2 = 0.3.

Reflection of a solitary wave
In a coastal engineering context it is important to determine how much of the energy of an incoming "soliton" is absorbed by a breakwater, here the three-layer porous media.

S. R. Pudjaprasetya: Three-layer porous breakwater
With this in mind, we examine the impact of a soliton propagating across the ocean (depth h 3 ) on a finite-length porous medium of the above type.The soliton is assumed to impact normally on the medium.Of particular interest is the dependence of the transmission coefficient K T on the length L of the porous medium with certain conductivities.The impacting wave is assumed to be long and of small amplitude such that the standard Boussinesq approximating equations are used in the free-water zone; see Whitham (1973).These equations are governing equations for bi-directional gravity waves written as variables surface elevation η and horizontal momentum q.In the porous domain, we implement the threelayer Boussinesq Eqs. ( 21) and ( 25) in physical variables; they read with physical parameters Ā ≡ h2 3 A, B ≡ h3 B, C ≡ h2 3 C.We simulate a solitary wave initially located upstream of the emerged porous medium.As time progresses, the wave propagates through the porous medium and further generates reflected and transmitted waves.Here we show that the numerical predictor-corrector method as explained in Sect. 4 follows through and can be used for solitary wave simulation.
The predictor-corrector method for Boussinesq equations in free areas (32, 33) is described below.We first rewrite Eq. (33) as The new variable ψ introduced here is just for computational needs.Its values at grid points relate to q values through a linear correspondence analogous to Eq. ( 30).The first step implements forward time forward space in Eqs. ( 32) and ( 36) for calculating predicted values, indicated by overlines, as follows: The correction step is then Table 2. Wave transmission coefficient for several length of porous medium, calculated using conductivity K1 = 0.7, K2 = 0.75, K3 = 0.7 cm/sec for a/ h3 = 0.1 − 0.4.
Here we show how the solitary wave evolves when it passes an emerged porous medium.The initial solitary wave profile of Boussinesq with amplitude a is as follows 3 290 which travels undisturbed in shape in free water region.Figure 5 presents simulation of an incoming solitary wave with amplitude a = 3 cm on a water depth h3 = 30 cm passing through a porous medium (indicated in grey) with length L = 60 cm.As the soliton passes through the porous media, reflected and transmitted waves occur.Amplitude ratio between these two waves depends on the precise detail of the emerged porous medium, i.e. its length L and its hydraulic conductivity Ki, i = 295 1,2,3.Smaller conductivity yields smaller amplitude of transmitted wave.For fixed conductivities, wave transmission coefficient are calculated for several length of porous medium, and the results are given in Table 2.
Length L of porous belt 20 cm 30 cm 40 cm 50 cm 60 cm 0.45-0.500.40-0.450.27-0.330.22-0.280.20-0.26 where we use forward space to compute q x | n j , η x | n j and backward space to compute q x | n+1 j , η x | n+1 j .Consider a computational domain with a porous domain in the middle, located on [0, L].The approximate Eqs. ( 32) and (33), and the approximate of three-layer model Eqs.( 21),(25) are combined.Vertical matching conditions come from the continuity of surface elevation where − and + denote the left-hand side and the right-hand side of the interface, respectively.In physical variables continuity of horizontal momentum is as follows: where the term K 1 Q denotes the reduced horizontal momentum in porous media.These matching conditions are commonly used for simulations of wave structure interaction, such as in Liu andWen (1997), Lynett et.al. (2000), Scarlatos and Singh (1987), and Wang and Li (2002).
Here we show how the solitary wave evolves when it passes an emerged porous medium.The initial solitary wave profile of Boussinesq with amplitude a is as follows: η( x, 0) =a sech 2 3a 4 h3

)
Along the left or right boundary, centre difference is replaced with one-sided difference.The overline indicates predicted values.When the predicted values n+1 j for j = 1, • • • , N x are known, the predicted values η n+1 j for j = 1, • • • , N x can be obtained by solving the linear relation Dη = , with

Fig. 2 .
Fig. 2.A periodic tidal wave in a porous medium as the result of Eq. (21) with significant wave reduction.

2406
Reflection of a solitary wave