🔢 mathematics

A Stabilized Mortar Method for Discontinuities in Geological Media with Non-Conforming Grids

This paper presents a stabilized mortar method using a traction-jump term and an automated macro-element scaling strategy to enforce contact constraints on non-conforming grids, thereby restoring inf-sup stability and enabling robust simulation of frictional fault mechanics in geological media.

Original authors: Daniele Moretto, Andrea Franceschini, Massimiliano Ferronato

Published 2026-03-20

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Daniele Moretto, Andrea Franceschini, Massimiliano Ferronato

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Earth's crust as a giant, complex puzzle made of different types of rocks, sand, and clay. Sometimes, this puzzle has cracks or "fault lines" running through it, like a broken eggshell. When we want to predict how these rocks will move—perhaps because we are pumping water out of an underground aquifer or storing carbon dioxide deep underground—we need a computer simulation.

The problem is that these cracks are messy. They don't line up perfectly with the grid lines we use to draw our computer models. It's like trying to fit a square peg into a round hole, or trying to stitch together two pieces of fabric where the patterns on one side don't match the patterns on the other.

This paper presents a new, smarter way to stitch these mismatched pieces together so the computer simulation doesn't fall apart or give crazy, wrong answers.

The Problem: The "Mismatched Puzzle"

In the past, scientists tried to force the computer grid to match the crack perfectly. But in real-world geology, the rock layers are often jagged and irregular (using what's called a "corner-point grid"). Forcing a perfect match is like trying to force a square peg into a round hole; it either breaks the model or requires so much computing power it becomes impossible.

So, scientists use a method called the Mortar Method. Think of this like a "bridge" or a "mortar" that holds two mismatched bricks together. Instead of forcing the bricks to align perfectly, the mortar allows them to sit at slightly different angles while still holding them tight.

However, there's a catch. When you use this "bridge" method with certain types of math (specifically, using simple, blocky numbers to represent the forces), the simulation can become unstable. It's like building a house of cards: if the wind blows (or if the math gets tricky), the whole thing wobbles and collapses into nonsense. The computer starts predicting that the rocks are vibrating wildly or that the forces are jumping up and down in a jagged, unrealistic way.

The Solution: The "Shock Absorber"

The authors of this paper, Daniele Moretto, Andrea Franceschini, and Massimiliano Ferronato, came up with a clever fix. They realized that the instability happens because the "bridge" is too rigid in some spots and too loose in others.

They introduced a Stabilization Term, which acts like a shock absorber in a car.

Without the shock absorber: When the car hits a bump (a complex geological feature), the ride is bumpy and the passengers (the simulation results) get thrown around.
With the shock absorber: The system smooths out the bumps. The car still moves, but the ride is steady and safe.

In their math, this "shock absorber" is a special formula that looks at the difference in forces between the mismatched grid pieces and gently smooths them out. Crucially, they figured out how to make this shock absorber automatic. You don't have to guess how strong it should be; the computer calculates the perfect amount of "damping" based on the local geometry, just like a modern car's suspension adjusts to the road instantly.

Why This Matters

This new method is a game-changer for three main reasons:

It handles the messy reality: It works perfectly with the "corner-point grids" that oil and gas companies, and geothermal engineers, already use. They don't have to rebuild their entire modeling software to use this.
It stops the "jitter": Previous methods would sometimes show the rocks vibrating or the pressure jumping wildly, which makes the simulation useless for safety checks. This new method produces smooth, realistic results, even when the grid is very fine on one side and coarse on the other.
It handles complex boundaries: It deals easily with tricky spots where faults cross each other or hit the edge of the model, which used to cause the math to break down.

The Real-World Test

To prove it works, the team ran a simulation of an underground aquifer (a water-bearing rock layer) being drained over 10 years. As the water was pumped out, the ground above it started to sink, and the fault line began to slide and open up.

Using their new "stabilized mortar" method, the simulation showed the fault opening up smoothly and realistically, just as nature would behave. It showed that the ground would settle and the crack would widen in a predictable way, without any of the weird, jagged errors that plague older methods.

The Bottom Line

Think of this paper as inventing a new, super-strong, self-adjusting glue for geological models. It allows engineers to simulate how the Earth's crust behaves under stress—whether for building safe underground storage or predicting earthquake risks—without getting bogged down by the messy, mismatched geometry of real-world rocks. It makes the computer models more reliable, safer, and ready for the complex challenges of the real world.

1. Problem Statement

The accurate numerical simulation of fault and fracture mechanics is essential for assessing the safety and performance of subsurface systems (e.g., geothermal energy, carbon sequestration, gas storage). A primary challenge lies in discretizing discontinuity surfaces and robustly simulating their frictional contact behavior, particularly when using non-conforming grids.

Context: Industrial geological applications often utilize corner-point grids, which frequently result in non-matching meshes on either side of a fault interface.
The Conflict: Standard Mortar Methods allow for non-conforming grids by enforcing contact constraints via Lagrange multipliers defined on one side of the interface (the "non-mortar" side). However, using piecewise constant (P0) multipliers with piecewise linear (Q1) displacements violates the inf-sup (Ladyzhenskaya-Babuška-Brezzi) stability condition.
Specific Instability: While P0 multipliers offer practical advantages (simplifying boundary condition enforcement and handling cross-points), they produce spurious oscillatory traction profiles, especially when the non-mortar side is finer than the mortar side ( $r_h < 1$ ). Existing stable formulations (e.g., dual nodal multipliers or bubble enrichments) often fail to recover stable traction profiles in these specific refinement scenarios.

2. Methodology

The authors propose a stabilized mortar formulation that combines piecewise linear displacements with piecewise constant Lagrange multipliers, augmented by a novel stabilization term.

A. Mathematical Formulation

Governing Equations: The problem models quasi-static frictional contact using the Coulomb friction law (with cohesion $c$ and friction angle $\phi$ ). The contact conditions are partitioned into stick, slip, and open regions.
Discretization:
- Displacement Space ( $V_h$ ): Piecewise trilinear ( $Q1$ ) elements.
- Multiplier Space ( $M_h$ ): Piecewise constant ( $P0$ ) elements defined on the non-mortar side.
Weak Formulation: The contact constraints are enforced via Lagrange multipliers acting as interface tractions. The system is solved using an active-set strategy to handle the non-linear inequality constraints (stick/slip transitions).

B. The Stabilization Strategy

To restore stability without sacrificing the benefits of P0 multipliers, the authors introduce a traction-jump stabilization term added to the constraint equations.

Stabilization Term: $s(\mu, t_h) = \sum_{E_{KL}} \langle J\mu \rangle_{E_{KL}} S_{E_{KL}} J t_h \rangle_{E_{KL}}$ , where $J\cdot\rangle$ represents the jump across internal edges.
Algorithmic Scaling (Key Innovation): Instead of requiring user-defined stabilization parameters, the paper derives an algorithm to automatically evaluate the scaling tensor $S_{E_{KL}}$ $S_{E_{K L}}$ .
- Macro-element Analysis: The method utilizes a local macro-element approach (inspired by [40]) to compute a local approximation of the Schur complement.
- Automatic Tuning: The stabilization matrix is constructed such that its entries have the same order of magnitude as the Schur complement of the unstabilized system. This ensures the removal of spurious modes without excessively smoothing the primary variables or introducing artificial stiffness.

3. Key Contributions

Novel Stabilized Formulation: Introduction of a traction-jump stabilization term specifically designed for P0 Lagrange multipliers in non-conforming mortar methods.
Parameter-Free Algorithm: Development of an algorithmic strategy based on local Schur complement estimates to automatically scale the stabilization term, eliminating the need for user-selected coefficients.
Restoration of Stability: Demonstration that this approach satisfies the inf-sup condition even when the non-mortar side is finer than the mortar side ( $r_h < 1$ ), a scenario where other stable formulations often fail.
Robustness in Complex Geometries: The method naturally handles cross-points (intersections of multiple faults) and Dirichlet boundary conditions on fault interfaces without the local modifications typically required for nodal multipliers.
Industrial Applicability: Successful application to corner-point grids, the industry standard for reservoir simulation, proving the method's viability for real-world geological modeling.

4. Numerical Results

The authors validate the method through five numerical tests:

Test 1 (Patch Test): Verifies the transfer of constant traction across a non-conforming interface with cross-points. The stabilized P0 formulation reproduces the exact constant traction, whereas standard nodal multipliers exhibit spurious oscillations.
Test 2 (Constant Sliding): A 2D sliding contact problem where the numerical solution exactly matches the analytical solution with minimal active-set updates.
Test 3 (Single Fracture under Compression): Validates the model against analytical solutions for normal traction and tangential sliding. The stabilization helps smooth oscillations near the fracture tip caused by the slip-stick transition.
Test 4 (Two Blocks under Compression/Shear):
- Scenario: A vertical crack with a finer non-mortar mesh ( $r_h = 0.5$ ).
- Comparison: The proposed method yields a smooth traction profile. In contrast, standard nodal multipliers and bubble-stabilized formulations (which are theoretically inf-sup stable) exhibit significant oscillations due to over-constraint in the $r_h < 1$ regime.
Test 5 (Aquifer Withdrawal in Faulted Domain):
- Scenario: A realistic 3D simulation of fluid withdrawal from a faulted geological formation using a corner-point grid.
- Outcome: The simulation successfully captures the evolution of the fault from a stick state to sliding and opening. The traction profiles remain smooth, and the method handles the complex, non-aligned grid geometry typical of industrial reservoir simulators.

5. Significance and Conclusion

This work addresses a critical gap in geomechanical modeling: the ability to use non-conforming grids (essential for complex geological structures and industrial workflows) while maintaining numerical stability and accuracy in contact tractions.

Practical Impact: By enabling the use of P0 multipliers with automatic stabilization, the method simplifies the implementation of boundary conditions and handles complex fault intersections (cross-points) more robustly than previous methods.
Theoretical Contribution: It provides a rigorous, parameter-free mechanism to satisfy the inf-sup condition in mixed formulations where the multiplier space is "too rich" relative to the displacement space on the non-mortar side.
Future Outlook: The authors suggest extending the method to more complex geometries and integrating more sophisticated contact algorithms beyond the standard active-set strategy.

In summary, the proposed traction-jump stabilized mortar method offers a robust, accurate, and industrially relevant solution for simulating fault mechanics in subsurface geological media with non-conforming discretizations.