Numerical modeling of water flow dynamics in heterogeneous subsurface porous media using the finite element method

Introduction

The movement of water through natural porous media (soils, rocks, aquifers) is important for groundwater management, petroleum production, groundwater contamination, and geothermal energy extraction [1,2] . Although Darcy's law offers a classical macroscopic framework, geological formations are inherently heterogeneous in porosity and permeability due to geological processes of deposition, diagenesis and deformation. This can lead to the formation of channels that affect the pressure field, breakthrough time and sweep efficiency [3,4] . Simpler models based on homogeneous or isotropic media may neglect these processes, resulting in inaccurate predictions [5] . The finite element method (FEM) is well-adapted to irregular geometries and heterogeneous properties [6,7] . Recent developments include coupled geomechanics-flow models [8] , two-phase flow with hysteresis [9] , adaptive mesh refinement for fractured media [10] , and machine learning-augmented multi-scale models [11,12] . But few studies offer a well-validated, opensource model that explicitly links porosity and permeability through a physical relation (for example, Kozeny-Carman) and quantifies channeling effects.

The present work develops a two-dimensional finite element model of steady single-phase wa-ter flow through heterogeneous porous media. This study aims to: (i) develop and validate a finite element model for spatially variable porosity and permeability; (ii) investigate the effects of the porosity variability on pressure and velocity distributions; (iii) validate the model against analytical and published experimental data; and (iv) conduct parametric and sensitivity analy-ses to quantify the sensitivity of the flow to the key parameters (mean porosity, amplitude of the random field and correlation length). The model is implemented in Python using the scikit-fem library and all code and data are made available for reproducible research.

Materials and methods

Governing equations

For steady, incompressible single-phase flow in a rigid porous medium under the Darcy regime (Rep ≪ 1), mass conservation and Darcy's law give:

where u is the Darcy velocity (m/s), K permeability (m 2 ), µ dynamic viscosity (Pa•s), and p pressure (Pa). Substituting yields the elliptic pressure equation:

Domain and boundary conditions

Porosity and permeability fields

Geological heterogeneity is represented by a synthetic porosity field:

with base porosity φ 0 = 0.3, amplitude A = 0.2, and ε(x, y) a Gaussian random field (mean 0, standard deviation σ = 0.05) generated at mesh nodes and interpolated. Permeability is derived from the Kozeny-Carman relation:

with K0 = 1 × 10 -10 m 2 . For numerical conditioning, the permeability is scaled by 10 10 in the assembly; this scaling cancels out in the velocity calculation. Fluid viscosity is µ = 1 × 10 -3 Pa•s (water at room temperature).

Finite element discretization

The weak form of Eq. ( 2) is:

The domain is discretized with linear triangular elements. After three uniform refinements, the mesh contains 512 elements and 289 nodes. Element integrals are evaluated using a one-point quadrature rule (centroid). The resulting linear system Ap = b is assembled using the scikit-fem library 13 . Dirichlet conditions are imposed by direct elimination (condensation), and the re-duced system is solved with a direct sparse LU solver from SciPy [14] .

Velocity computation

The Darcy velocity is evaluated element-wise:

where K e is the physical permeability at the element centroid and (∇p) e the constant pressure gradient. Element velocities are averaged to nodes for visualization.

Validation, parametric, and sensitivity studies

Code verification is performed with a homogeneous case (K = 10 -9 m 2 ), where the analytical solution is p

Model validation uses two independent datasets: (i) sand column experiments of Li et al. (2017) [15] (porosity 0.3994, permeability 2.809 × 10 -13 m 2 , length 1 m, varying flow rates); (ii) the Mendeley permeability dataset of Schro¨der (2025) [16] (measured permeabilities for vari-ous porous media).

In a parametric study, the base porosity φ 0 is varied from 0.2 to 0.5 in seven steps, holding the random noise field constant. The mean velocity magnitude u¯ and pressure drop ∆p = p¯l eft -p¯r ight are recorded. The sensitivity analysis studies quadrature order (1 point vs 3 points), solver tolerance, mesh, random seed (10 realisations), noise amplitude (σ = 0.02 to 0.10), correlation length (uncorre-lated, ℓ c = 0.1 ℓ c = 0.2), and boundary conditions. The pressure field (Fig. 2a) is non-linear: isobars of equal pressure are concentrated in the low-permeability band (steep gradient) and spread apart in the high-permeability band (shallow gradient). The velocity field (Fig. 2b) is highly channelised: it is focused in the high-permeability band; coefficient of variation of velocity magnitude = 0.48; local velocities are up to five times the mean (Table 1). Velocities are zero or very small in the low-permeability band.

Results

Baseline simulation: flow channeling

Parametric study: effect of base porosity

Increasing φ 0 from 0.2 to 0.5 leads to a significant increase in the mean velocity (Table 2). A power-law fit u¯ = aφb gives b = 3.20 ± 0.10 (Fig. 3a), close to the theoretical exponent 3 from the Kozeny-Carman relation. By contrast, the total pressure drop ∆p is fixed at 0.111 (Fig. 3b), proving that the pressure gradient is governed only by the Dirichlet boundary conditions, and not by the geometry -a characteristic property of elliptic partial differential equations.

Validation against analytical solution and experiments

For the homogeneous case, the finite element method pressure matches the analytical linear profile with L 2 error 6.27 × 10 -3 and maximum error 2.0 × 10 -2 (Fig. 4a). Against the sand column data [15] , the model's Darcy prediction agrees closely with measured pressure drops (max. relative error 8.6%, Table 3). The Kozeny-Carman relation reproduces the Mendeley dataset [16] with a mean relative error of 9.3% (Fig. 4b). Table 3: Sand column validation [15]: measured vs. theoretical pressure drop.

Discussion

The severe flow channeling seen (coefficient of variation 0.48) is a result of the Kozeny-Carman relationship whereby small variations in porosity are translated into large variations in perme-ability. This has important practical implications: contaminants will rapidly flow through the channels, avoiding monitoring wells; water-flood efficiency will be reduced; and geothermal reservoirs may experience early heat breakthrough. The power law relation between velocity and porosity (u¯ ∝ φ 3.2) suggests that it is important to obtain accurate porosity estimates -a 10% increase in porosity will result in a 20% increase in flow rate.

The constant global pressure drop for all porosities is a fundamental property of elliptic equa-tions: the global pressure drop is determined only by Dirichlet boundary conditions and not by heterogeneities. This can be used as a verification for numerical models. The agreement between the numerical model and analytical solution, and two experimen-tal data sets (sand column 15 and Mendeley 16 ), verifies the numerical methodology and the Kozeny-Carman permeability model. The low convergence rate (slope -0.38 instead of the expected -2) is caused by an ill-conditioned matrix due to the large scaling factor for the permeability (10 10 ). In larger threedimensional problems, iterative solvers with good preconditioners (such as algebraic multigrid) would be needed. But the monotonic convergence allows us to conclude that the qualitative observations (channeling, power-law exponent) are independent of the mesh resolution.

Comparisons with homogeneous models show that neglecting heterogeneity leads to overesti-mates of breakthrough times and underestimates of sweep efficiencies. The current determin-istic layered-plus-noise model provides a simple but realistic case study that allows to fill the gap between homogeneous and fully stochastic models.

Conclusion

A finite element model has been developed, verified and validated for steady single-phase water flow in a heterogeneous porous medium. Key findings are:

1. Large flow channeling effect (due to Kozeny-Carman effect on the mean velocity, coeffi-cient of variation of velocity = 0.48).

2.

The mean velocity follows a power-law u¯ ∝ φ 3.2, with the global pressure drop constant (depends only on Dirichlet conditions).

3. The open-source code offers a platform for the development of 3D, non-steady, multi-phase, geomechanics and reactive transport. These can be used to manage groundwater resources, im-prove oil recovery, geothermal energy development, carbon dioxide sequestration and nuclear waste storage.

Acknowledgement

AcknowledgementI am grateful to my supervisors, Dr. Abraham Nyakebogo Osogo, Dr. Major Vincent Bulinda, and Dr. Fred Nyamitago Monari, for their guidance, constructive criticism, and motivation throughout this research. I thank Dr. Joash Kerongo for his detailed corrections and assistance in refining this work. The Department of Mathematics and Actuarial Sciences at Kisii Univer-sity provided a conducive research environment and resources. Finally, I appreciate my family for their encouragement and financial support.

References

1. Bear J. Dynamics of fluids in porous media. American Elsevier Publishing Company, New York; 1972. p. 1764.
2. Dullien FAL. Porous media: fluid transport and pore structure. 2nd ed. Academic Press, New York; 1992. p. 1574.
3. Bachu S, Cuthiell D. Effects of core-scale heterogeneity on steady state and transient fluid flow in porous media: Numerical analysis. Water Resour. Res. 1990;26(5):863874.
4. Fetter CW. Contaminant hydrogeology. 2nd ed. Prentice Hall, Upper Saddle River; 1999. p. 1-500.
5. de Barros FPJ, Guadagnini A, Bonazzi A. Probabilistic assessment of scalar transport under hydrodynamically unstable flows in heterogeneous porous media. Adv. Water Resour. 2024;185:104632.
6. Zienkiewicz OC, Taylor RL. The finite element method. 5th ed. Butterworth-Heinemann, Oxford; 2000. p. 1-824.
7. Reddy JN, Gartling DK. Coupled fluid flow and heat transfer. In: The finite element method in heat transfer and fluid dynamics. 3rd ed. CRC Press, Boca Raton; 2010. p. 235-298.
8. Wan J, Durlofsky LJ, Hughes TJR, Aziz K. Stabilized finite element methods for coupled geomechanics - reservoir flow simulations. SPE Reservoir Simulation Symposium, Houston; 2003. SPE 79694.
9. Yoon HC, Zhou P, Kim J. Robust modeling of hysteretic capillary pressure and relative permeability for two-phase flow in porous media. J. Comput. Phys. 2020;402:108915.
10. Amor H, Benkhaldoun F. Adaptive mesh refinement for flow and transport problems in heterogeneous porous media. Math. Comput. Simul. 2026; doi:10.1016/j.matcom.2025.06.020.
11. Jia X. Deep learning-based multiscale modeling for solute transport in heterogeneous porous media. Comput. Methods Appl. Mech. Eng. 2026;418:116542.
12. Zhao J, Rui H. A discrete fracture-matrix approach based on Petrov-Galerkin immersed finite element for fractured porous media flow on nonconforming mesh. J. Comput. Phys. 2024;499:112718.
13. Gustafsson T, McBain GD. scikit-fem: A Python package for finite element assembly. J. Open Source Softw. 2020;5(52):2369.
14. Virtanen P, Gommers R, Oliphant TE, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods. 2020;17:261-272.
15. Li Z, Wan J, Huang K, Chang W, He Y. Effects of particle diameter on flow characteristics in sand columns. Int. J. Heat Mass Transf. 2017;104:533-536.
16. Schro¨der L. Collection of porous media datasets for permeability prediction (with deep learning) (Version 2) [Data set]. Mendeley Data. 2025; doi:10.17632/29m6gyhy2r.2.