groundHeatRecovery

Open-source coupled finite-element model for ground heat recovery with ice formation. We solve the Richards equation (variably saturated groundwater flow) coupled to an enthalpy-based heat balance with van Genuchten--Clapeyron freezing, for soils containing embedded pipe heat exchangers. The methodological signature is a conforming 1D--3D pipe coupling whose discrete trace selects soil unknowns exactly. Applications include ground-source heat pump sizing, seasonal latent-heat storage, and permafrost/frost dynamics.

The model is built on NGSolve. It conserves mass (machine precision in 1D via a dual-mixed formulation, a small bounded error in primal 2D/3D) and energy, and is solved by a robust monolithic Newton method across the freezing front in 1D, 2D, and 3D.

Associated publication

This code accompanies:

C. Tello Fachin, F. Schranz, A. Witzig. Coupled Finite Element Modeling of Ground Heat Recovery with Ice Formation. BauSIM 2026.

If you use this code, please cite the paper. The software itself is archived on Zenodo (concept DOI 10.5281/zenodo.21381614, always resolving to the latest release); see CITATION.cff.

Installation

We use mise and uv for reproducible environments (Python locked to 3.12):

mise install
uv sync

Without mise/uv:

pip install -e .

NGSolve is pulled from PyPI as a dependency; see the NGSolve documentation for platform notes. Optional extras: dev (adds pytest for the test suite).

Package layout

• src/local_ngsolve_utils/ (installed editable) is the core library:
  • materials.py soil parameter presets; common.py CoefficientFunction utilities.
  • constitutive_law_{sharp,smooth,test}.py phase-change constitutive laws (smooth is the production-recommended, C-infinity-regularised variant).
  • adaptive_timestepping.py LTE-based adaptive BDF1 controller.
  • point_monitoring.py transient point probes with CSV export.
  • mesh_pipe_in_cube.py, mesh_horizontal_helix_in_box.py conforming pure-netgen.occ meshes with a Glue()d 1D pipe/helix embedded along soil-tetrahedron edges.
• debug/ runs the numerical experiments (the solvers that produce the paper figures):
  • run_1d_seasonal.py 1D seasonal column with a phase-change twin.
  • run_2d_freezing.py 2D freezing front across two soils.
  • run_3d_helix_freezing.py 3D embedded horizontal helix; export_3d_helix_vtu*.py export the solution history to VTU for rendering.
• benchmarks/ performance scaling study (bench_runner.py, run_study.py).

A curated subset of the development notebooks is planned for a future release.

Reproducing the paper experiments

Each of the three numerical examples in the paper is produced by a runner in debug/:

Paper section	Experiment	Runner
3.1	2D freezing front across two soils	debug/run_2d_freezing.py
3.2	3D embedded horizontal helix	debug/run_3d_helix_freezing.py
3.3	1D seasonal column (phase-change twin)	debug/run_1d_seasonal.py

# 3.3 - 1D seasonal column; --smoke = 150 days, no spin-up (the quick variant)
uv run python debug/run_1d_seasonal.py --smoke
uv run python debug/run_1d_seasonal.py --no_pc --smoke     # phase-change-off twin

# 3.1 - 2D freezing front; positional args: t_end[h] h_conv dt_max_ice soil h_top label
uv run python debug/run_2d_freezing.py 4.0 500 3 sandyloam -2.0 sandyloam

# 3.2 - 3D embedded helix (all flags optional; defaults keep the paper geometry)
NGS_THREADS=10 uv run python debug/run_3d_helix_freezing.py --t_end_h 0.25 --label sandyloam

Each runner writes its .npz (and, for 3D, a .vol.gz mesh) artifacts to debug/output/ by default. Set GHR_OUTDIR to redirect (thread count via NGS_THREADS, default 4). For the 3D case, debug/export_3d_helix_vtu.py and export_3d_helix_vtu_full.py read those artifacts back and export VTU for rendering.

Runtimes scale with the simulated duration and mesh: the --smoke 1D column is on the order of ten to twenty minutes on four threads, whereas the full multi-year columns and longer 3D transients used for the paper figures take hours. debug/run_1d_decade_batch.sh orchestrates the full 1D decade batch.

Tests

uv sync --extra dev
uv run pytest tests/ -q

The suite is an import-and-assemble smoke check (build a coarse conforming mesh, evaluate the constitutive law); it runs on every push via GitHub Actions.

License

GNU Lesser General Public License v2.1 (LGPL-2.1). See LICENSE. This allows use in proprietary, academic, and commercial software without imposing copyleft on applications.