Impact of porosity and permeability variations on water flow dynamics in heterogeneous subsurface porous media: A finite element analysis

Introduction

Water flow in heterogeneous subsurface porous media is crucial in many fields such as: management of groundwater resources, enhanced oil recovery, geothermal energy extraction and for prediction of contaminant transport [2, 5]. Both porosity ϕ (the volume fraction of void spaces in porous material) and permeability K (the ease of flow under a pressure gradient) are spatially variable as a result of depositional processes, diagenesis, and tectonic deformations [6] . These variations create preferential flow paths (channeling) that can significantly change breakthrough times, sweep efficiency and mass transport [1] . The two properties are related empirically by the Kozeny-Carman relation K ∝ ϕ 3

/(1 -ϕ) 2 which converts small differences in porosity to orders-of-magnitude differences in permeability [3] . Many of the current numer-ical models, however, do not capture this physical relationship, nor the nonlinear implications, and either assume homogeneity or treat porosity and permeability as uncoupled random fields. The finite element method (FEM) is particularly suited for heterogeneous problems be-cause it can directly incorporate spatially varying material properties [16, 10]. In recent years, some research efforts have focused on the coupled thermo-hydro-mechanical processes [14, 8], multi-phase flow [15] , and fractured media [13] using FEM. However, a systematic quantita-tive investigation on the effect of base porosity and spatial structure on the velocity field and pressure distribution is still missing in a well controlled heterogeneous environment. The present study fills this gap by creating a 2D finite element model which explicitly cou-ples porosity to permeability using the Kozeny-Carman relation, introduces a synthetic hetero-geneity field (sinusoidal layering + random noise), and conducts a parametric study of base porosity from 0.2 to 0.5. The main objective is to determine the influence of porosity and per-meability changes on flow rate of water and pressure distribution and to present the benchmark and physical insights for subsurface flow modeling.

Methodology

Governing equations

For steady, incompressible, single-phase flow in the Darcy regime, mass conservation and Darcy's law give:

where u is the Darcy velocity (m/s), K is permeability (m 2

), µ is dynamic viscosity (Pa•s), and p is pressure (Pa). Combining yields the elliptic pressure equation:

Domain and boundary conditions

The domain is a unit square Ω = [0, 1] × [0, 1]. Dirichlet pressure conditions are imposed on the left and right boundaries:

No-flow (Neumann) boundary conditions on top and bottom: ∂p/∂n = 0. This pressure difference pumps flow from the left to the right.

Heterogeneous porosity and permeability fields

To account for geological heterogeneity, the porosity field is assumed to be a sinusoidal layering with some added random noise: 𝜙 ( 𝑥, 𝑦 ) = 𝜙 0 + 𝐴 sin ( 2𝜋𝑦 ) + 𝜖 ( 𝑥, 𝑦 )

The amplitude is A = 0.2, and the base porosity ϕ 0 is varied from 0.2 to 0.5, and ϵ is a Gaussian random field with zero mean and standard deviation σ = 0.05 (interpolated from nodal values). Using Kozeny-Carman relation, the permeability is calculated as:

where K 0 = 1 × 10 -10 m 2 . In the assembly, K is scaled by 10 10 for numerical stability, but this factor is cancelled in the calculation of the velocity. Viscosity is set to µ = 1 × 10 -3 Pa • s (water at 20 • C).

Finite element discretization

The weak form of Eq. ( 2) is:

The domain is discretised with a symmetric mesh which is refined three times, resulting in a linear triangular (P1) element set consisting of 512 elements and 289 nodes. Basis functions are linear between the nodes. The stiffness matrix A with entries A ij = Ω(K/µ)∇ϕ i • ∇ϕ j dΩ is assembled using one-point quadrature. Eliminate (condense) Dirichlet conditions. A direct sparse LU solver (scipy.sparse.linalg.spsolve) is used to solve the reduced system A F,F p F = b mod . Element-wise Darcy velocity is computed after the nodal pressures are obtained as u e = -(K e /µ)(∇p) e and the nodal velocity is averaged for visualisation.

Parametric study design

The baseline parameters are ϕ 0 = 0.3, A = 0.2, σ = 0.05. Seven values of ϕ 0 are examined: 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50. For each, the permeability field is recomputed via Eq. ( 5), and the following quantities are recorded: Average velocity magnitude: u¯ = 1 ΣNelem u e ; Pressure drop: ∆p = p¯l eft -p¯r ight , where p¯l eft and p¯r ight are the averaged (node-wise) left and right boundary pressures, respectively. A power-law fit u¯ = aϕ b is obtained by non-linear least squares. All the other parameters (such as the random noise pattern) are kept the same to isolate the influence of ϕ 0 .

Code verification and validation

The numerical implementation was verified against the analytical solution for a homoge-neous medium (K constant), which reduces to Laplace's equation with solution p(x) = 1 -x. The L 2 error was 6.27 × 10 -1

. The model was further validated against sand column experiments [7] and the Mendeley permeability dataset [12] ; the maximum relative error was 8.6% (see [11] for full details). A mesh convergence study (16, 64, 256, 1024 elements) showed monotonic error reduction, confirming solution convergence.

Results

Baseline heterogeneous flow: Channeling

Figure 2 shows the baseline porosity field (ϕ 0 = 0.3): high-porosity bands (yellow, ϕ ≈ 0.45) at y ≈ 0.25 and low-porosity bands (dark blue, ϕ ≈ 0.2) at y ≈ 0.75. The differences are further enhanced by the Kozeny-Carman relation which shows that the permeability of the high porosity band is approximately 2 orders of magnitude higher than the low porosity band (Fig. 3). Figure 5 displays the velocity magnitude (colour) and direction (white arrows). Flow is intensely channeled into the high-permeability layer at y ≈ 0.25, with velocities up to five times the mean. Large regions of the domain (dark blue) have near-zero velocity. The coefficient of variation (CV) of the velocity magnitude is 0.48 (Table 1), which quantitatively shows strong flow focusing. The pressure field (Fig. 4) shows distorted isobars: dense contours in low-permeability regions (steep gradients) and sparse contours in the channel (shallow gradients), consistent with flux conservation.

Parametric study: Effect of base porosity (ϕ 0 )

The values of the average velocity u¯ and the pressure drop ∆p are given in Table 2 for ϕ 0 ranging between 0.2 and 0.5. In Figure 6 the plot of u¯ against ϕ 0 is shown. A power-law fit yields:

‾‾ 𝑢 = 1.1 × 10 3 𝜙 3.20 , 𝑅 2 = 0.995 (7) with exponent b = 3.20 ± 0.10. This is slightly higher than the theoretical exponent of 3 from the Kozeny-Carman relation, due to the (1ϕ) 2 term and the random noise. In figure 7 it is seen that for all values of ϕ 0 , the value of ∆p is constant (0.111). The slight departure from the theoretical value of 1.0 is because of averaging at the boundaries, the invariance shows that the global pressure gradient depends only on the Dirichlet boundary conditions and is independent of the internal heterogeneity.

Sensitivity to noise amplitude and correlation length

Since the amplitude of the noise and its correlation length are crucial, the system's sen-sitivity to them is analyzed. Additional sensitivity analyses (not part of the main parametric study but reported for ro bustness) showed that increasing the amplitude of random noise from σ = 0.02 to 0.10 increases the CV of velocity from 0.35 to 0.62 (Fig. 8), and increasing the spa-tial correlation length produces more continuous channels. These results verify that sub-scale heterogeneity plays a direct role in regulating the degree of flow channeling.

Discussion

Flow channeling and its physical origin

The strong channeling effect observed (CV = 0.48) is directly related to the nonlinear na-ture of the Kozeny-Carman relationship: a slight increase in porosity from 0.3 to 0.45 increases the permeability by a factor of about 6, which is approximately equal to ≈ 0.45 3 /0.55 2 . This amplification generates permeability contrasts that force the flow to seek the path of least resis-tance. This behavior has been observed in heterogeneous aquifers and oil reservoirs [1, 6]. The coefficient of variation of velocity (0.48) is a simple measure to quantify channeling intensity for different heterogeneity patterns.

Nonlinear velocity-porosity relationship

The power-law exponent 3.2 is close to the theoretical value of 3 from the Kozeny-Carman model. The slight excess is explained by the additional (1 -ϕ) 2 term and the random noise, which introduces additional variability. Practically, this exponent means that a 10% increase in porosity (e.g., from 0.3 to 0.33) increases the average flow rate by approximately 1.10 3.2 ≈ 1.36 or 36%. It is therefore important to characterise the porosity accurately in order to predict the flow rates in such a heterogeneous medium.

Invariance of the global pressure drop

The constant pressure drop (Fig. 7) is something intrinsic to elliptic equations with Dirichlet boundary conditions: the total head loss is not affected by the internal permeability distribution, but rather by the boundary potentials imposed on the solution. This outcome gives confidence that numerical implementations are correct in conserving flux and that any discrepancies in the computed value of ∆p would be due to incorrect handling of boundary conditions or solver tolerances.

Comparison with existing literature

The present results agree qualitatively with stochastic and multi-scale studies that report strong channeling and non-linear flux-porosity relationships [4, 9]. However, the use of a deter-ministic layered structure with controlled random noise allows a cleaner isolation of the effect of ϕ 0 . The exponent 3.2 falls within the range reported for granular porous media (3.0-3.5) [3] . The constant pressure drop observation corroborates theoretical results from homogenisation theory [2] .

Implications for practice

These results have direct practical implications:

Groundwater management: Monitoring for contaminant transport in low-velocity areas may not capture the transport in highvelocity channels; network designs need to consider channeling.

Enhanced oil recovery: Sweep efficiency is strongly dependent on small changes in poros-ity as shown by power-law exponent; polymer gels or profile modification should be aimed at the high-permeability layers.

Geothermal energy: Channeling leads to early thermal breakthrough; model can be used to optimise injection rates to reduce short-circuiting.

Upscaling: The constant ∆p property can be used to test coarsescale models; effective permeability must be a function of the entire distribution, not just the mean porosity.

Limitations

The study considers only steady-state, 2D, single-phase Darcy flow with isotropic perme-ability and idealized boundary conditions. The following effects are not accounted for: transient effects, multi-phase flow, geomechanics, and reactive transport. There were no uncertainty bounds provided, since a full Monte Carlo ensemble was not performed, but only sensitivity tests for the random noise amplitude and correlation length. Nevertheless, the parametric study provides a wellcontrolled benchmark for more complex models.

Conclusion

This study analysed the impact of porosity and permeability variations on water flow dy-namics in a heterogeneous subsurface porous medium using a finite element model. The main findings are:

1. Orders of magnitude differences in permeability result from small variations in porosity which are exacerbated by the Kozeny-Carman cubic relation, resulting in very strong flow channeling. The coefficient of variation of velocity magnitude is 0.48, and local velocities reach up to five times the mean.

2. The average velocity is related by a power-law to the base porosity u¯ ∝ ϕ 3.2 . This expo-nent is close to the theoretical value of 3 from the Kozeny-Carman model and quantifies the strong nonlinear sensitivity of flow rate to porosity.

3. The pressure drop along the global streamline is constant (∆p = 0.111) for all porosities from 0.2 to 0.5, which is a property of the elliptic flow equations and is a result of the imposed Dirichlet boundary conditions.

4. The results give quantitative measures (CV = 0.48, exponent 3.2) to validate other numer-ical models and upscaling strategies in heterogeneous porous media.

These results highlights the need to explicitly couple porosity-permeability field and the spatial heterogeneity in predictive models of subsurface flow. Future work should be extended to 3D geometries, transient conditions and multi-phase flow with uncertainty quantification.

Acknowledgement

AcknowledgementI thank Kisii University, Department of Mathematics and Actuarial Science, for providing computational facilities and a supportive research environment. I am grateful to Dr. Joash Kerongo for valuable discussions on the finite element implementation and to Robertus Urvakis for encouragement and support. No specific funding was received for this work.

References

1. Bachu S, Cuthiell D. Effects of core-scale heterogeneity on steady state and transient fluid flow in porous media: Numerical analysis. Water Resources Research. 1990;26(5):863-874.
2. Bear J. Dynamics of fluids in porous media. Elsevier; 1972.
3. Chapuis RP, Aubertin M. On the use of the Kozeny-Carman equation to predict the hy-draulic conductivity of soils. Canadian Geotechnical Journal. 2003;40(3):616-628.
4. de Barros FPJ, Guadagnini A, Bonazzi A. Probabilistic assessment of scalar transport un-der hydrodynamically unstable flows in heterogeneous porous media. Advances in Water Resources. 2024;185:104632.
5. Dullien FAL. Porous media: fluid transport and pore structure. 2nd ed. Academic Press; 1992.
6. Fetter CW. Contaminant hydrogeology. 2nd ed. Prentice Hall; 1999.
7. Li Z, Wan J, Huang K, Chang W, He Y. Effects of particle diameter on flow characteristics in sand columns. International Journal of Heat and Mass Transfer. 2017;104:533-536.
8. Pandey SN, Chaudhuri A, Kelkar S, Sandeep VR, Rajaram H. Investigation of perme-ability alteration of fractured limestone reservoir due to geothermal heat extraction using three-dimensional thermo-hydro-chemical (THC) model. Geothermics. 2014;51:46-62.
9. Pashchenko S, Kolesnikov A, Khachaturov R. Finite element modeling of filtration flows in heterogeneous porous media with log-normal and fractal permeability distributions. Computational Geosciences. 2025;29(1):89-106.
10. Reddy JN, Gartling DK. The finite element method in heat transfer and fluid dynamics. 3rd ed. CRC Press; 2010.
11. Rwito EK, Osogo AN, Bulinda MV, Monari FN. Numerical modeling of water flow dy-namics in heterogeneous subsurface porous media using the finite element method. Inter-national Journal of Statistical and Applied Mathematics. 2026;11(5):01-07.
12. Schro¨der L. Collection of Porous Media Datasets for Permeability Prediction (with Deep Learning) (Version 2) [Data set]. Mendeley Data. 2025.
13. Wang C, Winterfeld P, Johnston B, Wu YS. An embedded 3D fracture modeling approach for simulating fracture-dominated fluid flow and heat transfer in geothermal reservoirs. Geothermics. 2020;86:101831.
14. Wan J, Durlofsky LJ, Hughes TJR, Aziz K. Stabilized finite element methods for coupled geomechanics -reservoir flow simulations. SPE Reservoir Simulation Symposium. 2003.
15. Yoon HC, Zhou P, Kim J. Robust modeling of hysteretic capillary pressure and rela-tive permeability for two phase flow in porous media. Journal of Computational Physics. 2020;402:108915.
16. Zienkiewicz OC, Taylor RL. The finite element method. 5th ed. Butterworth-Heinemann; 2000.
17. How to cite this article: Rwito EK, Kerongo J. Impact of porosity and permeability variations on water flow dynamics in heterogeneous subsurface porous media: A finite element analysis. International Journal of Statistics and Applied Mathematics. 2026;11(9):138-143. doi: https://www.doi.org/10.22271/maths.2 026.v11.i9b.2696