A finite element method using a bounded auxiliary variable for solving the Richards equation

Basic Information

• Paper ID: 2510.13012
• Title: A finite element method using a bounded auxiliary variable for solving the Richards equation
• Authors: Abderrahmane Benfanich (University of Ottawa), Yves Bourgault (University of Ottawa), Abdelaziz Beljadid (University of Ottawa & Mohammed VI Polytechnic University)
• Classification: math.NA cs.NA physics.comp-ph
• Publication Date: October 14, 2025 (arXiv preprint)
• Paper Link: https://arxiv.org/abs/2510.13012

Abstract

The Richards equation is a nonlinear elliptic-parabolic equation widely used for modeling infiltration in porous media. This paper develops a finite element method for solving the Richards equation by introducing a new bounded auxiliary variable to eliminate unbounded terms in the weak formulation. The method employs a semi-implicit discretization scheme and uses Newton's method to solve the resulting nonlinear systems. The approach eliminates the need for regularization techniques and offers advantages in handling dry and fully saturated regions. In the proposed technique, non-overlapping Schwarz domain decomposition methods are employed to simulate infiltration in layered soils. Numerical experiments validate the effectiveness of the proposed method, including flow modeling in completely dry fiber sheets, cases with both fully saturated and dry regions, and infiltration problems in layered soils.

Research Background and Motivation

Problem Definition

The Richards equation describes unsaturated flow in porous media and is a fundamental model in hydrogeology, environmental engineering, and soil science. The equation exists in multiple formulations:

1. Pressure head form (Ψ): Suitable for fully saturated and partially unsaturated flow, but numerically unstable as S→0
2. Saturation form (S): Performs better in dry soils but encounters numerical difficulties in saturated regions (S=1)
3. Mixed form (S,Ψ): Improves mass conservation properties but increases complexity

Limitations of Existing Methods

1. Primary variable switching methods: May produce non-physical solutions at saturated/unsaturated interfaces
2. Regularization techniques: Require additional modeling parameter selection
3. Traditional finite element methods: Encounter numerical difficulties when handling completely dry (S=0) and fully saturated (S=1) regions

Research Motivation

Develop a unified numerical method that can:

• Avoid variable switching and regularization
• Stably handle dry and saturated regions
• Maintain good mass conservation properties
• Apply to heterogeneous porous media

Core Contributions

1. Propose a new u-formulation: Introduce a bounded auxiliary variable u to eliminate unbounded terms in the weak formulation
2. Develop a stable numerical scheme: Combine semi-implicit time discretization with Newton iteration
3. Implement domain decomposition methods: Address heterogeneous media problems such as layered soils
4. Verify method effectiveness: Demonstrate stability and accuracy through multiple numerical experiments

Methodology Details

Problem Formulation

Solve the Richards equation: $\frac{\partial\theta}{\partial t} + \nabla \cdot q = 0$

where θ is the water content and q is the water flux according to the Darcy-Buckingham law: $q = -K_s(x)K_r(x,S)\nabla(\Psi(x,S) + z)$

Construction of the New Variable u

For all models, the derivative of the Leverett function can be written as: $J'(S) = C(x)S^{-a}(1-S^c)^{-b}$

Define the auxiliary variable: $u = \int_0^S (1-s^c)^{-b} ds$

Expressed using the incomplete Beta function: $u = \frac{1}{c}B(S^c, \frac{1}{c}, 1-b)$

Richards Equation in u-Form

In homogeneous media, the u-form is: $\phi \frac{\partial S(u)}{\partial t} - \nabla \cdot \left(K_sK_r(S(u))\left(h_{cap}CS(u)^{-a}\nabla u + e_z\right)\right) = 0$

Variational Formulation

Find u ∈ L²(I;H¹(Ω)) such that: $\int_\Omega \phi \frac{\partial S(u)}{\partial t} v dx + \int_\Omega K_sK_r(S(u))\left(h_{cap}S(u)^{-a}\nabla u + e_z\right) \cdot \nabla v dx + \text{boundary terms} = 0$

Time Discretization

Employ a semi-implicit Euler scheme: $\int_\Omega \phi \frac{S(u_h^{n+1}) - S(u_h^n)}{\Delta t} v_h dx + \int_\Omega K_sK_r(S(u_h^n))\left(h_{cap}S(u_h^n)^{-a}\nabla u_h^{n+1} + e_z\right) \cdot \nabla v_h dx = 0$

Newton Iteration

For the nonlinear term S(u^{n+1,m+1}), use first-order Taylor expansion: $S(u^{n+1,m+1}) = S(u^{n+1,m}) + \frac{\partial S}{\partial u}(u^{n+1,m})(u^{n+1,m+1} - u^{n+1,m}) + O((u^{n+1,m+1} - u^{n+1,m})^2)$

Domain Decomposition Method

For layered soils, employ the non-overlapping Schwarz method:

• Solve the Richards equation independently in each subdomain
• Couple through interface conditions: pressure head continuity and flux conservation
• Use Robin-type transmission conditions for iterative solving

Experimental Setup

Numerical Implementation

• Software platforms: FreeFEM++ and FENICSx
• Linear solvers: UMFPACK, PETSc, MUMPS
• Mesh: Uniform triangular mesh with linear finite elements (k=1)

Evaluation Metrics

L² error: $E_{h,\Delta t}(T)^2 = \|u(\cdot,T) - u_{h,\Delta t}(\cdot,T)\|_{L^2(\Omega)}^2$

Convergence rates:

• Spatial convergence rate: $p = \log_2(E_{h,\Delta t}(T)/E_{h/2,\Delta t}(T))$
• Temporal convergence rate: $q = \log_2(E_{h,\Delta t}(T)/E_{h,\Delta t/2}(T))$

Test Cases

1. 1D completely unsaturated fiber sheet: Verify stability under dry conditions
2. 1D saturated-unsaturated mixed region: Test capability to handle discontinuous initial values
3. Manufactured solution verification: Compute convergence rates
4. 2D fully saturated flow: Compare with Hydrus software
5. Layered soil problem: Verify domain decomposition method
6. Heterogeneous media inclusion: Test handling of complex geometries

Experimental Results

Main Results

1D Fiber Sheet Test

• Mesh: 500 uniform grid points
• Time step: Δt = 1s
• Results: u-form and S-form results are completely consistent and agree well with experimental data

Saturated-Unsaturated Mixed Region

• Mesh: 5000 uniform grid points
• Time step: Δt = 10⁻⁵ hours
• Newton convergence: Average 1-5 iterations
• Stability: Solution maintains monotonicity with no numerical oscillations

Manufactured Solution Convergence Analysis

Spatial Convergence (fixed Δt = 5×10⁻⁵):

Temporal Convergence (fixed h = 1×10⁻³):

2D Saturated Flow Comparison

• Mesh: 95,142 triangular elements, 47,972 vertices
• Time step: Δt = 0.001 hours
• Results: Highly consistent with Hydrus software results

Layered Soil Test

• Mesh: 89,306 triangular elements, 45,064 vertices
• Time step: Δt = 0.05 hours
• Convergence: Tolerance 10⁻³, maximum 10 iterations, λ = 25
• Results: Consistent with results in reference 20

Method Comparison

In the heterogeneous media inclusion test:

• u-form: Converges stably even when ε=0, average 4 iterations
• Ψ-form: Requires >1000 iterations when ε=10⁻⁵, residual tolerance still 0.04

Related Work

Main Research Directions

1. Variable switching methods: Diersch & Perrochet (1999), Forsyth et al. (1995)
2. Regularization techniques: Schweizer (2007), Pop & Schweizer (2011)
3. Mixed formulations: Celia et al. (1990), Maina & Ackerer (2017)
4. Domain decomposition: Lions (1990), Manzini & Ferraris (2004)

Advantages of This Work

1. Unified treatment: No need for variable switching or regularization
2. Numerical stability: Stable across the entire range S∈0,1
3. Physical consistency: Maintains pressure head continuity
4. Computational efficiency: Rapid Newton convergence

Conclusions and Discussion

Main Conclusions

1. Effectiveness of u-form: Successfully eliminates unbounded terms in traditional formulations
2. Numerical stability: Maintains stability in both completely dry and fully saturated regions
3. Convergence performance: Second-order spatial convergence, first-order temporal convergence
4. Practicality: Applicable to complex heterogeneous media problems

Limitations

1. Model dependency: Requires the condition $\lim_{S\to 0} K(S)S^{-a} < \infty$
2. Nonlinear solving: Requires Newton iteration in saturated regions
3. Interface handling: Heterogeneous media require domain decomposition techniques
4. Parameter sensitivity: Robin parameter λ in domain decomposition requires tuning

Future Directions

1. Theoretical analysis: Rigorous numerical analysis and convergence proofs
2. Multiphysics coupling: Extension to solute and heat transport
3. Adaptive techniques: Develop adaptive mesh and time-stepping strategies
4. Parallel computing: Large-scale parallel implementation

In-Depth Evaluation

Strengths

1. Strong innovation: Clever construction of the u-variable with solid theoretical foundation
2. Numerical stability: Avoids numerical difficulties of traditional methods
3. Comprehensive experiments: Covers cases from 1D to 2D, homogeneous to heterogeneous
4. Engineering practicality: Comparison with commercial software Hydrus enhances credibility

Weaknesses

1. Insufficient theoretical analysis: Lacks rigorous convergence and stability proofs
2. Incomplete computational complexity analysis: Insufficient detail on Newton iteration computational costs
3. Parameter selection guidance: Lacks theoretical guidance for Robin parameter λ selection
4. Missing 3D verification: Lacks verification on three-dimensional problems

Impact

1. Academic value: Provides new insights for numerical solution of Richards equation
2. Practical value: Broad application prospects in hydrogeology and soil science
3. Reproducibility: Detailed implementation details facilitate reproduction

Applicable Scenarios

1. Porous media flow: Groundwater flow, soil water movement
2. Engineering applications: Dam seepage, slope stability analysis
3. Environmental problems: Contaminant migration, groundwater remediation
4. Material science: Fluid infiltration in fibrous materials

References

Key References

1. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1(5), 318-333.
2. Celia, M.A., et al. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26(7), 1483-1496.
3. van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892-898.
4. Manzini, G., Ferraris, S. (2004). Mass-conservative finite volume methods on 2-d unstructured grids for the Richards' equation. Adv. Water Resour. 27, 1199-1215.

Overall Assessment: This is a paper of significant innovative importance in the field of numerical solution of the Richards equation. By introducing a bounded auxiliary variable u, it cleverly resolves the numerical difficulties encountered by traditional methods when handling dry and saturated regions. The method has a solid theoretical foundation and practical value, with numerical experiments thoroughly validating its effectiveness. Although theoretical analysis requires further strengthening, it is overall a high-quality research work.