Analytical Geothermal Model Simulates Ground Freezing Within Seconds With High Accuracy

A novel analytical modeling framework can simulate the complex transient freezing and thawing of soil surrounding shallow geothermal heat collectors in seconds, delivering temperature predictions within 1.4 Kelvin of intensive numerical simulations. Developed by researchers Adinda Van de Ven, Peter Bayer, and Roland Koenigsdorff and published in the journal Geothermal Energy, the method resolves a longstanding computational bottleneck in shallow geothermal engineering. Ground-source heat pumps offer substantial environmental advantages over air-source units, yet widespread deployment requires highly reliable, computationally accessible design tools that accurately reflect seasonal subsurface thermodynamics. Traditional planning software frequently relies on simplified analytical solutions developed for deep borehole heat exchangers, which operate under relatively stable temperatures deep underground where soil freezing is absent. In contrast, shallow ground heat collectors operate close to the surface, where severe seasonal winter chill and intensive thermal extraction routinely drop local temperatures below the freezing point of pore water, dramatically altering surrounding ground conditions.

Shallow ground heat collectors are closed-loop heat exchangers buried within the upper few meters of the ground, placing them in direct thermal contact with fluctuating atmospheric boundaries. Unlike deep vertical boreholes, shallow collectors interact with a significantly smaller surrounding soil volume and therefore possess a lower inherent ground heat storage capacity. During the peak winter heating season, continuous extraction of energy causes the operating temperatures of these shallow systems to fall below the freezing threshold of water for extended operational intervals. Under these conditions, a substantial portion of the thermal energy drawn by the heat pump is derived directly from the latent heat released as groundwater undergoes phase change into ice. Because seasonal regeneration relies on natural energy inputs from solar radiation, ambient air heat transfer, and rainfall infiltration, shallow collectors are typically installed beneath unsealed ground surfaces to ensure complete annual thermal recovery before the next heating season commences.

Engineering guidelines impose strict operational boundaries on shallow collectors, requiring that the collector fluid outlet temperature remains at or above minus five degrees Celsius and that frozen soil zones around adjacent pipes or trench plates do not merge. Maintaining unfrozen pathways between collector elements is essential to allow meltwater and rainwater percolation, which prevents surface mud accumulation and protects infrastructure from frost heave uplift and ground settlement. However, incorporating subsurface phase change into analytical design tools has historically presented formidable mathematical challenges because freezing and thawing introduce strong nonlinearities that violate the governing differential equations of linear potential theory. Numerical simulators using finite difference or finite element techniques can track phase change through enthalpy or apparent heat capacity formulations, but their extensive computational demands make them impractical for iterative design tasks, automated optimization, or routine engineering planning where multiple layout variants must be tested rapidly.

To overcome these limitations, the research team constructed an analytical architecture that divides the simulation domain into the ground heat collector source system and the surrounding subsurface media. The subsurface conduction model employs fundamental solutions of potential theory, utilizing an infinite line source model with an isothermal ground surface boundary condition for horizontal pipe collectors and a finite plane source formulation for planar trench collectors. Natural seasonal variations in background ground temperature are superposed based on depth-dependent cyclic equations accounting for annual surface temperature amplitude and thermal diffusivity. To handle realistic, time-varying heat extraction schedules, the model applies a temporal superposition algorithm operating across discrete equidistant time increments. The historical contribution of all past thermal extraction steps is calculated using the physical properties of unfrozen ground, while the thermal response of the immediate time step is evaluated according to the current local phase state of the soil.

The core innovation of the new methodology lies in a dynamic thermal power balance applied at the collector wall whenever subsurface temperatures drop to the freezing point. The total operational heat extraction rate is partitioned into a conductive heat flow rate from the ground and a latent heat flow rate generated by moisture freezing or thawing. When temperatures remain above freezing, the latent heat flow rate remains zero and heat transfer proceeds purely conductively. Once the collector wall falls below freezing, the conductive heat flux is split into an unfrozen fraction governed by the temperature difference between the natural ground and the freezing point, and a frozen fraction driven by the temperature difference across the expanding ice layer. By balancing these fluxes against the heat extraction rate, the algorithm tracks the mass and volume of frozen soil accumulated or melted during each discrete interval across four distinct operational cases covering onset, ongoing freezing, partial thawing, and complete thaw.

Converting the accumulated quantity of frozen soil into precise collector wall temperatures requires geometry-specific shape factors that describe heat conduction across the evolving frost layer. For horizontal pipe collectors, the authors established that freezing spreads with near radial symmetry under typical frost-free burial depths, allowing the radial thickness of the frozen envelope to be computed directly from pipe geometry and frozen area. The resulting temperature difference between the collector wall and the outer freezing isotherm is then calculated using the steady-state heat conduction equation for concentric circular cylinders. For planar trench collectors, freezing expands in three dimensions, prompting the researchers to establish an equivalent frozen volume that lengthens both the collector dimensions and the perpendicular frost front. The shape factor for trench walls is derived from planar conduction principles, while the maximum horizontal frost extent is extracted through linear interpolation of precalculated steady-state numerical reference curves.

A vital operational discovery highlighted by the researchers is the existence of strict lower limits for computational time step durations in the analytical formulation. Because the analytical shape factor relies on a steady-state approximation across the frozen shell, applying excessively short time increments destabilizes the calculation routine. An undersized step causes the model to overestimate initial ice growth, triggering an overcompensating thaw in the subsequent step that initiates an escalating cycle of numerical oscillations. Through rigorous stability analyses, the authors determined that a minimum time step of four hours is necessary to ensure stable convergence for horizontal pipe collectors, whereas planar trench collectors require a minimum step size of forty-eight hours. These discrete time intervals provide robust computational stability without degrading simulation accuracy, ensuring that thermal inertia and latent energy balances remain physically consistent throughout multi-month seasonal simulations.

To validate the analytical approach, the investigators developed comprehensive high-resolution two-dimensional and three-dimensional numerical benchmarks within COMSOL Multiphysics using the Heat Transfer in Porous Media module. Within the fully saturated numerical domain, pore water was assigned as a phase change material governed by the apparent heat capacity method operating over a two-Kelvin phase transition interval. Validation was executed across multiple twelve-month operating scenarios. In the first scenario, featuring constant heat extraction alongside a constant surface temperature, the analytical model for horizontal pipes exhibited remarkable agreement with numerical benchmarks, demonstrating a maximum collector wall temperature deviation of only 0.14 Kelvin and a maximum horizontal frost boundary discrepancy of one millimeter, representing a ten percent deviation. Under the second scenario, which incorporated annual surface temperature oscillations and a six-month heat extraction cycle followed by natural summer ground regeneration, wall temperature discrepancies remained below 0.60 Kelvin with frost extent errors limited to five millimeters.

The third and most realistic validation scenario subjected the horizontal pipe configuration to an actual Northern Hemisphere seasonal heating profile extending from October through April alongside fluctuating surface temperatures. Under these realistic conditions, the analytical model maintained excellent performance, exhibiting a maximum temperature deviation of 1.33 Kelvin and a maximum frost extent difference of 1.73 centimeters, or eleven percent. Crucially, simulations conducted without accounting for soil phase change produced wall temperatures more than two Kelvin lower than actual conditions, demonstrating that ignoring latent heat leads engineers to severely oversize ground collector installations. In terms of computational efficiency, the analytical pipe model solved the annual cycle in less than fourteen seconds on a standard desktop computer equipped with an Intel Core i7 processor, compared to one hour and two minutes required by the numerical model. For planar trench collectors, the analytical code solved annual performance in just 0.14 seconds, whereas the numerical model required over 13.5 hours on a high-performance computing cluster.

For planar trench collectors, the analytical framework predicted minimal collector wall temperatures within 0.7 Kelvin of numerical results, although the maximum horizontal frost extent exhibited a greater deviation of approximately ten centimeters, or 14.4 percent. The authors attributed this larger spatial discrepancy to pronounced temperature variations across the planar trench surface in numerical simulations, where local surface temperatures varied by more than six Kelvin and created non-uniform, irregular frost boundaries that deviate from idealized plane geometry. Nevertheless, because the analytical model slightly overpredicts horizontal ice extent while underestimating minimal wall temperatures, its predictions provide a conservative safety margin that inherently guards against underdimensioning collector arrays. While the current formulation focuses on isolated collectors and does not yet account for thermal interference between neighboring pipes or plates, it provides an exceptionally fast, highly accurate foundation that empowers engineers to optimize shallow geothermal heating installations with confidence.

Subject of Research: Analytical modeling of soil freezing and thawing dynamics around shallow geothermal collectors

Article Title: Analytical modelling approach of ground freezing and thawing around ground heat collectors

Article References: Van de Ven, A., Bayer, P., & Koenigsdorff, R. (2025). Analytical modelling approach of ground freezing and thawing around ground heat collectors. Geothermal Energy, 13(1), Article 43. https://doi.org/10.1186/s40517-025-00367-0

Image Credits: AI Generated

DOI: 10.1186/s40517-025-00367-0

Keywords: ground heat collectors, shallow geothermal energy, analytical modeling, soil freezing and thawing, latent heat transfer, infinite line source, finite plane source, planar trench collectors, horizontal pipe collectors, heat pump systems, thermal power balancing, temporal superposition, ground thermal regimes, frost isotherm propagation, COMSOL Multiphysics, apparent heat capacity method, geothermal design tools