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 (). 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.
