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

Abstract: Accurate numerical simulation of fault and fracture mechanics is critical for the performance and safety assessment of many subsurface systems. The discretized representation of discontinuity surfaces and the robust simulation of their frictional contact behavior still represent major challenges. In this work, we use the mortar method to enforce the contact constraints and allow for non-conformity around the discontinuity surface, with a set of Lagrange multipliers playing the role of interface tractions. The formulation combines piecewise linear displacements with piecewise constant multipliers defined on one side of the fault interface (the non-mortar side). This choice for the Lagrange multipliers has a number of advantages from practical and computational viewpoints, but violates the inf-sup stability constraint. In order to stabilize the proposed formulation, we develop a traction-jump stabilization term to be added to the constraint equations. We use a macro-element analysis to derive an algorithmic strategy that automatically evaluates the proper scaling of the stabilization, without requiring any additional user-selected parameter. Numerical experiments demonstrate that the proposed formulation not only restores the inf-sup stability condition, but also recovers stable traction profiles in the presence of finer non-mortar sides, where other inf-sup-stable formulations fail. The proposed method is finally used to simulate non-linear contact conditions in non-conforming corner-point grids typically used in industrial geological applications.

• The paper introduces a stabilized mortar method using Lagrange multipliers for simulating frictional contact in geological media, addressing challenges with non-conforming grids.
• The method combines piecewise linear displacements and constant multipliers, which violates the inf-sup condition but is stabilized through a traction-jump term.
• The stabilization ensures smooth traction fields, as demonstrated in various test cases, making it applicable to industrial reservoir simulations

Motivation and scope

The paper addresses the numerical simulation of frictional contact on faults and fractures in subsurface geomechanics, a problem central to the safety assessment of geothermal systems, gas storage, and carbon sequestration. The authors adopt a Discrete Fracture Model (DFM) framework in which the discontinuity is explicitly discretized as a lower-dimensional surface, and enforce contact constraints via Lagrange multipliers within a single-pass mortar formulation. The distinguishing feature of the work is the combination of piecewise linear ($\mathbb{Q}_1$) displacements with piecewise constant ($\mathbb{P}_0$) multipliers defined on the non-mortar side of the interface — a pairing that violates the inf-sup condition but offers substantial practical advantages — together with a traction-jump stabilization term whose scaling is computed automatically from macro-element theory.

The choice of $\mathbb{P}_0$ multipliers is motivated by two well-known implementation difficulties in mortar methods applied to geological problems. First, Dirichlet boundary conditions intersecting the fault boundary can be imposed without local modification of multiplier basis functions [puso2003mesh]. Second, solvability issues at cross-points of intersecting interfaces are inherently resolved, which matters because geological models frequently contain intricate fault networks. The price paid is the loss of discrete inf-sup stability, which the paper restores through stabilization rather than through space enrichment.

Problem setting and discretization

The model consists of two elastic bodies $\Omega^{(1)}$ and $\Omega^{(2)}$ in $\mathbb{R}^3$, separated by a possibly non-conforming interface partitioned into a non-mortar side $\Gamma_f^{(1)}$ (where all contact quantities are evaluated) and a mortar side $\Gamma_f^{(2)}$. Under infinitesimal strains and quasi-static conditions, the strong form comprises the elastic equilibrium equations coupled with unilateral normal contact conditions and a Coulomb friction law with limiting tangential traction $\tau_{lim}(t_N) = c - t_N\tan(\phi)$. The interface jump operator involves an orthogonal projection $\Pi$ mapping functions from the mortar to the non-mortar side, so that mortar coupling matrices require integration of products of basis functions living on unrelated grids — identified by the authors as the major algorithmic challenge of the method.

The variational inequality is resolved with an active-set strategy iterating on the KKT conditions until the stick/slip/open partition stabilizes, combined with Newton linearization at each active-set iteration. The resulting Jacobian has saddle-point structure $\mathbb{P}_0$0, reducing to the classical form $\mathbb{P}_0$1 when the entire interface is in stick. Nonlinearity enters through the derivatives of the limiting tangential traction with respect to both the displacement increment and the normal traction.

Stability analysis and stabilization

For the unstable pair $\mathbb{P}_0$2, the Schur complement $\mathbb{P}_0$3 of the stick-state system possesses a nontrivial kernel whose eigenmodes correspond to spurious traction profiles. The proposed remedy adds a bilinear traction-jump term to the constraint equation,

$\mathbb{P}_0$4

summed over internal edges of the non-mortar mesh, in analogy with pressure-jump stabilization for Stokes. The stabilized system satisfies $\mathbb{P}_0$5 by construction for any positive definite scaling tensor; the difficulty lies in choosing $\mathbb{P}_0$6 large enough to remove spurious modes but small enough not to over-smooth the tractions and corrupt the constraint enforcement.

The paper's main algorithmic contribution is a parameter-free scaling rule derived from macro-element analysis. For each internal mortar node, a local macroelement is assembled comprising the adjacent mortar cells and the overlapping non-mortar cells, and a local Schur complement approximation $\mathbb{P}_0$7 is formed using only the diagonal of the local stiffness blocks. The scaling tensor on each internal edge is then the average of the face-related diagonal submatrices of $\mathbb{P}_0$8. This extends prior conforming-grid work to the non-conforming case and eliminates any user-selected stabilization coefficient — a notable practical advantage over penalty-type or ad hoc stabilized schemes.

Stability is verified numerically via the standard discrete inf-sup test: $\mathbb{P}_0$9 is obtained as the square root of the smallest eigenvalue of $\mathbb{P}_0$0 on successively refined patches. Two results are reported: under uniform refinement at fixed refinement ratio $\mathbb{P}_0$1, the unstabilized constant collapses to zero while the stabilized one remains bounded away from zero; refining only the non-mortar side also yields an essentially constant $\mathbb{P}_0$2. It should be noted that the test relies on a weaker reformulation of the inf-sup condition (avoiding evaluation of the fractional dual norm), so it provides a necessary rather than sufficient certificate; no rigorous a priori error analysis for the stabilized non-conforming formulation is given in this paper.

Numerical validation

Five test cases assess accuracy, robustness, and applicability:

• Patch test: on a unit cube split into eight differently discretized blocks meeting at a cross-point, the stabilized $\mathbb{P}_0$3 formulation reproduces the exact constant normal traction $\mathbb{P}_0$4, whereas standard nodal multipliers — despite being theoretically inf-sup stable — exhibit spurious oscillations along the non-conforming interfaces. This directly demonstrates that inf-sup stability alone does not guarantee clean traction fields in the presence of cross-points.
• Constant sliding: the analytical fully-sliding solution with $\mathbb{P}_0$5 is reproduced exactly, converging in a single active-set update after the first elastic step.
• Single fracture under compression: good agreement with the closed-form solutions for normal traction and sliding profile; the known pathological error at the fault tip, caused by the abrupt slip-to-stick transition, persists but is smoothed by the combined effect of grid non-conformity and stabilization.
• Two-block compression/shear test with $\mathbb{P}_0$6: this is the key discriminating experiment. When the non-mortar side is finer, the interface is over-constrained and all compared inf-sup-stable alternatives — standard nodal multipliers, dual multipliers [popp2010dual], and bubble-enriched formulations [hauret2007discontinuous] — produce oscillatory tractions, whereas the proposed traction-jump stabilization yields a smooth profile closely matching a fine conforming reference. Since oscillatory $\mathbb{P}_0$7 feeds directly into the Coulomb limit $\mathbb{P}_0$8, these oscillations would compromise active-set robustness; their elimination is therefore consequential for the nonlinear solver, not merely cosmetic. The test also illustrates how boundary-condition enforcement with nodal multipliers propagates oscillations along the interface due to their non-local support.
• Aquifer withdrawal in a faulted domain: a ten-year simulation of fluid withdrawal inducing fault reactivation, discretized with a non-conforming corner-point grid generated with MRST — the industrial-standard format whose misaligned vertices motivate the entire methodology. Sliding initiates near the top of the domain and propagates downward until nearly the entire fracture opens, with smoothly evolving contact states throughout. The resulting fault gap is of the same order as the maximum settlement near the wells, indicating mechanically significant reactivation.

Limitations and open questions

Several caveats accompany the results. The stability argument is empirical: the inf-sup verification uses a weaker necessary condition and is performed on a canonical patch, without a formal proof of uniform stability or convergence rates for the stabilized non-conforming three-dimensional formulation. The derivation assumes the entire interface in stick for the saddle-point analysis, and the extension to mixed stick/slip/open configurations rests on the active-set decomposition. The scaling strategy requires the projection operator $\mathbb{P}_0$9 to be sufficiently smooth, an assumption stated but not analyzed in detail. The numerical examples use hexahedral elements exclusively; extension to tetrahedra is claimed to be straightforward but is not demonstrated. Finally, the solver technology is limited to a standard active-set strategy, and the authors themselves identify more sophisticated contact algorithms and more complex geometrical scenarios — including intersecting fracture networks at scale — as open directions.

Conclusion

This paper contributes a parameter-free, macro-element-scaled traction-jump stabilization that renders the practically convenient $\Omega^{(1)}$0 single-pass mortar pairing viable for three-dimensional frictional contact on non-conforming grids. Its most substantive claim — supported by direct comparison against dual and bubble-stabilized multipliers — is that restoring nominal inf-sup stability does not suffice when the non-mortar side is finer, whereas the proposed stabilization recovers smooth, reference-matching traction profiles in exactly that regime. The successful application to a corner-point-grid aquifer compaction problem establishes relevance to industrial reservoir simulation workflows, where non-conforming fault representations are unavoidable.