Numerical simulation of steady state groundwater flow field in discretely fractured rocks using the Green element method with the concept of multiple interacting normal flux

Numerical simulation is an effective tool for estimating the groundwater flow field in discretely fractured rocks (DFR). Unlike most numerical simulation methods that require the discretization of the model domain, boundary element method (BEM) is renowned of waiving the spatial discretization task but focusing on solving the integral form of the governing groundwater flow equation. However, for groundwater flow simulation in DFR, the solution obtained by BEM tends to have large error in the vicinity of fracture intersection. Therefore, a new numerical scheme, the green element method (GEM) is adopted in this study. GEM is built on the same mathematical background as BEM but turns the domain discretization back on as a necessary task. Using the second Green’s identity, GEM produces a general equation that applies to each grid block by integrating the governing equation. By making use of the singular characteristic of the Green’s function, GEM transforms the integral equation into a discretized system of equations with nodal head or nodal head gradient as unknowns. The cost of discretizing the model domain is compensated by the convenience of handling the heterogeneity of the medium. Conventional GEM manages the normal flux across a boundary segment by differentiating head values from 2 nodes in an individual grid block. This approximation overlooks the mechanism of normal flux as the exchange of fluid mass between grid blocks. To take this mechanism into consideration, a modified model of normal flux is proposed if the fracture plane is discretized into triangular elements. This model expresses the normal flux across a grid boundary segment in terms of the difference of head values in two grid blocks that are connected to this segment. For convenience, the head value at the centroid of a triangular element is used to calculate the normal flux. In other words, the unknowns of a triangular element are three nodal heads plus one centroidal head. Thus, the modified normal flux will be able to consider the interaction of all grid blocks that are connected to a target grid block. More importantly, the resulting global coefficient matrix is a square one and the system of equations is closed. The solution obtained from the closed system of equations will be exact but not a least-square approximated one. This modified GEM will be applied to simulate the steady state groundwater flow field in discretely fractured rocks.