A 2D vertical model for simulating surface and subsurface flows using finite element–finite volume methods

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.

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.

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.

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.

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).

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.

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).

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.

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.

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 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).

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).

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.

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.

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).

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.

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).

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.

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.

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