EXTENDED BOUSSINESQ EQUATIONS FOR WAVES IN TWO POROUS LAYERS

In this study we continue the work of Vu et al. (2018) [Coastal Eng. 139, 85-97] to develop an extended Boussinesq model that predicts the propagation of water waves in two porous layers. The first and second layer can be a water layer or a porous layer. The inertial and drag resistances are considered in the developed model. After being successfully validated against the analytical solutions, the model is used to simulate waves propagating over a submerged triangular porous bar. The numerical results show good agreement with the physical experimental data of Hsiao et al. (1998) [Proc. Royal Society of London A 458, 1291-1322].


INTRODUCTION
There exist multiple porous layers in coastal area. Submerged breakwaters are composed of two layers, i.e., upper non-porous water layer and lower porous layer. Rubble mound breakwaters are composed of two layers with different porosities. Sand beach is composed of one porous layer on the land side from the coastline. Also, it is composed of two layers with upper water and lower porous layers on the sea side. We need to have a wave model for such multiple porous layers. Until now, the Boussinesq equations have been developed for two layers with upper non-porous water layer and lower porous layer (Cruz et al., 1997;Hsiao et al., 2002). Liu and Wen (1997) developed conventional Boussinesq equations in shallow water for waves inside porous media including drag resistance but neglecting inertial resistance. These models cannot be applied for rubble mound breakwaters or sand beach. Lee et al. (2014) and Vu et al. (2018) developed extended Boussinesq equations for waves in porous media which considers drag and inertial resistances. The numerical solutions of their model are well verified by comparison with exact solutions and physical experiment data. In this study, we continue the work of Lee et al. and Vu et al. to develop a model for waves in two porous layers which can solve the aforementioned problems. The validity of the developed model is examined by comparing numerical results with analytical solution and experimental data.

DEVELOPMENT OF GOVERNING EQUATIONS
Setting-up boundary value problem In order to develop the extended Boussinesq equations for waves inside porous layers we first set up boundary value problems for waves inside multiple porous layers. The whole domain is vertically divided into the 1st, 2nd,…, and j th  layers which are numbered from the top to the bottom layers with different porosities, as shown in Fig. 1. The free surface is located at the 1st layer. Since the porosity is uniform, the continuity equation inside the j th  permeable layer is given by where j p is the pore pressure, j D is the permeable drag resistance term, and j I is the inertial resistance term. It should be noted that the continuity and momentum equations (1) and (2), respectively, are expressed in terms of the seepage velocity of the pore water. Several people defined the drag resistance term differently. Ergun (1952) define the drag resistance term in the Forchheimer (1901) type using a volume-averaged discharge velocity In our study, we use Ergun's definition of D in terms of the seepage velocity instead as where l  and t  are coefficients which represent the laminar and turbulent flow resistances, respectively,  is the kinematic viscosity of water, and d is the solid size. The mathematical form of Eq.
(3) including the parameters  ,  , and d can be derived from the Navier-Stokes equation (Irmay, 1958) or from the Reynolds equation (Burcharth and Andersen, 1995). The inertial resistance term j I is given by where  is porosity, and  is added mass coefficient. In unsteady flow, the inertial resistance term is necessary to consider the divergence and convergence of streamlines in the presence of the solid material. The resistance term should consider both the local and convective accelerations. In Eq. (4), the value unity which is added to the added mass coefficient is to consider the inertial resistance of the water with the volume of solid material, and the added mass coefficient is to consider the inertial resistance in view of geometrical smoothness of the solid material. When the layer is filled with water, then 1   and the inertial resistance term becomes zero. When the layer is fully filled with the solid material, then 0   and the inertial resistance term becomes infinitely large. Several people proposed different momentum equations including the drag and inertial resistance terms. Sollitt and Cross' (1972) momentum equation is the same as the present momentum equation (2) except that the convective acceleration is neglected to get a linear solution. van Gent's (1994) momentum equation is the same as the present momentum equation (2). However, he defined the inertial resistance term only in the local acceleration as where the value unity, which is omitted, should be added to the added mass coefficient. Cruz et al. (1997) proposed the momentum equation as where the drag resistance term is defined in terms of the seepage velocity as where ' l  and ' t  are coefficients which represent the laminar and turbulent flow resistances, respectively, and the inertial resistance term is defined as At the interface between the layers of different porosities, the pore pressure should be continuous in terms of the local pore pressure instead of the volume-averaged pore pressure. The difference of pressure forces between the layers will act on the solid material in the layer of smaller porosity. To the contrary, the normal flux at the interface should be continuous in terms of the volume-averaged discharge velocity. Thus, at the interface between the j th  and   1 j th   layers, both the pore pressures at and normal fluxes through the interface are continuous as 1 , We get the Boussinesq equations by specifying the boundary value problem with a governing equation and boundary conditions. The seepage velocity potential is defined as Then, the continuity equation (1) becomes At the free surface, the dynamic and kinematic boundary conditions (12) and (13) become, respectively, At the interface, the continuities of the pore pressure and normal flux (15) and (16) become, respectively, The variables are normalized using the relevant characteristic length and time as where l is the wavelength, 0 h is the maximum water depth, a is the maximum amplitude of the water surface elevation, t is the wave period. When the normalization is applied, the terms in the governing equation and boundary conditions will group according to two non-dimensional quantities where  is the nonlinearity parameter and  is the dispersivity parameter. Omitting the primes for convenience, the continuity equation (18) and boundary conditions (19)-(23) become, respectively, The velocity potential can be expressed as a power series in the vertical coordinate given by Substitution of Eqs. (33) and (35) into the bottom boundary condition (29) gives And, the velocity potential functions in the lowest layer can be expressed in terms of ,0 J  using Eqs. (37) and (38) Thus, the velocity potential in the bottom layer J  can be expressed to the order of  

Development of Boussinesq equations for waves in two porous layers
For two permeable layers, all the subscripts j in the variables are 1, 2 j  . Thus, from Eq. (41), the velocity potential in the bottom layer 2  becomes The velocity potentials 1  and 2  are related through the interface boundary conditions (30) or (31).
We obtain the velocity potential function in the 1st layer 1,1  by applying the interface boundary condition (31) And, the velocity potential functions in the 1st layer 1,2 1,3 , ,...

 
can be expressed in terms of 1,0  and 2,0  using Eqs. (37) and (43) The momentum equation in the 1 st layer is obtained by substituting Eq. (46) into the dynamic freesurface boundary condition (27) where 2,0 is the seepage velocity in the 2 nd layer. The momentum equation in the 2 nd layer is obtained by substituting Eqs. (42) and (46) into the interface continuity of pore pressure (30) Or, using Eq. (47) Here, we use the depth-averaged velocity defined as After substituting Eq. (42) into Eq. (50) with 2 j  , the seepage velocity at the lower boundary of the 2 nd layer can be expressed in terms of the depth-averaged velocity as Substitution of Eqs. (51) and (52) Substitution of Eqs. (51) and (52) The continuity equation is obtained by substituting Eq. (46) into the kinematic free-surface boundary condition (28) as Eqs. (53)-(55) are the set of Boussinesq equations for waves in two permeable layers.
If the 1st layer is not in porous media but in clean water, then 1 1   , 1 1   and 1 0   , and thus Eqs.
(53)-(55) become the Boussinesq equations of Cruz et al. (1997) which are for waves on a single permeable layer. If the 1st layer is in clean water and the interface between the 1st and 2nd layers is rough, then 1 1   , 1 1   and 1 1 where f is the friction factor related to the shear stress on the rough interface. Further, if waves are very long, the resulting equations (53)-(55) are similar to the shallow-water equations of Kobayashi (1986).

Extended Boussinesq equations for waves inside two permeable layers
The momentum equation (56) From the lowest-order momentum equation in the 2 nd layer, we can get the following relation Using this relation, the last term in the momentum equation (59) The corrector step is iterated until total relative error between successive results is within a certain limit.

Linear waves inside two porous layers
The computational domains for waves propagate in 2 porous layers are given in Figs. 1 (a) and ( ). In this case, the wave amplitude is negligibly small and thus, we can use the analytical solution for linear waves given by k is imaginary part of the complex wavenumber. The numerical solutions of surface elevation and wave amplitude are almost the same as the analytical solution for all the cases. In deep water with the first non-porous water layer and the second porous layer, waves propagate with almost no energy dissipation. Waves in shallow water are damped more than in deep water. , d=19 mm .The computational domain with 9 wave gauges is given in Fig. 3.  Fig. 4 shows that the numerical solutions of water surface elevations are close to the experimental data at all gauges. It should be noted that using the present model we do not need to use any matching conditions at the interfaces of the submerged porous breakwater. It should be noted that gauge G1 was installed in front of the toe of the submerged breakwater while gauge G7 was at the crest of the breakwater. This paper derived an extended Boussinesq model for waves in two porous layers considering inertial as well as drag resistances. For linear waves in two porous layers, numerical solutions are close to analytical solutions. For waves above a submerged triangular porous breakwater, numerical solutions are well compared with the experimental data of Hsiao et al. (2002).