Abstract

A numerical model is presented for simulation of hydrodynamics of a 2D vertical free surface domain consisting of an arbitrary partitioned porous and non-porous area. To this end, modified Navier–Stokes equations are considered which could be applied in surface water and in subsurface flows, simultaneously. A wide range of Reynolds number has been considered, from which non-Darcy effects have also been taken into account. A fractional step method has been used in the time discretization procedure, where the convection and diffusion terms are separated from the pressure term while solving the momentum equations. To include the variation of surface elevation in computation, the domain has been divided into two parts, namely, ‘interior subdomain’, which never gets dry during the simulation period, covered by fixed unstructured triangular grids and ‘top layer’ with only a one layer structured grid, the position of which varies with the water surface. The validation of the proposed model has been achieved by comparison of its results with both theoretical and experimental data reported in the literature.

INTRODUCTION

Modeling of flow in porous media is a classic problem in civil engineering, and it plays a critical role in analysis and design of civil engineering structures (Wang & Ren 1993; Liu et al. 1995). Because of its mathematical challenges, it is also interesting from a mathematical point of view. Thus, a considerable amount of literature has been published on flow in porous media considering the notion of the phenomenon, numerical modeling and analytical solution (Biot 1956, 1962, 1995; Zienkiewicz et al. 1980; Sahebkar & Eskandari-Ghadi 2016, 2017). In traditional studies, Darcy's equation (Darcy's law) is used and which represents a linear relation between the gradient of pressure and the velocity vector (Biot 1941; Bejan 1984; Zienkiewicz & Shiomi 1984; Griffiths & Fenton 1997). This assumption is valid as long as the Reynolds number is very low, while it is not a reasonable assumption in many other engineering boundary value problems, where the grain sizes of porous media are not fine and the flow velocity is considerable. A number of studies have been made to modify Darcy's law (Forchheimer 1901; Balho et al. 2007). The idea of including velocity by a power of two (or more) in the Darcy equation was introduced by Joseph et al. (1982) and known as Darcy–Forchheimer law, converting the linear Darcy's equation to a nonlinear one. Brinkman (1947) suggested another method to modify the Darcy equation in porous media with large porosity. The logic of the mentioned modification was the force due to the viscosity of fluid (water) to be taken into account in large porosities. This force has been considered by the Brinkman additional term. Several attempts have been made to use some combination of Forchheimer and Brinkman terms to form the momentum equation (Wooding 1957; Vafai & Tien 1981; Nield 1991; Chen & Chen 1992; Chen et al. 1998; Wu et al. 2008). Recently, investigators have used Navier–Stokes (N-S) type equations in modeling of flow in porous media (Akbari & Namin 2013; Li et al. 2016). One of the challenges of modeling of flow can be seen in the problems where the domain of interest consists of both porous and non-porous media (Konga et al. 2010; Yuan et al. 2011; Akbari & Namin 2013; Li et al. 2016). These two parts are traditionally solved in two different procedures, and the interaction of these two domains is imposed by using the result of the model in the previous time step on each part for setting the boundary condition of the other part, and vice versa. Recently, a few researchers have tried to solve the coupled equations for surface water and subsurface flow to obtain more accurate results (Konga et al. 2010; Larese et al. 2014; Li et al. 2016). For example, Akbari & Namin (2013) used this approach in a smoothed particle hydrodynamics (SPH) numerical method in wave interaction with structures. Their developed model solves porous and pure fluid (water) flows simultaneously by means of an equation that is equivalent to the unsteady 2D N-S equations for the flows outside the porous media, while the extended Forchheimer equation has been used for the flows in the porous media. In their model, the interface boundary between pure fluid (water) and porous media is effectively taken into account by the SPH integration technique (Akbari & Namin 2013).

One of the main challenges in flow modeling, in the context of this paper, is the tracking of free surface variations. The hydrostatic assumption has been used by many authors in previous researches (Heniche et al. 2000; Casulli 2017). However, the solutions based on the simple hydrostatic assumption are not very accurate in most engineering problems such as short wave test cases, where the vertical velocities and accelerations are not ignorable. For this reason, detailed investigation based on the hydrodynamic treatment of the phenomenon is required where the vertical plane of the model area has to be divided into a finite number of elements. Two different approaches have been followed based on hydrodynamic assumption, which are updated mesh and fixed mesh approaches. Studies with updated mesh mostly use an arbitrary Lagrangian Eulerian (ALE) technique (Braess & Wriggers 2000; You & Bathe 2015). When using the ALE technique in simulations, the mesh of the domain has to be regenerated and the nodes on the boundaries of the domain have to move along with the materials to precisely track the boundaries accordingly. In the fixed mesh (Eulerian) approach, on the other hand, the level set method usually is used by which the free surface is represented by a signed distance function that takes zero on the free surface, negative values in the fluid domain, and positive values outside the fluid domain (Lin et al. 2005; Larese et al. 2014; Cruchaga et al. 2016). In most of the applications, the above two mentioned methods could track surface evolution accurately; however, the main disadvantages of these methods are their high cost and complexities (Braess & Wriggers 2000; Larese et al. 2014; You & Bathe 2015; Cruchaga et al. 2016). The improved accuracies due to the mentioned methods do not justify suffering the high computation expenses and model complexities in some of the real applications.

In this paper, a new approach for tracking the surface of combined surface water and subsurface flows is presented, where a semi-implicit numerical method with unstructured grid is used. In this new approach, the flow in clear fluid (non-porous) and porous media are described by the N-S type of equations. Moreover, resistance forces and porosity in the porous media have been taken into account, without being limited to low Reynolds number. For free surface tracking, the domain has been divided into the interior and the top layer subdomains. The 2D interior subdomain is discretized with unstructured triangular elements, which are fixed at all time steps, and fully N-S type of equations have been employed to obtain the discrete equations in this subdomain. On the other hand, 1D linear elements are used in the top layer, which contains the free surface variations. Shallow water type equations have been applied to define the flow conditions in this layer, which means that the pressure distribution is obtained based on hydrostatic assumptions in the top layer subdomain. Therefore, a combination of hydrostatic–hydrodynamic methods has been applied to track the free surface variations. This approach reduces the computational cost and could be applied in a wide range of engineering applications. On the other hand, fractional step method (Kim & Moin 1985; Perot 1993; Despotis & Tsangaris 1995) has been used for temporal discretization while the finite volume–finite element method has been applied for space discretization of the equations.

It is worth mentioning that in the presented model, the main reason for using fully N-S equations among the shallow water equations, is the required accuracy of pressure distribution in the vertical direction. In 2D-plane models, vertical velocities and accelerations are not taken into account and hydrostatic assumption is instead assumed for pressure distribution in the vertical direction. However, in the presented model, vertical velocities and accelerations are considered in the calculations and thus the pressure distribution is calculated with more precision. One of the desired applications for the presented model is numerical simulation of flow at a vertical plane passing through the inlet and outlet. Also, this model is suitable for any short wave test case. A similar combination of hydrostatic and hydrodynamic method has also been used by Chegini & Namin (2012) in free surface water using finite volume approach. It should be mentioned that their model uses a z-level Cartesian grid system, but in this model, an unstructured triangular gridding is used. Moreover, one of the applications of this model is in modeling of surface water–groundwater interaction which was not a concern in their study.

In summary, in the present model, unlike hydrostatic assumption, vertical velocities and accelerations are taken into account. Moreover, surface water–groundwater interaction is considered automatically using unified equations, which also considers non-Darcy effects. A combination of finite element–finite volume is used to solve the governing equations.

In this paper, some experimental and theoretical test cases are used to evaluate the accuracy and performance of the proposed method. The tests have been selected so that all numerical properties of the proposed method can be assessed. Test cases mainly focus on surface–subsurface interaction modeling, non-Darcy effects, and non-hydrostatic flow conditions.

GOVERNING EQUATIONS

The domain of interest in this paper is a 2D vertical plane consisting of surface water and subsurface flow (porous media). A Cartesian coordinate system, whose z-axis is positive in the opposite direction of gravity acceleration and its x-axis is in the horizontal direction, has been applied (Figure 1). The domain includes both clear fluid and saturated porous media. Time-dependent flow condition imposes the free surface variations to the problem, resulting in a time varying interface in the upper side of the domain as moving boundary condition. Other boundaries of the model domain have been considered as some fixed defined geometries in which the boundary conditions could be either specified pressure (p) or flux (q) type. In this paper, a modified version of N-S equations have been used for both clear fluid and porous media, by which the flow in porous media has been described more accurately when compared with the classic Darcy law since parameters such as viscosity and inertial force have additionally been taken into account. It should be mentioned that Darcy's equation is extracted for viscous fluid at low velocities and thus it is applicable for low Reynolds numbers (Bejan 1984) and, as previously explained, the current research is free of such a restriction. Since the governing equations for the porous media and the clear fluid have been solved simultaneously, the effects of these two different domains on each other are considered automatically in a precise manner. Hydrodynamic assumption has been used for pressure change in depth in both clear water and porous media (Figure 1) except for specific height of the domain in the upper side of the model, called ‘top layer’. The top layer includes the height in which free surface evolutions may take place (Figure 1). Thus, as can be seen in Figure 1, the computational domain including porous and non-porous area, has been divided into two different zones, namely, ‘top layer’ and ‘interior subdomain’, which are separated by a horizontal divider line. A specific set of equations has been employed for the two mentioned subdomains. For the top layer, shallow water type equations have been used, while N-S type equations have been used for the interior subdomain. These equations are presented in the following sections.

Figure 1

A schematic presentation of the domain of interest consisting of an interior part overlaid by a top layer.

Figure 1

A schematic presentation of the domain of interest consisting of an interior part overlaid by a top layer.

Interior subdomain

The continuity equation for incompressible fluid flow over a unit width inside the porous (or non-porous) media may be expressed as follows  
formula
(1)
where is the Darcian or averaged velocity which is equal to the intrinsic fluid particle velocity () multiplied by the porosity () of the porous media (). The clear water will take unit as its porosity. The momentum equation, using extended Darcy's law, without being limited to porous media, can be derived as the following equation (Akbari & Namin 2013; Akbari 2014):  
formula
(2)
where p and are the fluid pressure and the velocity vector, respectively. The left-hand side of this equation is the inertia forces per unit mass caused by displacement acceleration and moment acceleration. On the other hand, the first term on the right-hand side of the momentum equation is the force caused by the pressure gradient, and the second term represents the force due to fluid viscosity. This term plays a major role in a domain with high porosities, providing the basis to expand porous media flow equations to N-S type equations proposed by Brinkman (1947). The third and the fourth terms represent the drag forces exerted on the fluid by the solid skeleton and eventually the last term () is body forces such as gravitational acceleration. in Equation (2), is the density of fluid (water), while is the effective viscosity by which the real viscosity () is modified since the intrinsic velocity has been replaced by the averaged velocity in the formulation (Li et al. 2016). In addition, and are empirical coefficients of linear and nonlinear forces in porous media, which have been discussed in the literature (Irmay 1958; Amao 2007; Teng & Zhao 2000). For example, Teng & Zhao (2000) proposed coefficients from experimental results for different materials:  
formula
(3)
where and k are dynamic viscosity and intrinsic permeability, respectively. is a parameter defined in terms of intrinsic permeability and porosity as (Teng & Zhao 2000):  
formula
(4)
where is a parameter depending on porosity that has to be obtained experimentally (Teng & Zhao 2000). The intrinsic permeability may be determined by the following equation if no experimental measurement is available (Urumovic & Urumovic 2014):  
formula
(5)
is the average diameter of grains forming the porous media. It should be mentioned that when the porosity value approaches unity, a clear water condition is presented and by neglecting the drag forces in Equation (2), this equation together with Equation (1) converts to the well-known N-S equations.

Top layer

The top layer, as shown in Figure 1, is a subdomain located over the horizontal divider line and bounded by the free surface, the physics of which have been expressed by shallow water type equations. The equations are derived from depth integrating of the modified N-S equations presented in the previous section. In the case where the horizontal length size is much larger than the vertical one, it can be indicated from the momentum equation that the vertical component of pressure gradient is almost hydrostatic. Thus, the momentum equation in z-direction is reduced to a gravity equation, implying that the horizontal velocity field is constant throughout the depth of the fluid. In this way, the fluid equations are obtained as:  
formula
(6)
 
formula
(7)
where is the free surface height varying in time and position (x), and is the component of the velocity in x direction and is the velocity component crossing the divider line and varying in time and x direction. The other parameters have been defined in the previous section. As for the interior subdomain, similar equations have been considered for both the porous and non-porous area.

NUMERICAL METHODS

As has been mentioned, the domain of interest has been divided into two subdomains (Figure 1). The first one is the top layer containing free surface variations and the second one is the interior subdomain. Since the divider line is chosen in such a way that the change of free surface always takes place above the divider line, the mesh of the interior subdomain is not affected by free surface fluctuations. It is worth mentioning that the numerical formulations are derived in an Eulerian framework in which individual fluid particles are not identified. Instead, a fixed mesh is defined, as shown in Figure 2. Pressure, velocity and all other flow properties are described as fields within the fixed mesh. In other words, each property is expressed as a function of space and time. For numerical solution, an unstructured grid is employed for the interior subdomain and a staggered layout has been selected so that in the interior subdomain the velocities are positioned in the centroid of each triangular element and the pressures are defined on the vertices of triangular elements, as shown in Figure 2. 1D-type linear elements have been used for the top layer where the pressure variables are located at two end nodes of the bottom faces of the elements and velocities are considered to be at the middle of the bottom face (divider line) as shown in Figure 2. It must be recalled that the horizontal velocity field is constant throughout the depth of the fluid in the top layer subdomain. It means the free surface tracking is conducted by a combination of hydrostatic–hydrodynamic methods. Numerical discretization of the equations is presented for two subdomains in the following two separate subsections. Fractional step method (Kim & Moin 1985; Perot 1993; Despotis & Tsangaris 1995) has been applied for temporal discretization of momentum equations and finite volume method has been employed to solve advection and diffusion terms by which the intermediate velocities, to be defined later, are obtained. Then, by combining the continuity and the momentum equations with omitted advection and diffusion terms, the Poisson equations are formed to be solved via finite element method to obtain pressure values. Since the equations for both clear fluid and porous media are solved simultaneously, their interaction is taken into account automatically in both domains.

Figure 2

(a) Schematic view of domain filled by elements; (b) interior triangular elements; (c) top layer; (d) variables locations in triangular elements; (e) variables locations in top layer elements.

Figure 2

(a) Schematic view of domain filled by elements; (b) interior triangular elements; (c) top layer; (d) variables locations in triangular elements; (e) variables locations in top layer elements.

Regarding the application of fractional step method, it has to be explained that different terms in N-S equations have different physical natures. By virtue of fractional step method, one can use an appropriate numerical method for each term. For example, we use explicit method for advection terms, but implicit methods for diffusion and Poisson equations. This helps get stability in numerical modeling. This is done with no more complexity and computational cost. It should be mentioned that it is possible to use different time steps for different terms during use of the fractional step method; however, in the presented model, a unique time step is used for all terms.

Interior elements

In the following, the vector form of momentum equation has been applied for brevity. Momentum Equation (2) is a time-dependent one, therefore the variables have to be computed in different time steps with a time interval of in a progressive manner. To do so, fractional step method is employed for the temporal discretization of these equations. In the numerical procedure, the intermediate velocities denoted as () are first calculated by splitting the momentum equation as:  
formula
(8)
 
formula
(9)
Equations (8) and (9), respectively, take the advection and diffusion terms into account, and are solved by finite volume method with cell-centered triangular elements. The Fromm (Fromm 1968; Venkatakrishnan & Barth 1989) second-order upwind explicit scheme and an implicit cell-centered method (Eymard et al. 2005; Breil & Maire 2006) are used to solve the advection and diffusion stages, respectively. Stability for the advection equation can be achieved by selecting a time step that satisfies the Courant condition. Besides, implicit methods are unconditionally stable and do not have time step limitations (Casulli & Cattani 1994; Vitousek & Fringer 2013). The obtained intermediate velocities are then used in the momentum equation (Equation (2)) to determine the velocity at time , denoted by as:  
formula
(10)
where is the velocity vector at time , and is the magnitude of the velocity at the same time. One may define in terms of yet to be calculated pressure values as:  
formula
(11)
The next step is to discretize the continuity equation (Equation (1)). Applying the standard Galerkin procedure (Zeinkiewicz & Taylor 2000) and the Gauss theorem, and with the use of shape functions as the weight functions, one may write:  
formula
(12)
where is the shape function (interpolation function) for triangular elements of the interior subdomain and is the normal vector to surface . One may write the special variation of in the form of , where the overline denotes the appropriate nodal values of pressure and considering constant shape function for velocities , they are written as , where is the value of the velocity located at the centroid of each triangular element. For brevity, , and are, respectively, summarized to and . Inserting the velocity at time from Equation (11) into Equation (12) results in the following equation for the time  
formula
(13)
Rearranging this equation results in:  
formula
(14)
which may be written in the matrix form as:  
formula
(15)
where and are the coefficient and force matrices, respectively, as:  
formula
(16)
 
formula
(17)

It is worth mentioning that in the presented model, shape functions for three-node triangular element for pressure and constant shape function for velocities are used in spatial discretization of Poisson equation. However, any n-node plane element can be utilized in the finite element procedure. As seen in Equation (11), the velocity is proportional to pressure gradient; because of this, the shape functions of pressure are chosen one order higher than the velocity shape functions.

Additional details about the standard Galerkin procedure can be found in Zeinkiewicz & Taylor (2000).

Top layer elements

With the use of fractional step method, the discretization of momentum equation in x direction (Equation (7)) for linear elements is similar to the interior part. As the first step of the fractional method, the following two equations have been solved to calculate the intermediate velocities at the top layer using finite volume method where linear elements based on cell-centered assumption is used:  
formula
(18)
 
formula
(19)
As it has already been explained, these intermediate velocities are used in the momentum equation with omitted advection and diffusion terms to get the velocities at time denoted as in terms of pressures as:  
formula
(20)
The next step is to discretize Equation (6). As for the interior subdomain, the finite element method has been applied here. Applying Gauss integration rule and using as weight function, which is the shape function for two-node linear elements of the top layer subdomain, the following equations result in:  
formula
(21)
 
formula
(22)
where is the length of the linear elements and are its two end coordinates. Hydrostatic assumption for this part results in the following formula for constant density:  
formula
(23)
Variations of the pressure values in the element can be defined in terms of the nodal values as , where the overline denotes the appropriate nodal values of pressure and the velocities, by considering constant shape function, are written as and , where and are, respectively, the velocities at the middle of the linear elements and at the divider line. By inserting these values in Equation (23) and applying Equation (24), the following relation may be derived:  
formula
(24)
and have been summarized to and , respectively. Replacing the velocity from Equation (21) in this equation, one gets:  
formula
(25)
from which the following equation may be inferred for the flux that crosses the divider line passing volume from top to bottom layer, or vice versa:  
formula
(26)

In fact, Equation (27) has taken into account the effect of the moving surface as a boundary condition to the interior subdomain that affects coefficient and force matrices. This flux equation is actually a combination of Dirichlet and Neumann boundary condition at the divider line for the interior subdomain, since neither pressure nor its derivation are known, but Equation (27) defines a relation between them.

The time derivative of the pressure that exists in Equation (27) has been discretized in a finite difference manner as:  
formula
(27)
The flux has then been written in the matrix form as:  
formula
(28)
where  
formula
(29)
 
formula
(30)
 
formula
(31)
 
formula
(32)
 
formula
(33)

Once the pressure values are computed, the velocities on interior subdomain and top layer can be obtained from Equation (11) and Equation (21), respectively, and finally the free surface elevation is determined with the use of Equation (24).

It should be mentioned that a stability check for time steps is done during simulation. In the developed model, the advection equations (Equations (8) and (18)) are solved explicitly, whereas implicit methods are employed to solve the diffusion (Equations (9) and (19)) and Poisson equation (Equation (15)). Since implicit methods are unconditionally stable, we need to check just the stability of advection terms. Therefore, restricting the time step to less than the length interval divided by the velocity would guarantee the stability of the model (Casulli & Cattani 1994; Vitousek & Fringer 2013).

As well, sparse data structure is used in numerical solution of diffusion and Poisson equations to reduce the cost of computation and memory storage consumption. To store the sparse matrices row, indexed sparse storage scheme is used and also conjugate gradient method is used to solve sparse linear systems. Details of the conjugate gradient algorithm and the storage scheme can be found in the book entitled Numerical Recipes in Fortran by Press et al. (1992).

Free surface tracking algorithm

Based on the formulations explained in the previous sections, the following steps have to be taken for the algorithmic procedure:

  1. Solving Equations (8) and (9) in the interior subdomain and obtaining the intermediate velocities

  2. Solving Equations (18) and (19) for the top layer and obtaining the related intermediate velocities

  3. Creating stiffness and force matrices for each element of the interior subdomain from Equations (16) and (17), respectively

  4. Repeating step 3 for the top layer using Equations (30) and (34)

  5. Assembling all matrices from steps (3) and (4) to one total coefficient matrix and one total force vector

  6. Solving the derived matrix equation for obtaining the nodal pressure values

  7. Calculating velocities in the interior subdomain with the use of Equation (11)

  8. Calculating velocities of the top layer with the use of Equation (21)

  9. Calculating the free surface level with the use of pressure values on the divider line.

MODEL EVALUATIONS

Traveling of solitary waves

A solitary wave entering from the left boundary of the domain has been tested to prove the validity of the formulations and procedure proposed in this paper. The wave amplitude used in this test has been considered as 0.5 m, while the water depth has been assumed as 5 m and the initial water elevation is set to zero (Figure 3). The analytical solution for this test case is (Namin et al. 2001):  
formula
(34)
Figure 3

Dimensions and boundary conditions of traveling solitary wave test (not to scale).

Figure 3

Dimensions and boundary conditions of traveling solitary wave test (not to scale).

Velocities from the left boundary have been imposed to the system with the following values:  
formula
(35)
where a, , , and are the wave amplitude, the mean water depth, the wave celerity, and the phase, respectively. The following formulations as analytical solution are presented for the wave celerity and the phase for this case:  
formula
(36)

Porosity in this example has been set to one. Mesh length has been considered as 0.5 m, which is the triangle side length in the interior subdomain. A schematic view of mesh plot is shown in Figure 4. Free surface (water) elevations have been obtained from the theoretical solution (Namin et al. 2001) and the present numerical model solution, which are shown in Figure 5. The comparison shows that the calculated wave surfaces are in a good agreement with the analytical dynamic solution, which proves both the validity of the formulations and the accuracy of the numerical solutions presented in this paper.

Figure 4

Schematic view of the solitary wave domain mesh plot (not to scale).

Figure 4

Schematic view of the solitary wave domain mesh plot (not to scale).

Figure 5

Comparison of numerical results with theoretical solution.

Figure 5

Comparison of numerical results with theoretical solution.

Furthermore, some small oscillations are observed in Figure 5 for the present solution, which are not related to numerical stability. In fact, most numerical methods with the order of accuracy higher than first, experience some oscillations. Whereas the order of accuracy in this model is from the second order, and therefore we see these oscillations. There have been several studies discussing these types of oscillation and attempting to remove them using different schemes like TVD, ENO and WENO (Osher & Chakravarthy 1986; Harten et al. 1987; Balsara & Shu 2000). However, capturing these oscillations are not a concern of the presented model. The values of oscillations are not significant in this example. The root mean square error for the present solution in comparison with theoretical results is 0.02 m.

Water table fluctuation in a laboratory model aquifer

Water table oscillations in porous media applied by a piston was experimentally studied by Parlange et al. (1984). Figure 6 shows a schematic view of the experiment. The average water depth in the reservoir is 27.6 cm. A sinusoidal water table elevation, with an amplitude of 9 cm and a period of 35 s has been outlined at the left boundary while the right and the bottom boundaries have been defined as no flow boundary conditions. Porosity and intrinsic permeability for this test have been set to 0.3 and , respectively. It should be noted that Reynolds number range is low in this example, so no nonlinear drag terms were required. This example has been solved in a number of previous researches (Liu & Wen 1997; Yuan & Lin 2009; Li et al. 2016). Further details of the example can be found in Parlange et al. (1984).

Figure 6

Schematic view of experimental design for a sinusoidal fluctuation designed by Parlange et al. (1984) (not to scale).

Figure 6

Schematic view of experimental design for a sinusoidal fluctuation designed by Parlange et al. (1984) (not to scale).

In the numerical model, the mesh length has been considered as 1 cm in a triangular shape (triangle side length is 1 cm). A schematic view of the mesh plot is shown in Figure 7. The time step in the present solution has been taken as 0.005 s, which is the same as the time step considered by Li et al. (2016). The results of the present simulation have been compared with the experimental data (Parlange et al. 1984) and theoretical solution (Liu & Wen 1997). Moreover, the results have been compared with the numerical solution of Li et al. (2016) at four points A, B, C, and D at distances of 0, 5, 18, and 33.5 cm from the left boundary. The comparisons are illustrated in Figure 8(a)–8(d). As can be seen, the agreement between results of the present solution, analytical solution, and those of Li et al. are very good. In addition, agreement between experimental results and present solutions are satisfactory. The difference between experimental and theoretical solutions could be relevant to measurements. Although the experiment was well conducted, the water level and frequency may not be very accurate, which was also commented on by Liu & Wen (1997), Yuan & Lin (2009), and Li et al. (2016).

Figure 7

Mesh plot for test domain (not to scale).

Figure 7

Mesh plot for test domain (not to scale).

Figure 8

Comparisons of experimental, theoretical, and Li et al. (2016) results with the current research results at four distances from the left boundary: (a) x = 0, (b) x = 5, (c) x = 18, (d) x = 33.5.

Figure 8

Comparisons of experimental, theoretical, and Li et al. (2016) results with the current research results at four distances from the left boundary: (a) x = 0, (b) x = 5, (c) x = 18, (d) x = 33.5.

Root mean square errors at four points A, B, C, and D are presented in Table 1 for the current model and for Li et al. (2016) results. It should be mentioned that Li et al. (2016) used a depth integrated model for surface flow, while hydrodynamic assumption is used for flow modeling in the present solution. Since shallow water assumption is valid in this test case, both Li et al. (2016) and the current model present good results. On the other hand, in some test cases, such as in the example in the section ‘Standing wave model’, where the hydrostatic assumption is not valid, the accuracy of the present solution is good enough.

Table 1

Root mean square errors (RMSE) for the present solution and Li et al. (2016) solution compared to experimental and theoretical results (cm)

 Point A, x = 0 Point B, x = 5 Point C, x = 18 Point D, x = 33.3 
RMSE (present solution compared to theoretical result)     
RMSE (present solution compared to experimental result)     
RMSE (Li et al. (2016) solution compared to theoretical result)     
RMSE (Li et al. (2016) solution compared to experimental result)     
 Point A, x = 0 Point B, x = 5 Point C, x = 18 Point D, x = 33.3 
RMSE (present solution compared to theoretical result)     
RMSE (present solution compared to experimental result)     
RMSE (Li et al. (2016) solution compared to theoretical result)     
RMSE (Li et al. (2016) solution compared to experimental result)     

Pressure drop in non-Darcy flow

In this example, uniform gas flow in a glass tube has been numerically studied. In the experimental model, a glass tube with 1-meter internal diameter and 10-meter length was considered. The incompressible gas flow has been inserted from one end of the tube while the boundary condition at the other end has been set as fixed pressure. Because of the 1D flow nature of this problem, the computing environment has been considered as a rectangle with width and thickness of 1 meter and a length of 20 meters in the present model. A schematic view of the model with a mesh length of 0.1 m accompanied with the boundary conditions is presented in Figure 9. The density and viscosity of the gas have been considered as and , respectively. Different values for flux from 1 to 8 have been inserted from the left boundary and pressure drop measurements were made in the tube. Nonlinear effects may be significant on flows in porous media if the Reynolds number is larger than 10. In this case, the nonlinear effects should be taken into account by adding the nonlinear terms to the linear Darcy's law. The Reynolds number in porous media has been calculated as following Equation (12):  
formula
(37)
Figure 9

Mesh plot and boundary conditions for pressure drop test.

Figure 9

Mesh plot and boundary conditions for pressure drop test.

By which, the Reynolds number has been estimated to be greater than 400 for different packing materials. Thus, linear Darcy equation was not valid in this problem and nonlinear drag forces were considered. Supported by experimental data of Ergun & Orning (1949), k and have been specified in a porosity variety from 0.37 to 0.76. The values of k and are presented in Table 2 for different packing materials.

Table 2

Material parameters for different packing materials

Packing material D50 (cm) Porosity (m2) (Permeability)  
1. Celite cylinders 0.64 0.37 1.37 × 10−8 0.0265 
2. Celite sphere 0.64 0.38 9.10 × 10−9 0.0251 
3. Celite cylinders 0.64 0.46 4.29 × 10−8 0.0386 
4. Celite spheres 0.64 0.47 2.31 × 10−8 0.0454 
5. Rasching rings 0.95 0.56 4.06 × 10−7 0.2898 
6. Rasching rings 0.95 0.62 5.69 × 10−8 0.0891 
7. Berl saddles 1.27 0.71 2.50 × 10−7 0.2775 
8. Berl saddles 1.27 0.76 4.33 × 10−7 0.3202 
Packing material D50 (cm) Porosity (m2) (Permeability)  
1. Celite cylinders 0.64 0.37 1.37 × 10−8 0.0265 
2. Celite sphere 0.64 0.38 9.10 × 10−9 0.0251 
3. Celite cylinders 0.64 0.46 4.29 × 10−8 0.0386 
4. Celite spheres 0.64 0.47 2.31 × 10−8 0.0454 
5. Rasching rings 0.95 0.56 4.06 × 10−7 0.2898 
6. Rasching rings 0.95 0.62 5.69 × 10−8 0.0891 
7. Berl saddles 1.27 0.71 2.50 × 10−7 0.2775 
8. Berl saddles 1.27 0.76 4.33 × 10−7 0.3202 
Pressure drop in this example has been calculated based on the following formula:  
formula
(38)
where is the difference between pressures of two ends of distance L, and u is the velocity at the middle of the two points where pressures have been measured.
Mass fluxes that enters from the left boundary were calculated according to the equation:  
formula
(39)
where is the total mass flux that enters from one end. and A are mass and the area, respectively. The area has been considered as 1.0 m2 in this test and is the time interval that has been considered to be 0.02 seconds. Figure 10 illustrates a comparison between the results from the present procedure (shown by lines) and that reported from the measurements of Oman & Watson (1944) (shown by symbols). As can be seen, the predictions are in good agreement with respect to the experimental results except that there are some small differences in packing 1. The dispersion in the measured data for packing 1 could be an indication for experimental error. Root mean square errors for different packing materials are presented in Table 3.
Table 3

Root mean square errors for pressure drop in different packing materials

Packing material 2. Celite sphere 1. Celite cylinders 4. Celite spheres 5. Rasching rings 3. Celite cylinders 6. Rasching rings 7. Berl saddles 8. Berl saddles 
RMSE 0.233 0.749 0.384 0.3 0.201 0.224 0.07 0.108 
Packing material 2. Celite sphere 1. Celite cylinders 4. Celite spheres 5. Rasching rings 3. Celite cylinders 6. Rasching rings 7. Berl saddles 8. Berl saddles 
RMSE 0.233 0.749 0.384 0.3 0.201 0.224 0.07 0.108 
Figure 10

Predicted and measured pressure drop vs. mass-flux relationships at high flow rates (lines show present solution and symbols show the measurements).

Figure 10

Predicted and measured pressure drop vs. mass-flux relationships at high flow rates (lines show present solution and symbols show the measurements).

Reflection and transmission of a solitary wave by a porous breakwater

In this example, a solitary wave hits a breakwater constructed with porous materials. The width of the breakwater is 0.5 m and its height is larger than the depth of the water plus the wave amplitude. Some part of the wave transmits through the barrier, some slowly damps, and the remaining reflects from the breakwater. Total domain length is 20 m, and the breakwater is located in the middle of the domain. The depth of the domain is 40 cm and the amplitude of the solitary wave is 3.8 cm before crossing the breakwater. The left boundary condition is specified as known elevations and the right boundary as closed boundary with zero velocities. A schematic view of the test case and its boundary conditions are shown in Figure 11. Transmission coefficient () is specified as the wave height ratio passed through the breakwater to the entire wave height and reflection coefficient () is defined as the ratio of the return wave height to the total wave height. Madson (1974) presented an analytical solution for this problem. In addition, Li et al. (2016) solved the problem numerically. More details may be found in the studies of Li et al. (2016) and Madson (1974).

Figure 11

Schematic view of porous breakwater and boundary conditions (not to scale).

Figure 11

Schematic view of porous breakwater and boundary conditions (not to scale).

Breakwater porosity is 0.5 and different values have been considered for mean grain diameter . Intrinsic permeability at each diameter has been calculated using Equation (6) with known value of porosity. Reynolds number in this problem is about 100, that is more than the range which is valid for linear assumption made by Darcy. Thus, non-Darcy terms should be considered in this case. Mesh length for this test has been considered as 4 cm and a schematic view of the mesh plot and dimensions are shown in Figure 12. Time intervals for the present solutions have been considered to be 0.02 s. The results determined in the present study are reported at five different and are compared in Figure 13 with the theoretical results and also the numerical results of Li et al. It should be mentioned that the numerical results of Li et al. (2016) were determined based on control volume method and non-Darcy shallow water assumption. As shown, the results are a very good match with the theoretical one. Root mean square errors between the present solution and analytical results are 0.015 and 0.014 for transmission and reflection coefficients, respectively. In addition, RMSE values for the results reported by Li et al. (2016) for the theoretical solution are 0.043 and 0.023 for transmission and reflection, respectively, which are larger than RMSE values obtained from the present model. The difference between the present solution and the results of Li et al. (2016) could be related to their assumption about flow in the vertical direction. They used hydrostatic assumption, which does not consider vertical velocities and accelerations, while in this model a hydrodynamic assumption is considered in the whole domain except for the part above the divider line.

Figure 12

Mesh plot of breakwater domain (not to scale).

Figure 12

Mesh plot of breakwater domain (not to scale).

Figure 13

Comparisons between transmission () and reflection () coefficients for the present solution, theoretical solution, and Li et al. (2016) numerical solution.

Figure 13

Comparisons between transmission () and reflection () coefficients for the present solution, theoretical solution, and Li et al. (2016) numerical solution.

Tidal oscillations in a lagoon system

This test is based on Ebrahimi et al. (2007) experimental set-up reported on the transmission of the tidal oscillations in a lagoon. In this test, a trapezoidal barrier built of non-cohesive sand separates a closed lagoon from the sea (Figure 14). Porosity and intrinsic permeability of the barrier are 0.3 and (Konga et al. 2010; Yuan et al. 2011), respectively. The left-hand side boundary condition is open sea that oscillates with amplitude of 60 mm with a period of 355 s, while right and bottom boundaries are closed boundaries with zero velocities. Figure 14 illustrates the dimensions of the experimental set-up and its boundary conditions. Points A, B, and C illustrated in this figure are the locations, where the solutions are presented and compared with the experimental observations. Point A and point B display water level fluctuations in the lagoon and in the open sea, respectively. On the other hand, point C is used to illustrate the velocity fluctuations in the open sea. More detailed descriptions of the measurements and data may be found in Ebrahimi et al. (2007) and Yuan et al. (2011). Konga et al. (2010) and Li et al. (2016) studied this experiment using finite volume–finite difference and finite volume methods, respectively. The current proposed numerical procedure has been applied to simulate the mentioned test case in which the grid size has been taken as 2 cm in triangular shape. Figures 15 and 16 demonstrate the comparisons of water surface elevation of measured values against the numerical results for point A and B, respectively, where an excellent agreement is shown in the comparisons. Moreover, the flow velocity that has been predicted by the present model for point C is compared with the experimental observation in Figure 17, where again an excellent match is achieved. Root mean square errors at points A, B, and C are calculated and presented in Table 4. It is worth mentioning that the Reynolds number in this problem is less than 10, and Darcy assumption is acceptable to a large extent (Bejan 1984).

Table 4

Root mean square errors at points A, B, and C

 Point A elev. Point B elev. Point C elev. 
RMSE 0.501 cm 0.181 cm 0.518 cm 
 Point A elev. Point B elev. Point C elev. 
RMSE 0.501 cm 0.181 cm 0.518 cm 
Figure 14

Experimental set-up for tidal induced fluctuation in a lagoon (Ebrahimi et al. 2007).

Figure 14

Experimental set-up for tidal induced fluctuation in a lagoon (Ebrahimi et al. 2007).

Figure 15

Comparison of the experimental results and the numerical evaluations of the current model at point A (see Figure 14).

Figure 15

Comparison of the experimental results and the numerical evaluations of the current model at point A (see Figure 14).

Figure 16

Comparison of the experimental results and the numerical evaluations of the current model at point B (see Figure 14).

Figure 16

Comparison of the experimental results and the numerical evaluations of the current model at point B (see Figure 14).

Figure 17

Comparison of the experimental results and the numerical velocities of the current model at point C (see Figure 14).

Figure 17

Comparison of the experimental results and the numerical velocities of the current model at point C (see Figure 14).

Comparison between the numerical solution determined based on the procedure proposed in this paper and the experimental data shows that the proposed FEM model based on the modified N-S equations provides reliable results in any complicated domain that consists of arbitrary distribution of porous and non-porous media in the presence of free surface fluctuations.

Standing wave model

A standing wave oscillation is discussed in this test case to emphasize the capability of the present model in taking the vertical velocities and accelerations into account. A rectangular basin with 10 m length and 10 m depth is considered as the domain of the model. First initial condition for the free surface is defined by:  
formula
where A is the wave amplitude, which is set to 0.1 m. The mesh length is considered as 0.1 m and the time step as 0.001 s, in this simulation. The density of fluid is set to 1,000 kg/m3 and viscosity is set to zero. Klingbeil & Burchard (2013) presented an analytical solution for this test considering an infinite number of modes rather than just one mode. The results of the present model, analytical solution, and hydrostatic model are compared in Figure 18. As can be observed, the hydrostatic assumption is not valid in this case and its result shows a large phase difference with the hydrodynamic solution. The main reason for this difference is that in hydrostatic assumption, vertical velocities and accelerations are not taken into account, and this test emphasizes the advantages of the present model among shallow water models. Moreover, the results of the present model fit well with theoretical results with root mean square error of m.
Figure 18

Comparison of water elevation for the present solution, analytical solution, and hydrostatic solution.

Figure 18

Comparison of water elevation for the present solution, analytical solution, and hydrostatic solution.

SUMMARY AND CONCLUSIONS

A numerical model has been developed to simulate 2D vertical surface water/subsurface flows, simultaneously. The governing equations have been derived in such a manner to be employed to both surface and subsurface domains facilitating application of a wider range of Reynolds numbers, so that the non-Darcy terms to solve the porous media are presented in the general form of equations. Fractional step method has been used for discretization in time, while finite volume–finite element methods have been applied for the space discretization of the governing equations. The whole domain has been divided into two parts, namely, the interior and the top layer subdomains. Modified Navier–Stokes equations have been considered for the interior subdomain, while shallow water type equations have been applied for the top layer. Triangular and linear gridding are used for the interior and top layer subdomains, respectively, where the equations are supposed to be solved in a coupled form, numerically. The procedure for the solution contains three major steps; first, the momentum equations are solved without the pressure terms and intermediate velocities are calculated; second, the Poisson equation which has been derived from continuity equation, is solved to obtain the pressure values; and third, using pressures of the top layer and using hydrostatic assumption, free surface elevations are calculated. In this way, free surface tracking has been achieved using the combination of hydrodynamic and hydrostatic formulation.

The results obtained from the model have been compared with some theoretical and experimental already published test cases to show the ability of this new procedure in modeling simultaneous surface and subsurface flows. The first two test cases were wave modeling in clear water and in porous media, respectively. Comparisons of the model results for these test cases have shown proper agreement, particularly in presenting phase property of the model. The third and the fourth examples were about non-Darcy flows with which the model results have been in very good agreement. The fourth test case together with the fifth one were about coupled surface and subsurface flows, in which the model results have shown to be reliable when an interaction of porous and non-porous domain with free surface is concerned. The last test was about emphasizing the difference of the present model and hydrostatic models. Results show the capability of the present model in the conditions that vertical velocities and accelerations are not ignorable and should be taken into account, properly.

ACKNOWLEDGEMENTS

The authors acknowledge valuable comments from two anonymous reviewers, which have led to improvement of the paper.

REFERENCES

REFERENCES
Amao
M.
2007
A Mathematical Model for Darcy Forchheimer Flow with Applications to Well Performance Analysis
.
MSc thesis
,
Texas Tech University
,
Lubbock, TX
,
USA
.
Balho
M.
,
Mikelić
A.
&
Wheeler
M. F.
2007
Polynomial filteration laws for law Reynolds number flows through porous medium
.
Transport in Porous Media
81
(
1
),
35
60
.
Bejan
A.
1984
Convection Heat Transfer
,
4th edn
.
John Wiley & Sons
,
New York
,
USA
.
Biot
M. A.
1941
General theory of three-dimensional consolidation
.
Journal of Applied Physics
12
(
2
),
155
164
.
Biot
M. A.
1956
Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range
.
The Journal of the Acoustical Society of America
28
(
2
),
168
178
.
Biot
M. A.
1962
Mechanics of deformation and acoustic propagation in porous media
.
Journal of Applied Physics
3
(
4
),
1482
1498
.
Braess
H.
&
Wriggers
P.
2000
Arbitraray Lagrangian Eulerian finite element analysis of free surface flow
.
Computer Methods in Applied Mechanics and Engineering
190
,
95
109
.
Breil
J.
&
Maire
P. H.
2006
A cell-centered diffusion scheme on two-dimensional unstructured meshes
.
Journal of Computational Physics
224
,
785
823
.
Brinkman
H. C.
1947
A calculation of viscous force extending by a flowing fluid on a dense swarm of particles
.
Applied Scientific Research
1
,
27
34
.
Casulli
V.
&
Cattani
E.
1994
Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow
.
Computers and Mathematics with Applications
27
(
4
),
99
112
.
Chegini
F.
&
Namin
M. M.
2012
A new approach to solving Poisson system for free surface nonhydrostatic flow simulations
.
International Journal for Numerical Methods in Fluids
70
,
562
577
.
Chen
F.
&
Chen
C. F.
1992
Convection in superposed fluid and porous layers
.
Journal of Fluid Mechanics
234
,
97
119
.
Chen
L.
,
Li
Y.
&
Thorpe
G. R.
1998
High-Rayleigh-number natural convection in an enclosure containing a porous layer
. In:
Proceedings of 11th Heat Transfer Conference
.
Korean Society of Mechanical Engineers
,
Seoul
, pp.
423
428
.
Cruchaga
M.
,
Battaglia
L.
,
Storti
M.
&
D'Elía
J.
2016
Numerical modeling and experimental validation of free surface flow problems
.
Computational Methods in Engineering
23
,
139
169
.
Ebrahimi
K.
,
Falconer
R. A.
&
Lin
B.
2007
Flow and solute fluxes in integrated wetland and coastal systems
.
Environmental Modeling and Software
22
,
1337
1348
.
Ergun
S.
&
Orning
A. A.
1949
Fluid flow through randomly packed columns and fluidized beds
.
Industrial and Engineering Chemistry
41
,
1179
1184
.
Forchheimer
P.
1901
Wasserbewegung durch Boden (Water movement through the ground)
.
Zeitschrift des Vereins Deutscher Ingenieure
45
,
1782
1788
.
Fromm
J. E.
1968
A method for reducing dispersion in convective difference schemes
.
Journal of Computational Physics
3
(
2
),
176
189
.
Griffiths
D. V. V.
&
Fenton
G. A.
1997
Three-dimensional seepage through spatially random soil
.
Journal of Geotechnical and Geoenvironmental Engineering
123
(
2
),
153
160
.
Harten
A.
,
Engquist
B.
,
Osher
S.
,
Chakravarthy
S.
1987
Uniformly high order accurate essentially non-oscillatory schemes, III
. In:
Upwind and High-Resolution Schemes
(
Hussaini
M. Y.
,
van Leer
B.
&
Van Rosendale
J.
, eds).
Springer-Verlag
,
Berlin, Heidelberg
,
Germany
, pp.
218
290
.
Heniche
M.
,
Secretan
Y.
,
Boudreau
P.
&
Leclerc
M.
2000
A two-dimensional finite element drying-wetting shallow water model for rivers and estuaries
.
Advances in Water Resources
23
,
359
372
.
Irmay
S.
1958
On the theoretical derivation of Darcy and Forchheimer formula
.
Eos, Transactions American Geophysical Union
39
,
702
707
.
Joseph
D. D.
,
Nield
D. A.
&
Papanicolaou
G.
1982
Nonlinear equation governing flow in a saturated porous medium
.
Water Resources Research
18
(
4
),
1049
1052
.
Kim
J.
&
Moin
P.
1985
Application of a fractional-step method to incompressible Navier-Stokes equations
.
Journal of Computational Physics
59
(
2
),
308
323
.
Larese
A.
,
Rossi
R.
&
Oñate
E.
2014
Finite Element Modeling of Free Surface Flow in Variable Porosity Media
.
Archives of Computational Methods in Engineering
,
Barcelona
,
Spain
.
Li
Y.
,
Yuan
D.
,
Lin
B.
&
Teo
F.
2016
A fully coupled depth-integrated model for surface water and groundwater flow
.
Journal of Hydrology
542
,
172
184
.
Lin
C. L.
,
Lee
H.
,
Lee
T.
&
Weber
L. J.
2005
A level set characteristic Galerkin fnite element method for free surface flows
.
Numerical Methods in Fluids
49
,
521
547
.
Liu
P. L. F.
&
Wen
J.
1997
Nonlinear diffusive surface waves in porous media
.
Journal of Fluid Mechanics
347
,
119
139
.
Liu
P.
,
Lin
P.
,
Chang
K. A.
&
Sakakiyama
T.
1995
Numerical modeling of wave interaction with porous structures
.
Journal of Waterway Port Coastal and Ocean Enginearing
125
(
6
),
322
330
.
Madson
O.
1974
Wave transmission through porous structures
.
Journal of the Waterways, Harbors and Coastal Engineering Division
100
,
169
188
.
Namin
M. M.
,
Lin
B.
&
Falconer
R. A.
2001
An implicit numerical algorithm for solving non-hydrostatic free-surface flow problems
.
International Journal for Numerical Methods in Fluids
35
,
341
356
.
Oman
A. O.
&
Watson
K. M.
1944
Pressure drops in granular beds
.
National Petroleum News
36
,
R795
R802
.
Osher
S.
,
Chakravarthy
S.
1986
Very high order accurate TVD schemes
. In:
Oscillation Theory, Computation and Methods of Compensated Compactness
(
Dafermos
C.
&
Ericksen
J. L.
, eds).
Springer
,
New York
,
USA
, pp.
229
274
.
Parlange
J. Y.
,
Stagnitti
F.
,
Starr
J. L.
&
Braddock
R. D.
1984
Free-surface flow in porous media and periodic solution of the shallow-flow approximation
.
Journal of Hydrology
70
,
251
263
.
Perot
J. B.
1993
An analysis of the fractional step method
.
Journal of Computational Physics
108
(
1
),
51
58
.
Press
W. H.
,
Teukolsky
S. A.
,
Vetterling
W. T.
&
Flannery
B. P.
1992
Numerical Recipes in Fortran 77
, Vol.
1
,
2nd edn
.
University of Cambridge
,
Cambridge
,
UK
.
Sahebkar
K.
&
Eskandari-Ghadi
M.
2017
Displacement ring load Green's functions for saturated porous transversely isotropic tri-material full-space
.
International Journal for Numerical and Analytical Methods in Geomechanics
41
(
3
),
359
381
.
Teng
H.
&
Zhao
T. S.
2000
An extension of Darcy's law to non-Stokes flow in porous media
.
Chemichal Engineering Science
55
,
2727
2735
.
Urumovic
K.
&
Urumovic
S. K.
2014
The effective porosity and grain size relations in permeability functions
.
Hydrology and Earth System System Sciences
11
,
6675
6714
.
Vafai
K.
&
Tien
C. L.
1981
Boundary and ineria effects on flow and heat transfer in porous media
.
International Journal of Heat and Mass Transfer
24
,
195
243
.
Venkatakrishnan
V.
&
Barth
T. J.
1989
Application of direct solvers to unstructured meshes for the Euler and Navier-Stokes equations using upwind schemes
. In:
27th Aerospace Sciences Meeting, AIAA
,
Reno, NV, USA
, pp.
89
0364
.
Vitousek
S.
&
Fringer
O. B.
2013
Stability and consistency of nonhydrostatic free-surface models using the semi-implicit θ-method
.
International Journal for Numerical Methods in Fluids
72
,
550
582
.
Wang
K. H.
&
Ren
X.
1993
Water waves on flexible and porous breakwaters
.
Journal of Engineering Mechanics
119
(
5
),
1025
1047
.
Wooding
R. A.
1957
Steady state free thermal convection of liquid in a saturated permeable medium
.
Physics of Fluids
12
,
273
285
.
Wu
J.
,
Yu
B.
&
Yun
M.
2008
A resistance model for flow through porous media
.
Transport in Porous Media
71
,
331
343
.
You
S.
&
Bathe
K. J.
2015
Transient solution of 3D free surface flows using large time steps
.
Computers and Structures
158
,
346
354
.
Yuan
L.
,
Xin
P.
,
Kong
J.
,
Li
L.
&
Lockington
D.
2011
A coupled model for simulating surface water and groundwater interactions in coastal wetlands
.
Hydrological Processes
25
,
3533
3546
.
Zeinkiewicz
O. C.
&
Taylor
R. L.
2000
The Finite Element Method: Solid Mechanics
, Vol.
2
,
5th edn
.
Butterworth-Heinemann, Oxford
.
Zienkiewicz
O. C.
&
Shiomi
T.
1984
Dynamic behavior of saturated porous media: the generalized biot formulation and its numerical solution
.
International Journal for Numerical and Analytical Methods in Geomechanics
8
(
1
),
71
96
.
Zienkiewicz
O. C.
,
Chang
C. T.
&
Bettess
P.
1980
Drained, undrained, consolidating and dynamic behaviour assumptions in soils
.
Géotechnique
30
(
4
),
385
395
.