The study of groundwater flow in fractured media can be distinct from that of flow in porous media due to the complex geologic configurations of fractured media. In modelling studies of flow in fractured media, the discrete fracture approach is preferable to the equivalent continuum approach because it has the potentiality to describe the fractures in more detail. However, the difficulty of using the discrete fracture approach lies in the explicit representation of the geometry of fractured media, which, in numerical modelling, refers to the discretization of the fractured media into a computational mesh. The composite element method (CEM), a novel numerical method developed from the finite element method, has an outstanding advantage in the discretization, so it is well suited for developing the numerical model of groundwater flow in fractured media.

The numerical model of groundwater flow in fractured media developed using the composite element method has the following features: (i) The fractures do not need to be discretized into specific elements, but are inserted into the porous matrix elements according to their geometric positions, resulting in the composite elements containing sub-elements of fracture and porous matrix; (ii) In the composite elements, the computational formulas are established for the flow through the fracture and porous matrix sub-elements respectively, which can be further integrated into the unified composite element computational formula through the connection at the fracture-matrix interface; (iii) The composite element computational formula will be automatically degenerated into the conventional finite element one if there is no inserted fracture in the element. Since the composite element computational formula has the same format as the finite element one, the finite element numerical techniques used to solve groundwater flow problems in the past can be carried over into the composite element method.

By defining multiple sets of hydraulic heads at the nodes of the composite element, the expression for estimating the hydraulic head in each sub-element is first established. Then, according to the groundwater flow differential equation and variational principle, the governing equations of flow in each sub-element are derived and subsequently assembled into the composite element governing equation. Finally, the developed numerical model is applied into different groundwater flow problems, including confined or unconfined flow, saturated or variably-saturated flow, steady-state or transient flow, different fracture distributions and boundary conditions of flow domain, in order to verify the accuracy and applicability of the model. The simulation results provide more theoretical basis for further understanding of groundwater flow characteristics in fractured aquifers.