A monolithic finite element formulation for coupled thermo-hydro-mechanical analysis with two-phase flow: Formulation and verification

A monolithic finite element formulation for fully coupled thermo-hydro-mechanical (THM) problems with two-phase flow in deformable porous media is presented, implemented in OpenGeoSys 5 with displacement, capillary pressure, non-wetting phase pressure, and temperature as primary variables and quadratic-linear mixed-order interpolation. The core methodological contribution is a modified Newton–Raphson solver in which the displacement-scalar coupling residual from the ɛ b S α ∂ ɛ v / ∂ t mass-balance terms is assembled in closed form and added to the scalar equation right-hand side at every iteration. This residual recovers Newton-like convergence from the Picard tangent matrix using only the existing displacement–pressure coupling block and the current displacement increment; no derivative of the constitutive laws is required, making it portable to frameworks with or without automatic differentiation. On the Kueper DNAPL benchmark (45,765 DOFs, MUMPS direct solver), this residual reduces the mean nonlinear iterations per step from 8.0 to 6.0 (25%, 645 fewer iterations across 320 steps) and the wall-clock time from 803 to 598 s (26%), relative to the same solver without the coupling term. Verification uses four benchmarks: (i) Booker & Savvidou point heat source for single-phase THM, in 2D axisymmetric and 3D geometries against the corrected analytical solution of Chaudhry et al. (2019); (ii) Kueper DNAPL infiltration, cross-verified against an in-framework segregated multiphase solver; (iii) a new thermal extension of the Kueper problem, verified against the analytical complementary error function thermal-diffusion solution; and (iv) the Liakopoulos drainage problem, verified against Grunwald et al. (2022). Mesh and time-step convergence on Liakopoulos, a pressure-smoothness diagnostic across the suite, and a lumped-versus-consistent mass comparison establish the spatial, temporal, and discrete-stability properties; the solver remains stable across all schedules tested, including a 100-fold time-step increase on the thermal Kueper benchmark.