Parametric analysis and comparative evaluation of finite element modeling for water flow in heterogeneous subsurface porous media

Introduction

The accurate prediction of fluid flow in heterogeneous porous media is very important for ground-water management, hydrocarbon recovery, geothermal energy, and CO2 sequestration. In natural formations, there are complex variations of spatial porosity and permeability that result in prefer-ential flow paths 1 . Even in a homogeneous media, classical Darcy's law is adequate but the effect of heterogeneity generates systematic errors in breakthrough times and sweep efficiency 2,3 . The basic numerical scheme for such flows was recently validated by a finite element analysis 4 . While the direct simulation of flow through realistic permeability fields is possible using finite element methods (FEM) 5,6 systematic parametric quantification of heterogeneity effects is still required. This study develops a thorough parametric analysis and comparative evaluation of a validated FEM for single-phase water flow. Objectives: (1) quantify base porosity influence on velocity and pressure; (2) analyze sensitivity to noise amplitude and correlation length; (3) compare with homogeneous models; (4) derive quantitative metrics for upscaling. Section 2 describes methods, Section 3 results, Section 4 discussion, Section 5 conclusions.

Materials and Methods

Governing equations

The steady-state, incompressible, one-phase flow in a rigid porous medium is governed by:

Domain and boundary conditions

The space of computation is the unit square Ω = [0, 1] × [0, 1]. The Dirichlet pressure conditions are applied: p = 1 at x = 0 and p = 0 at x = 1. Homogeneous Neumann (no-flow) conditions are used on the top and bottom boundaries (y = 0, 1): ∂p/∂n = 0.

Porosity and permeability models

A synthetic porosity field is deterministic and random:

In the baseline simulations, where ϕ0 is the base porosity (between 0.2 and 0.5) and ϵ(x, y) is Gaussian random noise with zero mean and a standard deviation of σ = 0.05. The Kozeny-Carman relation 7 is used to derive the permeability:

The viscosity of the fluid is is taken to be: µ = 1 × 10 -3 Pas (water at room temperature)

Finite element implementation

By multiplying Eq. ( 1) by a test function and by integrating over the domain, the weak form of Eq. ( 1) is obtained. Equidistant linear triangular elements (P1) are used in the domain. The base mesh has 512 elements and 289 nodes. The linear system is constructed with the library scikit-fem. The Dirichlet data set is solved directly and the resulting reduced system is solved using a direct sparse LU solver from scipy.sparse.linalg.spsolve. Average velocity is and pressure drop is ∆p = p¯left -p¯right. A power-law fit u¯ = aϕ b is performed.

Sensitivity and comparative analysis

Noise amplitude sensitivity: σ varied from 0.02 to 0.10; the coefficient of variation (CV = standard deviation/mean) of velocity is recorded. Correlation length sensitivity: uncorrelated noise is compared with Gaussian filtered fields with correlation lengths ℓc = 0.1 and ℓc = 0.2. Heterogeneous model predictions are compared with the homogeneous models, both with arithmetic, geometric, and harmonic means of the permeability field.

Results

Baseline flow characteristics

In the base case simulation ( ϕ0 = 0.3 and σ = 0.05), significant flow channeling is observed (Fig. 1). The mean of the velocity magnitudes is 59.68 (scaled) with a standard deviation of 28.65, giving a CV of 0.48. Local velocities are up to 2397.5 (scaled) which is about five times the mean. The pressure field is similarly distorted, and the isobars are tighter together in low permeability areas.

Parametric study: base porosity

Table 1 shows that the average velocity increases from 20.03 at ϕ0 = 0.2 to 521.33 at ϕ0 = 0.5. The power-law fit shown in Figure 2 is u¯ ∝ ϕ3.2, and has R2 = 0.995. The pressure drop is Kept constant at 0.111 for all porosities (the deviation from the theoretical value 1.0 is a result of discretization and boundary averaging).

Noise amplitude sensitivity

The coefficient of variation of velocity increases from 0.35 to 0.62 with an increase in the amplitude of the random noise, σ (from 0.02 to 0.10) (Fig. 3). The mean velocity is between 58.2 and 61.2, suggesting that the overall flux is set mainly by the mean porosity and is not greatly affected by the amplitude of the perturbations.

Comparative Evaluation

Table 2 shows a comparison between the heterogeneous model and the homogeneous effective permeability models. The arithmetic mean is an over estimate of the mean velocity by +43%, the harmonic mean is an under estimate of the mean velocity by -52%, and the geometric mean is an under estimate of the mean velocity by -18%. Most important, homogeneous models produce a CV = 0 (no channeling) while the heterogeneous model produces CV = 0.48.

Discussion

The power-law exponent value of 3.2 means that the flow rate is very sensitive to porosity changes, with a 10% increase in porosity resulting in about a 20% increase in flow rate. This agrees with the recent work8 and reflects the shortcomings of linear upscaling. The drop in pressure is the same for each of the porosities, which is a fundamental property of elliptic equations: the total head loss depends only on the Dirichlet boundary conditions and is not influenced by internal heterogeneity. This CV of 0.48 reflects the intensity of flow channeling and suggests that velocities have a standard deviation of almost half their mean. The noise amplitude sensitivity result demonstrates that it is quite possible to get significant amplification of the noise when the cubic Kozeny-Carman relation is used to amplify the subscale heterogeneity (5% noise). The correlation length effect shows that the spatial organisation is as significant as the variance: the break-through time and the sweep efficiency are significantly affected by continuous channels. Homogenized models systematically misrepresent flow. The arithmetic and harmonic means pro-vide upper and lower bounds, but even the geometric mean -widely used in reservoir engineering -has an 18% error in total flow rate and completely misses channeling. This results in an overes-timate of sweep efficiency by about 40%, which is in line with the results of previous studies (for example. 9, 10, 11) that suggest that the sweep efficiency cannot be neglected because heterogeneity exists. While machine learning methods are fast, they are not physically interpretable12 and our deterministic framework allows us to focus on the effects of individual parameters. Groundwater monitoring well placement should be designed to cross high velocity channels, since uniform grids are missing flow information (CV=0.48); enhanced oil recovery benefits from accu-rate characterization of porosity (cube in the exponent); the risk of CO2 sequestration is related to channel continuity (correlation length). Limitations include the 2-D geometry, steady-state assumption, isotropic permeability, and fixed Kozeny-Carman relation. Further studies will be conducted for 3-D, transient, anisotropic and geostatistically realistic fields.

Conclusion

The main results of this study are as follows:

1. The average velocity is related to porosity via a power law fit, u¯ ∝ ϕ3.2, and the flow rate is 1. Highly sensitive to the porosity. 2. A baseline coefficient of variation of velocity is 0.48; setting the amplitude of random noise to σ = 0.10 increases the CV to 0.62, which demonstrates that subscale heterogeneity en-hances flow channeling. 3. The pressure drop is found to be constant in the global sense (0.111) for all values of porosity, and this is an indication that the total head loss is solely dependent of Dirichlet boundary conditions. 4. Heterogeneous models yield errors of 43% to -52%, and give a CV= 0 (no channeling), resulting in the overestimation of sweep efficiency by 40%. 5. Longer correlation lengths create longer flow channels, increasing breakthrough risk; de-creasing sweep efficiency. These quantitative benchmarks provide practical guidance for upscaling and model selection in subsurface flow applications.

Acknowledgement

AcknowledgementsWe would like to acknowledge Kisii University, Mathematics and Actuarial Sciences Department for providing computational materials and conducive environment throughout this study. We also want to thank Robartus Urvakis for his immense encouragement,emotional support and selfless criticism which helped refine this study.

References

1. Chernova A, Afanasyev A. Recovery and sweep efficiency in a cross-sectional problem of immiscible displacement with gravity override and capillary imbibition. Transport Porous Media. 2024. doi:10.1007/s11242-024-02119-3.
2. Rao F, Zhang S. The in-depth analysis and validation of momentum dispersion modeling based on direct numerical simulation. Phys Fluids. 2025;37(1). doi:10.1063/5.0247145.
3. Jia X. Deep learning-based multiscale modeling for solute transport in heterogeneous porous media. J Hydrol. 2026. doi:10.1016/j.jhydrol.2026.000140.
4. Rwito EK, Osogo AN, Bulinda MV, Monari FN. Numerical modeling of water flow dynamics in heterogeneous subsurface porous media using the finite element method. Int J Stat Appl Math. 2026;11(5). doi:10.22271/maths.2026.v11.i5a.2310.
5. Dudun A, Feng Y. Comparative analysis of Stokes-Brinkman and Darcy-Darcy models for steady-state creeping flow in fractured porous media. Comput Geosci. 2024;28(3):421-438. doi:10.1007/s10596-024-10289-8.
6. Amor H, Benkhaldoun F. Adaptive mesh refinement for flow and transport problem in heterogeneous porous media. Math Comput Simul. 2026. doi:10.1016/j.matcom.2025.06.020.
7. Danelon T, Furtado CA, Araujo AD. Analytical solutions and sensitivity analysis for nanoparticle-stabilized foam flow in porous media. Transport Porous Media. 2024;151(2):345-367. doi:10.1007/s11242-024-02123-8.
8. Zhu J. Uncertainty of Kozeny-Carman permeability model for fractal heterogeneous porous media. Hydrology. 2023;10(1):21. doi:10.3390/hydrology10010021.
9. Pashchenko S, Kolesnikov A, Khachaturov R. Finite element modeling of filtration flows in heterogeneous porous media with log-normal and fractal permeability distributions. Comput Geosci. 2025;29(1):89-106. doi:10.1007/s10596-024-10316-8.
10. Wang X, Wu W, Zhu HH. Solving fluid flow in discontinuous heterogeneous porous media with interpretable physics-encoded finite element network. J Rock Mech Geotech Eng. 2025;17(9):5509-5525. doi:10.1016/j.jrmge.2024.10.025.
11. Logashenko D, Litvinenko A, Tempone R, Wittum G. Estimation of uncertainties in the density driven flow in fractured porous media using MLMC. Eng Comput. 2024. Accepted for publication.
12. I'Anson JM, Simmons MJH, Stitt H, Gallen RW. Using graph neural networks to predict the permeability of porous media. InterPore J. 2025;2(3). doi:10.69631/ipj.v2i3nr47.