Auto-Stabilized Weak Galerkin Finite Element Methods for Biot's consolidation model on Non-Convex Polytopal Meshes

Abstract: This paper presents an auto-stabilized weak Galerkin (WG) finite element method for the Biot's consolidation model within the classical displacement-pressure two-field formulation. Unlike traditional WG approaches, the proposed scheme achieves numerical stability without the requirement of traditional stabilizers. Spatial discretization is performed using weak Galerkin finite elements for both displacement and pressure approximations, while a backward Euler scheme is employed for temporal discretization to ensure a fully implicit and stable formulation. We establish the well-posedness of the resulting linear system at each time step and provide a rigorous error analysis, deriving optimal-order convergence. A significant merit of this WG scheme is its flexibility on general shape-regular polytopal meshes, including those with non-convex geometries. By utilizing bubble functions as a primary analytical tool, the method produces stable, oscillation-free pressure approximations without specialized treatment. Numerical experiments are presented to validate the theoretical convergence rates and demonstrate the computational efficiency and robustness of the auto-stabilized formulation.

• The paper develops an auto-stabilized weak Galerkin finite element method for the quasi-static Biot consolidation model, enhancing numerical stability without traditional stabilizer terms.
• The method achieves well-posedness, optimal-order convergence, and oscillation-free pressure approximations on general shape-regular polytopal meshes, including non-convex elements
• Experiments validate the theory, showing locking-free behavior, optimal convergence rates, and robustness against low-permeability scenarios

Overview and contribution

This paper develops a weak Galerkin (WG) finite element method for the quasi-static Biot consolidation model in its classical two-field displacement–pressure formulation, discretized in time by backward Euler. The distinguishing feature is that the scheme is auto-stabilized: numerical stability is achieved without the traditional WG stabilizer terms, which are instead replaced by computing discrete weak differential operators against polynomial spaces of elevated degree. Specifically, for an element $T$ with $N$ faces, the degree of the operator space is taken as $r = k-1+N$ on convex elements and $r = k-1+2N$ on non-convex elements, where $k$ is the degree of the solution space (2603.27704). This strategy preserves the global sparsity and size of the stiffness matrix while removing stabilizer-dependent implementation complexity.

The paper claims three principal results: (i) well-posedness of the fully discrete linear system at each time step via a Babuška inf-sup argument; (ii) optimal-order convergence under standard regularity assumptions; and (iii) oscillation-free pressure approximations on general shape-regular polytopal meshes—including non-convex elements—without mass lumping or tuned stabilization parameters. The analysis leans heavily on bubble functions as the primary tool for handling non-convex geometry, extending a program previously applied to Poisson, biharmonic, elasticity, Maxwell, and Stokes problems by the same authors.

Discrete formulation

The spatial discretization uses equal-order weak finite element spaces: displacement approximated by $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$ with $\mathbf{v}_0\in[P_k(T)]^d$ and boundary component $\mathbf{v}_b\in[P_k(e)]^d$, and pressure analogously with $q_0\in P_k(T)$, $q_b\in P_k(\partial T)$. The bilinear forms are built from the discrete weak strain tensor $N$0, weak divergence $N$1, and weak gradient $N$2, defined elementwise through integration-by-parts identities tested only against polynomials of degree $N$3. Because $N$4 exceeds what would be needed for consistency alone, the resulting operators automatically encode sufficient stability—the essence of auto-stabilization.

The fully discrete backward Euler scheme at time $N$5 reads: find $N$6 such that

$N$7

$N$8

with $N$9, $r = k-1+N$0, and $r = k-1+N$1. Notably, the displacement space inherits the locking-free property established for primal-formulation WG elasticity (2603.27704), so the pair remains locking-free as $r = k-1+N$2 grows.

Well-posedness

Well-posedness rests on two norm equivalences—between the energy norms $r = k-1+N$3, $r = k-1+N$4 and the discrete $r = k-1+N$5 semi-norms—which hold because the elevated-degree operators control the jump terms $r = k-1+N$6. Combined with Korn-type arguments (rigid-body motions killed by clamped boundary conditions on $r = k-1+N$7) and the pressure constraint on $r = k-1+N$8, these establish that the energy semi-norms are genuine norms on the constrained spaces. A discrete inf-sup condition,

$r = k-1+N$9

is invoked from prior work, and the coupled form $r = k-1+2N$0 is shown to satisfy an inf-sup bound with constant $r = k-1+2N$1 independent of both $r = k-1+2N$2 and $r = k-1+2N$3, using the test choice $r = k-1+2N$4, $r = k-1+2N$5. Consequently the linear system at each time step is uniquely solvable by Babuška theory. The independence of $r = k-1+2N$6 from $r = k-1+2N$7 is important: it means well-posedness does not degrade in the small-time-step regime where pressure oscillations typically arise.

Error analysis

The error analysis follows the elliptic-projection splitting technique of Thomée. With $r = k-1+2N$8 denoting the $r = k-1+2N$9 projection onto the piecewise-$k$0 weak spaces and $k$1 the projection onto the degree-$k$2 operator space, the key commutativity properties $k$3, $k$4, and $k$5 reduce consistency errors to boundary functionals $k$6 involving the jumps $k$7 and the projection residuals $k$8 applied to $k$9, $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$0, $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$1, and $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$2. These are bounded as

$\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$3

using trace inequalities and the approximation power of the degree-$\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$4 projection. The main theorem then gives, assuming $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$5, $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$6, and analogous regularity for $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$7:

$\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$8

Thus the method delivers first-order convergence in time (as expected from backward Euler) and optimal order in space. Two caveats should be noted plainly: the analysis assumes the initial data satisfies the discrete divergence-free condition $\mathbf{v}=\{\mathbf{v}_0,\mathbf{v}_b\}$9 (the authors state alternative initializations do not alter the analysis but do not prove this), and the temporal estimate involves $\mathbf{v}_0\in[P_k(T)]^d$0 rather than a stronger norm, consistent with but not sharper than standard parabolic theory.

Numerical validation

Two experiment sets support the theory. In the first, a manufactured solution $\mathbf{v}_0\in[P_k(T)]^d$1, $\mathbf{v}_0\in[P_k(T)]^d$2 is approximated on families of genuinely non-convex polygonal grids using $\mathbf{v}_0\in[P_k(T)]^d$3-$\mathbf{v}_0\in[P_k(T)]^d$4 ($\mathbf{v}_0\in[P_k(T)]^d$5) and $\mathbf{v}_0\in[P_k(T)]^d$6-$\mathbf{v}_0\in[P_k(T)]^d$7 elements, at Poisson ratios $\mathbf{v}_0\in[P_k(T)]^d$8 and $\mathbf{v}_0\in[P_k(T)]^d$9. Representative convergence rates include:

Element	Quantity	Observed order
$\mathbf{v}_b\in[P_k(e)]^d$0-$\mathbf{v}_b\in[P_k(e)]^d$1	$\mathbf{v}_b\in[P_k(e)]^d$2	2.0
$\mathbf{v}_b\in[P_k(e)]^d$3-$\mathbf{v}_b\in[P_k(e)]^d$4	$\mathbf{v}_b\in[P_k(e)]^d$5	2.0
$\mathbf{v}_b\in[P_k(e)]^d$6-$\mathbf{v}_b\in[P_k(e)]^d$7	$\mathbf{v}_b\in[P_k(e)]^d$8	up to 6.4 at coarse levels, settling to ~4

Optimal-order convergence is observed in every configuration, and—significantly—the rates are insensitive to $\mathbf{v}_b\in[P_k(e)]^d$9, confirming the locking-free claim empirically even near the incompressible limit. Some pre-asymptotic orders exceed the theoretical rate (e.g., 7.1 for the $q_0\in P_k(T)$0 displacement $q_0\in P_k(T)$1 error on one grid), which the authors report without explanation; these likely reflect superconvergence effects not analyzed here.

The second experiment computes the steady state of a layered-medium problem with hydraulic conductivity dropping to $q_0\in P_k(T)$2 in the interior strip—a low-permeability regime where conventional discretizations produce spurious pressure oscillations. The computed pressure exhibits a sharp internal layer at the inflow interface $q_0\in P_k(T)$3 with no visible oscillation or interface smearing, providing direct evidence for the central stability claim in exactly the regime (low permeability, effectively small effective time steps) where inf-sup violation manifests.

Limitations and open questions

Several limitations deserve explicit statement. First, the theoretical error analysis is carried out under full $q_0\in P_k(T)$4 regularity assumptions on the exact solution; no results are given for low-regularity or interface solutions, despite the second numerical example involving discontinuous coefficients. Second, the monotonicity/stability argument for oscillation-free pressures is empirical: the paper demonstrates absence of oscillations numerically but does not prove a discrete maximum principle or a formal monotonicity result, leaving open whether the auto-stabilized scheme is provably oscillation-free for arbitrary low-permeability contrasts and time-step sizes. Third, the degree elevation $q_0\in P_k(T)$5 or $q_0\in P_k(T)$6 increases local operator cost per element; the paper asserts preservation of global matrix size and sparsity but provides no operation counts or timing comparisons against stabilized alternatives. Finally, the extension to three-field formulations (with flux as an independent unknown) and to nonlinear poroelasticity is not addressed.

Conclusion

The paper establishes a stabilizer-free WG finite element method for Biot's consolidation model with rigorous well-posedness and optimal-order error estimates, valid on shape-regular polytopal meshes including non-convex elements. Numerical experiments confirm the predicted rates across polynomial degrees one through three, robustness with respect to the Poisson ratio up to $q_0\in P_k(T)$7, and oscillation-free pressure in a severe low-permeability setting. The main open issues are a formal proof of pressure monotonicity, error analysis under reduced regularity, and quantitative assessment of the computational overhead introduced by the elevated-degree discrete operators.