Coupled thermal-hydraulic model of the continental subsurface: theory, regimes, and reduced solver
Abstract
The thermal regime of the continental subsurface is shaped by the interplay of atmospheric forcing, groundwater flow, and solid Earth heat. While many existing models couple groundwater flow and heat transport, and in some cases incorporate freeze-thaw or thermo-hydro-mechanical processes, representing the full spectrum of two-way feedbacks often remains challenging. This paper has two components and keeps them distinct throughout. The first is theoretical: a continuum formulation that explicitly couples two-way thermal-hydraulic feedback with a smooth freeze-thaw treatment. From first principles, we derive the governing equations for mass and energy and perform a scaling analysis that yields nine dimensionless numbers, including a newly identified thermal-hydraulic coupling metric (ΓN α ). These numbers quantify the competition between pressure diffusion, saturation storage, thermal expansion, conduction, advection, and latent heat, and they define quantitative regime diagrams for the three interfaces (land surface, subsurface, crust-mantle). We reduce the full three-dimensional system to a closed one-dimensional vertical column suitable for basin-scale applications. The second component is a reduced thermal-hydraulic numerical implementation, narrower than the theory. This open-source finite-volume solver (implicit time stepping, upwind advection, apparent-heat-capacity freeze-thaw) carries the thermal-expansion feedback, Darcy advection and freeze-thaw, but omits the phase-change mass source, the poroelastic porosity evolution, the skeleton momentum balance, and the gas phase.
The solver reports its own mass-and energy-budget residuals, and those diagnostics show that the omitted source is a first-order wherever the freezing point is crossed. The solver is cross-verified against an independent finite-difference Crank-Nicolson code and passes six analytical benchmarks (Carslaw-Jaeger, Stefan, Theis, Terzaghi, undrained thermo-poroelastic ratio, Bonacina), most of which are operator verifications. Four canonical 1-D case studies (permafrost, geothermal reservoir, arid basin, thermo-poro coupled) demonstrate consistency between the designed inputs and the proposed theoretical classification, and a 2-D extension explores lateral basin recharge-discharge and permafrost edge effects. The open-source solver, within its applicability domain, equipped with a test suite and conservation diagnostics, provides a computational tool for attributing subsurface thermal anomalies to their physical drivers and for projecting future thermal evolution under changing environmental conditions.
Related
- Gliding through marine heatwaves: Subsurface biogeochemical characteristics on the Australian continental shelf
- Thermal and Hydrological Controls on Soil Respiration and Subsurface Gas Transport in Arctic Tundra Ecosystems.
- Rheology-controlled hydraulic selection in fracture-matrix heat transport: Mechanisms and thermal signatures
- Rheology-controlled hydraulic selection in fracture-matrix heat transport: Mechanisms and thermal signatures
- Geochemical evaluation of subsurface hydrogen-brine-mineral interactions and hydrogen loss
BibTeX
@article{rostami2026coupled,
title = {{Coupled thermal-hydraulic model of the continental subsurface: theory, regimes, and reduced solver}},
author = {Masoud Rostami and Bijan Fallah and Farahnaz Fazel-Rastgar and Li-Yun Fu},
journal = {Journal of Hydrology},
year = {2026},
month = sep,
doi = {10.1016/j.jhydrol.2026.136412},
url = {https://enpc.hal.science/hal-05748493v1},
note = {CC BY-NC-ND 4.0},
}
RIS
TY - JOUR
TI - Coupled thermal-hydraulic model of the continental subsurface: theory, regimes, and reduced solver
AU - Masoud Rostami
AU - Bijan Fallah
AU - Farahnaz Fazel-Rastgar
AU - Li-Yun Fu
PY - 2026
DA - 2026/09/14
JO - Journal of Hydrology
DO - 10.1016/j.jhydrol.2026.136412
UR - https://enpc.hal.science/hal-05748493v1
AB - The thermal regime of the continental subsurface is shaped by the interplay of atmospheric forcing, groundwater flow, and solid Earth heat. While many existing models couple groundwater flow and heat transport, and in some cases incorporate freeze-thaw or thermo-hydro-mechanical processes, representing the full spectrum of two-way feedbacks often remains challenging. This paper has two components and keeps them distinct throughout. The first is theoretical: a continuum formulation that explicitly couples two-way thermal-hydraulic feedback with a smooth freeze-thaw treatment. From first principles, we derive the governing equations for mass and energy and perform a scaling analysis that yields nine dimensionless numbers, including a newly identified thermal-hydraulic coupling metric (ΓN α ). These numbers quantify the competition between pressure diffusion, saturation storage, thermal expansion, conduction, advection, and latent heat, and they define quantitative regime diagrams for the three interfaces (land surface, subsurface, crust-mantle). We reduce the full three-dimensional system to a closed one-dimensional vertical column suitable for basin-scale applications. The second component is a reduced thermal-hydraulic numerical implementation, narrower than the theory. This open-source finite-volume solver (implicit time stepping, upwind advection, apparent-heat-capacity freeze-thaw) carries the thermal-expansion feedback, Darcy advection and freeze-thaw, but omits the phase-change mass source, the poroelastic porosity evolution, the skeleton momentum balance, and the gas phase. The solver reports its own mass-and energy-budget residuals, and those diagnostics show that the omitted source is a first-order wherever the freezing point is crossed. The solver is cross-verified against an independent finite-difference Crank-Nicolson code and passes six analytical benchmarks (Carslaw-Jaeger, Stefan, Theis, Terzaghi, undrained thermo-poroelastic ratio, Bonacina), most of which are operator verifications. Four canonical 1-D case studies (permafrost, geothermal reservoir, arid basin, thermo-poro coupled) demonstrate consistency between the designed inputs and the proposed theoretical classification, and a 2-D extension explores lateral basin recharge-discharge and permafrost edge effects. The open-source solver, within its applicability domain, equipped with a test suite and conservation diagnostics, provides a computational tool for attributing subsurface thermal anomalies to their physical drivers and for projecting future thermal evolution under changing environmental conditions.
SN - 0022-1694
SN - 1879-2707
ER -