A hybrid analytical-PINN model for subsurface simulation of geothermal heat exchangers in heterogeneous underground

# A hybrid analytical-PINN model for subsurface simulation of geothermal heat exchangers in heterogeneous underground
Source: [https://arxiv.org/html/2607.12271](https://arxiv.org/html/2607.12271)
Moke Rao, Thomas Hamacher, Smajil Halilovic Technical University of Munich, Chair of Renewable and Sustainable Energy Systems, Germany

Abstract

In this paper, a parametric physics\-informed neural network for solving the heterogeneous soil thermal problem with borehole heat exchangers \(BHEs\) as singular sources is developed\. There are three novel features in the present framework; namely, \(i\) the singularity is naturally removed by using analytical line source models; \(ii\) using the explicit formulation for gradient thermal conductivity enables physics\-informed learning of the parametrization featuring the conductivity; \(iii\) the learned correction is utilized as an efficient universal corrector via superposition principles\. We first introduce the decomposition of the temperature change and transform the approximation of the entire heterogeneous response to the correction compensating the difference between the practical solution and idealized homogeneous approximation\. In such a way, the delta function singularity is excluded and the bulk heat transfer is captured for the sake of facilitating the effective training of the neural network\. The original problem is then reformulated as a governing correction diffusion or advection\-diffusion equation subject to a homogeneous initial condition\. The linearly varying thermal conductivity is used to model the soil heterogeneity\. We propose a physics\-informed neural network to approximate a universal corrector with respect to a single borehole with unit heat extraction rate\. As a result, the network is trained by minimizing the physics\-informed and data\-anchored loss function that is evaluated for sampled conductivity parameters on adaptively selected training points\. In addition, we include the location indicator function regarding the source as a feature input of network and find that it helps the network to process the local information\. We perform numerical tests to exhibit the effectiveness of the proposed method based on the infinite line source \(ILS\), finite line source \(FLS\) and moving finite line source \(MFLS\) models\.

Keywords:Shallow geothermal energy, Borehole heat exchanger, Physics\-informed neural networks, Analytical correction, Heterogeneous media

## 1Introduction

Decarbonizing heating and cooling is a central component of the transition toward low\-carbon energy systems\[[7](https://arxiv.org/html/2607.12271#bib.bib26),[14](https://arxiv.org/html/2607.12271#bib.bib29)\]\. Shallow geothermal energy, commonly utilized through ground\-source heat\-pump systems, provides a locally available and seasonally stable thermal resource for building thermal conditioning and underground thermal\-energy storage\[[28](https://arxiv.org/html/2607.12271#bib.bib28),[13](https://arxiv.org/html/2607.12271#bib.bib32),[27](https://arxiv.org/html/2607.12271#bib.bib27)\]\. Among the available ground heat exchanger configurations, vertical borehole heat exchangers \(BHE\) are particularly attractive\. A BHE circulates a heat\-carrier fluid through closed pipes installed in a vertical borehole and exchanges heat with the surrounding ground without directly extracting groundwater\[[26](https://arxiv.org/html/2607.12271#bib.bib25)\]\. Its long\-term performance depends strongly on the transient thermal response of the subsurface and on the thermal interaction among neighboring boreholes\[[3](https://arxiv.org/html/2607.12271#bib.bib24),[16](https://arxiv.org/html/2607.12271#bib.bib35)\]\.

Accurate and repeated forward simulation of the ground temperature field is therefore essential for the design and operation of BHE systems\[[31](https://arxiv.org/html/2607.12271#bib.bib17),[4](https://arxiv.org/html/2607.12271#bib.bib18)\]\. A broad range of analytical and semi\-analytical models have been developed for BHE simulation, including infinite line source, finite line source and moving finite line source models\[[5](https://arxiv.org/html/2607.12271#bib.bib3),[30](https://arxiv.org/html/2607.12271#bib.bib6),[24](https://arxiv.org/html/2607.12271#bib.bib7)\]\. They represent a borehole as a line source and provide physically interpretable descriptions of the radial transient response\. However, their computation efficiency is obtained at the expense of restrictive assumptions, such as idealized homogeneous and infinite ground domain, simplified boundary conditions and spatially uniform groundwater velocity\[[12](https://arxiv.org/html/2607.12271#bib.bib31)\]\. On the other hand, numerical methods provide greater flexibility and can resolve more complex soil properties, geometries and subsurface transport conditions\[[1](https://arxiv.org/html/2607.12271#bib.bib22),[23](https://arxiv.org/html/2607.12271#bib.bib23),[10](https://arxiv.org/html/2607.12271#bib.bib33)\]\. Over the last decades, a variety of specialized numerical solvers for BHE simulation have been developed by dimensional reduction and hybrid formulation, including TRNSYS\[[17](https://arxiv.org/html/2607.12271#bib.bib19)\], OpenGeoSys\[[29](https://arxiv.org/html/2607.12271#bib.bib21)\]and FEFLOW\[[6](https://arxiv.org/html/2607.12271#bib.bib20)\]\. However, long\-term three\-dimensional numerical simulations remain computationally prohibitive, especially for repeated evaluations under varying configurations and operating conditions\[[11](https://arxiv.org/html/2607.12271#bib.bib30)\]\.

Applying deep neural networks to repeated BHE and subsurface thermal simulation has thus drawn much attention recently\. A trained neural surrogate can predict temperature fields at a substantially lower online cost compared to numerical simulation and especially benefits any procedure requiring numerous simulation runs\. In terms of training mechanism, the data\-driven and physics\-informed neural networks define the loss based on different sources of information\. The quality and quantity of training samples largely determine the performance of the purely data\-driven models\. In contrast, physics\-informed neural networks \(PINNs\) augment the learning process by introducing governing partial differential equations \(PDEs\), initial conditions and boundary conditions into the training objective\[[25](https://arxiv.org/html/2607.12271#bib.bib16)\]\. Consequently, the use of physical constraints can reduce the expensive data generation burden, allow differentiable approximations over the space\-time domain and enhance the generalization ability outside the sampled parameter range\. These properties make PINNs attractive as reusable thermal simulators for BHE systems in heterogeneous media\. However, their utilization in the field is still lacking due to the challenging nature of the singular \(BHE\) sources\[[15](https://arxiv.org/html/2607.12271#bib.bib34)\]\.

Solving PDEs with singular sources is challenging for standard DNNs and vanilla PINNs since they lack the expressive power to effectively approximate the solution containing a highly localized singularity\. Therefore, several developments for physics\-informed neural solvers have been proposed to deal with the nonsmoothness of PDE solution\. Gao et al\.\[[9](https://arxiv.org/html/2607.12271#bib.bib9)\]proposed a failure\-informed adaptive sampling strategy that places additional collocation points in regions with large residual\-based failure probabilities\. This is further combined with re\-sampling and subset simulation to improve the detection of localized high\-error regions in\[[8](https://arxiv.org/html/2607.12271#bib.bib10)\]\. Adaptive sampling is also adapted in\[[21](https://arxiv.org/html/2607.12271#bib.bib4)\]where the localized source is represented by a normalized Gaussian profile\. Whereas in terms of approximating point sources, Huang et al\. use a continuous probability density and combine this regularization with an uncertainty\-based loss weighting strategy and a practical multiscale neural architecture in\[[20](https://arxiv.org/html/2607.12271#bib.bib12)\]\. Making use of the separation of the singularity from the network approximation target,Huet al\.\[[18](https://arxiv.org/html/2607.12271#bib.bib11)\]proposed a singularity\-splitting deep Ritz method that represents the known singular component analytically and trains the neural network only for the regular remainder\. Furthermore, Hu et al\.\[[19](https://arxiv.org/html/2607.12271#bib.bib13)\]developed a singularity\-enriched PINN that explicitly incorporates explicit singular basic functions for Poisson problems to construct a more expressive approximation space of PINN\. Later this framework was extended to solving inverse problems by jointly identifying the locations and intensities of unknown point sources\.Laiet al\.\[[22](https://arxiv.org/html/2607.12271#bib.bib15)\]handle elliptic problems with singular sources on an interface by formulating the variational problem and use level set function to give the network additional information about the interface\.

In this paper, adopting the idea that localized singularities can be extracted by analytical solution, we propose a parametric hybrid PINN method for solving the soil temperature simulation problem with singular line sources \(BHEs\)\. The contributions of the proposed method are three\-fold\. First, we exclude the delta function singularity embedded in the source terms by using analytical solutions as base models\. Second, we target the linearly changing thermal conductivity cases, which allows training the network with physics\-informed loss and generalization to unseen gradient conductivity scenarios\. Third, the parametric neural network is trained on various sampled conductivity conditions to approximate the temperature correction induced by a singular borehole at the origin with unit heat transfer rate, which could be employed as a computationally efficient universal corrector by applying superposition principle\. This helps avoid training different networks for diverse arrangements of BHEs and source load profiles\. In addition, we rewrite the governing equation into a formulation on relative coordinates with respect to the borehole location, and then the network learns the output only considering the distance from the source under a certain thermal conductivity in an open\-boundary domain\. We thus include the distance from the borehole as input feature of the network that effectively improves the simulation accuracy near the source\.

The rest of the paper is organized as follows\. Section[2](https://arxiv.org/html/2607.12271#S2)introduces the classical formulations for the global soil temperature model problem\. In Section[3](https://arxiv.org/html/2607.12271#S3), the parametric hybrid\-PINN universal corrector method to simulate the models is presented\. In Section[4](https://arxiv.org/html/2607.12271#S4), we test and compare several numerical examples concerning ILS, FLS and MFLS analytical models, followed by conclusions of the present method and future work directions in Section[5](https://arxiv.org/html/2607.12271#S5)\.

## 2Problem Formulation

Looking into the transient propagation of temperature perturbations in the subsurface induced by a field of BHEs, we consider the global thermal problem as the evolution of the soil temperature field around the BHEs\. AssumeΩs\\Omega\_\{s\}denotes the soil domain andI=\(0,T\)I=\\left\(0,T\\right\)is the time interval of interest\. The soil temperatureT​\(𝐱,t\)T\\left\(\\mathbf\{x\},t\\right\),𝐱∈Ωs\\mathbf\{x\}\\in\\Omega\_\{s\},t∈It\\in Iin a heterogeneous porous medium is governed by the combined effects of thermal diffusion and groundwater advection\. We first consider the conduction\-dominated case

ρs​Cs​∂T∂t−∇⋅\(λ​\(𝐱\)​∇T\)=f​\(𝐱,t\),\\displaystyle\\rho\_\{s\}C\_\{s\}\\frac\{\\partial T\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda\(\\mathbf\{x\}\)\\nabla T\\right\)=f\\left\(\\mathbf\{x\},t\\right\),\(1\)whereρs​Cs\\rho\_\{s\}C\_\{s\}is the volumetric heat capacity of the soil andλ​\(𝐱\)\\lambda\(\\mathbf\{x\}\)denotes the thermal conductivity of the soil\. Hence the temperature changeT~​\(𝐱,t\)=T​\(𝐱,t\)−Tu\\widetilde\{T\}\\left\(\\mathbf\{x\},t\\right\)=T\\left\(\\mathbf\{x\},t\\right\)\-T\_\{u\}relative to the undisturbed constant ground stateTuT\_\{u\}follows the same mechanism

ρs​Cs​∂T~∂t−∇⋅\(λ​\(𝐱\)​∇T~\)=f​\(𝐱,t\)\.\\displaystyle\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda\(\\mathbf\{x\}\)\\nabla\\widetilde\{T\}\\right\)=f\\left\(\\mathbf\{x\},t\\right\)\.\(2\)Taking the groundwater flow effect into account, BHE\-induced temperature change satisfies the advection\-diffusion equation

ρ​C​∂T~∂t\+ρf​Cf​𝐯⋅∇T~−∇⋅\(λ​\(𝐱\)​∇T~\)=f​\(𝐱,t\),\\rho C\\frac\{\\partial\\widetilde\{T\}\}\{\\partial t\}\+\\rho\_\{f\}C\_\{f\}\\mathbf\{v\}\\cdot\\nabla\\widetilde\{T\}\-\\nabla\\cdot\\left\(\\lambda\(\\mathbf\{x\}\)\\nabla\\widetilde\{T\}\\right\)=f\(\\mathbf\{x\},t\),\(3\)whereρ​C\\rho Cis the combined volumetric heat capacity of the solid and liquid phase, i\.e\.,ρ​C=ϵ​ρf​Cf\+\(1−ϵ\)​ρs​Cs\\rho C=\\epsilon\\rho\_\{f\}C\_\{f\}\+\\left\(1\-\\epsilon\\right\)\\rho\_\{s\}C\_\{s\},ρf​cf\\rho\_\{f\}c\_\{f\}is the volumetric heat capacity of groundwater and𝐯\\mathbf\{v\}is the Darcy velocity of the groundwater flow\.

Furthermore,f​\(𝐱,t\)f\\left\(\\mathbf\{x\},t\\right\)represents the heat source or sink induced by BHE operation\. For a field ofKKvertical boreholes, the heat source term can be written as a superposition of individual borehole sources,

f​\(𝐱,t\)=∑k=1Kfk​\(𝐱,t\),f\(\\mathbf\{x\},t\)=\\sum\_\{k=1\}^\{K\}f\_\{k\}\(\\mathbf\{x\},t\),wherefkf\_\{k\}denotes the source term associated with thekk\-th borehole\. Let𝐱k=\(xB​H​Ek,yB​H​Ek\)\\mathbf\{x\}^\{k\}=\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\)be the location of thekk\-th borehole in the horizontal plane andqk​\(t\)q\_\{k\}\(t\)be the heat injection or extraction rate per unit length\. The infinite line source idealizes the borehole as an infinitely long vertical source and writes

fkILS​\(𝐱,t\)=qk​\(t\)​δ​\(x−xB​H​Ek\)​δ​\(y−yB​H​Ek\),f\_\{k\}^\{\\mathrm\{ILS\}\}\(\\mathbf\{x\},t\)=q\_\{k\}\(t\)\\,\\delta\(x\-x\_\{BHE\}^\{k\}\)\\delta\(y\-y\_\{BHE\}^\{k\}\),\(4\)whereδ​\(⋅\)\\delta\(\\cdot\)is the Dirac delta distribution\.

A more realistic representation accounts for the finite borehole lengthLL\. If the borehole extends from the depthz=Dz=Dto the depthz=D\+Lz=D\+L, the finite line source term reads

fkFLS​\(𝐱,t\)=qk​\(t\)​δ​\(x−xB​H​Ek\)​δ​\(y−yB​H​Ek\)​χ\[D,D\+L\]​\(z\),f\_\{k\}^\{\\mathrm\{FLS\}\}\(\\mathbf\{x\},t\)=q\_\{k\}\(t\)\\,\\delta\(x\-x\_\{BHE\}^\{k\}\)\\delta\(y\-y\_\{BHE\}^\{k\}\)\\chi\_\{\[D,D\+L\]\}\(z\),\(5\)whereχ\[0,L\]​\(z\)\\chi\_\{\[0,L\]\}\(z\)is the indicator function along the borehole depth\.

The temperature field is initialized from the undisturbed state

T​\(𝐱,0\)=Tu,𝐱∈Ωs,T\(\\mathbf\{x\},0\)=T\_\{u\},\\quad\\mathbf\{x\}\\in\\Omega\_\{s\},and here we assume a constant temperature Dirichlet boundary condition is imposed on∂Ωs\\partial\\Omega\_\{s\}for \([1](https://arxiv.org/html/2607.12271#S2.E1)\),

T​\(𝐱,t\)=Tu,𝐱∈∂Ωs,t∈I\.T\(\\mathbf\{x\},t\)=T\_\{u\},\\quad\\mathbf\{x\}\\in\\partial\\Omega\_\{s\},\\ t\\in I\.Equivalently, for the temperature changeT~\\widetilde\{T\}, these conditions are enforced as

T~​\(𝐱,0\)=0,𝐱∈Ωs,\\widetilde\{T\}\(\\mathbf\{x\},0\)=0,\\quad\\mathbf\{x\}\\in\\Omega\_\{s\},\(6\)and for \([2](https://arxiv.org/html/2607.12271#S2.E2)\)

T~​\(𝐱,t\)=0,𝐱∈∂Ωs,t∈I\.\\widetilde\{T\}\(\\mathbf\{x\},t\)=0,\\quad\\mathbf\{x\}\\in\\partial\\Omega\_\{s\},\\ t\\in I\.\(7\)
In this work, we focus on a heterogeneous setting in which the soil thermal conductivity is modeled as a linear gradient field in the horizontal plane

λ​\(𝐱\)=λc\+λx​x\+λy​y,\\lambda\(\\mathbf\{x\}\)=\\lambda\_\{c\}\+\\lambda\_\{x\}x\+\\lambda\_\{y\}y,whereλc\\lambda\_\{c\}denotes the conductivity value at the horizontal coordinate origin\(x,y\)=\(0,0\)\\left\(x,y\\right\)=\\left\(0,0\\right\)and\(λx,λy\)\\left\(\\lambda\_\{x\},\\lambda\_\{y\}\\right\)is the constant gradient of spatially varyingλ​\(𝐱\)\\lambda\\left\(\\mathbf\{x\}\\right\)\. The conductivity is assumed to be invariant along the vertical direction\. Fully resolved subsurface conductivity fields are rarely available in practical shallow geothermal applications, while standard thermal response tests usually provide only effective or averaged thermal properties\. Therefore, the gradient conductivity field is adopted as a controlled representation of commonly encountered BHE system\-scale subsurface heterogeneity\. Compared with the homogeneous assumption, it captures more realistic subsurface heat transfer characteristics while remaining sufficiently simple for systematic analysis\.

## 3Hybird\-PINN Universal Corrector Method

### 3\.1Hybrid PINN Corrector Based on Analytical Models

We propose a hybrid PINN corrector method for problem \([2](https://arxiv.org/html/2607.12271#S2.E2)\) subject to conditions \([6](https://arxiv.org/html/2607.12271#S2.E6)\), \([7](https://arxiv.org/html/2607.12271#S2.E7)\) and for problem \([3](https://arxiv.org/html/2607.12271#S2.E3)\) with initial condition \([6](https://arxiv.org/html/2607.12271#S2.E6)\) by applying analytical models with respect to source terms \([4](https://arxiv.org/html/2607.12271#S2.E4)\) and \([5](https://arxiv.org/html/2607.12271#S2.E5)\)\. In this section, we discuss the case with a single source located inΩs\\Omega\_\{s\}with a unity load rate\. The case of multiple boreholes can be handled via linear superposition\.

Instead of approximating the entire BHE\-induced temperature field, we decompose the temperature change into an analytical base component and a learned correction component

T~​\(𝐱,t\)=T~m​\(𝐱,t\)\+T~cm​\(𝐱,t\),m∈\{ILS,FLS,MFLS\}\.\\widetilde\{T\}\(\\mathbf\{x\},t\)=\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t\)\+\\widetilde\{T\}\_\{c\}^\{m\}\(\\mathbf\{x\},t\),\\qquad m\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}\.Here,T~m\\widetilde\{T\}^\{m\}denotes the analytical prediction given by the selected line source model, andT~Cm\\widetilde\{T\}\_\{C\}^\{m\}is the correction field to be approximated by a neural network\. By incorporating analytical knowledge of the underlying solution to the problem of interest, the singularity due to the line source could be captured explicitly\. The correction component is then introduced to represent the discrepancy between the analytical reference approximation and global heterogeneous solution\. We consider three analytical base models depending on whether the groundwater flow and finite borehole length are taken into account\.

#### Infinite line source \(ILS\) model\.

For a two\-dimensional infinite homogeneous domain, ILS model assumes heat conduction without groundwater advection in the underground and expresses the radial temperature distribution caused by a single BHE subjected to a piecewise constant heat load overmmtime steps as

T~ILS​\(x,y,t\)=∑l=0m−1ql−ql−14​π​λ0​E1​\(r24​α0​\(t−tl\)\),\\widetilde\{T\}^\{\\mathrm\{ILS\}\}\(x,y,t\)=\\sum\_\{l=0\}^\{m\-1\}\\frac\{q\_\{l\}\-q\_\{l\-1\}\}\{4\\pi\\lambda\_\{0\}\}\\mathrm\{E\_\{1\}\}\\\!\\left\(\\frac\{r^\{2\}\}\{4\\alpha\_\{0\}\\left\(t\-t\_\{l\}\\right\)\}\\right\),\(8\)whereE1​\(⋅\)\\mathrm\{E\_\{1\}\}\\left\(\\cdot\\right\)is the exponential integral,r=\(x−xB​H​E\)2\+\(y−yB​H​E\)2r=\\sqrt\{\(x\-x\_\{BHE\}\)^\{2\}\+\(y\-y\_\{BHE\}\)^\{2\}\}is the distance from the BHE centered at\(xB​H​E,yB​H​E\)\\left\(x\_\{BHE\},y\_\{BHE\}\\right\),λ0\\lambda\_\{0\}is the bulk thermal conductivity,α0=λ0/\(ρs​Cs\)\\alpha\_\{0\}=\\lambda\_\{0\}/\\left\(\\rho\_\{s\}C\_\{s\}\\right\)denotes the thermal diffusivity andqlq\_\{l\}is the load during time step fromtl−1t\_\{l\-1\}totlt\_\{l\}, withq−1=0q\_\{\-1\}=0andt0=0t\_\{0\}=0\[[5](https://arxiv.org/html/2607.12271#bib.bib3)\]\. The corresponding homogeneous heat diffusion equation with respect to ILS solution is

ρs​Cs​∂T~ILS∂t−λ0​Δ​T~ILS=q​\(t\)​δ​\(x−xB​H​E\)​δ​\(y−yB​H​E\)\.\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}^\{\\mathrm\{ILS\}\}\}\{\\partial t\}\-\\lambda\_\{0\}\\Delta\\widetilde\{T\}^\{\\mathrm\{ILS\}\}=q\(t\)\\delta\\left\(x\-x\_\{BHE\}\\right\)\\delta\\left\(y\-y\_\{BHE\}\\right\)\.\(9\)

#### Finite line source \(FLS\) model\.

To evaluate the thermal behavior in a three\-dimensional homogeneous domain, FLS model represents the boreholes as line sources of finite lengths and assumes a uniform heat rate per unit length along the borehole, an idealized geometry with zero radius and no groundwater effect on the temperature distribution in the underground\[[30](https://arxiv.org/html/2607.12271#bib.bib6)\]\. If the borehole extends from the ground surfacez=Dz=Dto the depthz=D\+Lz=D\+L,

T~FLS​\(x,y,z,t\)=∑l=0m−1ql−ql−14​π​λ0​∫DD\+Lerfc⁡\[R​\(z′\)2​α0​\(t−tl\)\]R​\(z′\)​dz′,\\widetilde\{T\}^\{\\mathrm\{FLS\}\}\(x,y,z,t\)=\\sum\_\{l=0\}^\{m\-1\}\\frac\{q\_\{l\}\-q\_\{l\-1\}\}\{4\\pi\\lambda\_\{0\}\}\\int\_\{D\}^\{D\+L\}\\frac\{\\operatorname\{erfc\}\\\!\\left\[\\dfrac\{R\(z^\{\\prime\}\)\}\{2\\sqrt\{\\alpha\_\{0\}\\left\(t\-t\_\{l\}\\right\)\}\}\\right\]\}\{R\(z^\{\\prime\}\)\}\\,\\mathrm\{d\}z^\{\\prime\},\(10\)whereerfc​\(⋅\)\\mathrm\{erfc\}\\left\(\\cdot\\right\)denotes the complementary error function andR​\(z′\)=r2\+\(z−z′\)2R\(z^\{\\prime\}\)=\\sqrt\{r^\{2\}\+\(z\-z^\{\\prime\}\)^\{2\}\}\. The corresponding homogeneous equation regarding FLS solution is

ρs​Cs​∂T~FLS∂t−λ0​Δ​T~FLS=q​\(t\)​δ​\(x−xB​H​E\)​δ​\(y−yB​H​E\)​χ\[D,D\+L\]​\(z\)\.\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}^\{\\mathrm\{FLS\}\}\}\{\\partial t\}\-\\lambda\_\{0\}\\Delta\\widetilde\{T\}^\{\\mathrm\{FLS\}\}=q\(t\)\\delta\(x\-x\_\{BHE\}\)\\delta\(y\-y\_\{BHE\}\)\\chi\_\{\[D,D\+L\]\}\(z\)\.

#### Moving finite line source \(MFLS\) model\.

The MFLS model extends the FLS solution to a homogeneous medium with groundwater advection governed by a uniform Darcy velocity𝐯=\(vx,vy,0\)\\mathbf\{v\}=\\left\(v\_\{x\},v\_\{y\},0\\right\), which reads

T~MFLS​\(x,y,z,t\)=∑l=0m−1ql−ql−18​π​λ0​∫0t−tl1τ​exp⁡\[−\(x−vx​τ\)2\+\(y−vy​τ\)24​α0′​τ\]​𝒵​\(z,τ\)​dτ,\\widetilde\{T\}^\{\\mathrm\{MFLS\}\}\(x,y,z,t\)=\\sum\_\{l=0\}^\{m\-1\}\\frac\{q\_\{l\}\-q\_\{l\-1\}\}\{8\\pi\\lambda\_\{0\}\}\\int\_\{0\}^\{t\-t\_\{l\}\}\\frac\{1\}\{\\tau\}\\exp\\\!\\left\[\-\\frac\{\(x\-v\_\{x\}\\tau\)^\{2\}\+\(y\-v\_\{y\}\\tau\)^\{2\}\}\{4\\alpha^\{\\prime\}\_\{0\}\\tau\}\\right\]\\mathcal\{Z\}\(z,\\tau\)\\,\\mathrm\{d\}\\tau,\(11\)whereα0′=λ0/\(ρ​C\)\\alpha^\{\\prime\}\_\{0\}=\\lambda\_\{0\}/\\left\(\\rho C\\right\),

𝒵​\(z,τ\)=erf⁡\[D\+L−z2​α0​τ\]−erf⁡\[D−z2​α0​τ\],\\mathcal\{Z\}\(z,\\tau\)=\\operatorname\{erf\}\\\!\\left\[\\frac\{D\+L\-z\}\{2\\sqrt\{\\alpha\_\{0\}\\tau\}\}\\right\]\-\\operatorname\{erf\}\\\!\\left\[\\frac\{D\-z\}\{2\\sqrt\{\\alpha\_\{0\}\\tau\}\}\\right\],anderf​\(⋅\)\\mathrm\{erf\}\\left\(\\cdot\\right\)denotes the error function\[[24](https://arxiv.org/html/2607.12271#bib.bib7)\]\. The corresponding homogeneous advection\-diffusion equation concerning MFLS solution is

ρ​C​∂T~MFLS∂t\+ρf​Cf​𝐯⋅∇T~MFLS−λ0​Δ​T~MFLS=q​\(t\)​δ​\(x−xB​H​E\)​δ​\(y−yB​H​E\)​χ\[D,D\+L\]​\(z\)\.\\rho C\\frac\{\\partial\\widetilde\{T\}^\{\\mathrm\{MFLS\}\}\}\{\\partial t\}\+\\rho\_\{f\}C\_\{f\}\\mathbf\{v\}\\cdot\\nabla\\widetilde\{T\}^\{\\mathrm\{MFLS\}\}\-\\lambda\_\{0\}\\Delta\\widetilde\{T\}^\{\\mathrm\{MFLS\}\}=q\(t\)\\delta\(x\-x\_\{BHE\}\)\\delta\(y\-y\_\{BHE\}\)\\chi\_\{\[D,D\+L\]\}\(z\)\.
All three analytical models are used here as reference solutions under simplified homogeneous assumptions\. They play a key role in removing the dominant singular response so that the remaining correction field can be learned more efficiently by the physics\-informed neural network\. We now derive the governing equations forT~cm\\widetilde\{T\}\_\{c\}^\{m\},m=\{ILS,FLS,MFLS\}m=\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}\. We illustrate the procedure usingT~cILS\\widetilde\{T\}\_\{c\}^\{\\mathrm\{ILS\}\}as an example, and the deduction forT~cFLS\\widetilde\{T\}\_\{c\}^\{\\mathrm\{FLS\}\}andT~cMFLS\\widetilde\{T\}\_\{c\}^\{\\mathrm\{MFLS\}\}can be handled similarly\.

Settingf​\(𝐱,t\)=q​\(t\)​δ​\(x−xB​H​E\)​δ​\(y−yB​H​E\)f\\left\(\\mathbf\{x\},t\\right\)=q\(t\)\\delta\\left\(x\-x\_\{BHE\}\\right\)\\delta\\left\(y\-y\_\{BHE\}\\right\)in \([2](https://arxiv.org/html/2607.12271#S2.E2)\) and subtracting \([9](https://arxiv.org/html/2607.12271#S3.E9)\) from it, we have the correction equation

ρs​Cs​∂T~cI​L​S∂t−∇⋅\(λ​\(x,y\)​∇T~cI​L​S\)=∇⋅\(\(λ​\(x,y\)−λ0\)​∇T~I​L​S\),\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda\\left\(x,y\\right\)\\nabla\\widetilde\{T\}\_\{c\}^\{ILS\}\\right\)=\\nabla\\cdot\\left\(\\left\(\\lambda\\left\(x,y\\right\)\-\\lambda\_\{0\}\\right\)\\nabla\\widetilde\{T\}^\{ILS\}\\right\),\(12\)where the left\-hand side describes the transport of temperature correction and the forcing term represents the heat flow error introduced by plugging the homogeneous solution into the heterogeneous physicsλ​\(𝐱\)\\lambda\\left\(\\mathbf\{x\}\\right\)\. By rewritingλ​\(𝐱\)\\lambda\\left\(\\mathbf\{x\}\\right\)in \([12](https://arxiv.org/html/2607.12271#S3.E12)\) asλ​\(𝐱\)=λB\+λx​xr​e​l\+λy​yr​e​l\\lambda\\left\(\\mathbf\{x\}\\right\)=\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}, whereλB=λ​\(xB​H​E,yB​H​E\)\\lambda\_\{B\}=\\lambda\\left\(x\_\{BHE\},y\_\{BHE\}\\right\),xr​e​l=x−xB​H​Ex\_\{rel\}=x\-x\_\{BHE\}andyr​e​l=y−yB​H​Ey\_\{rel\}=y\-y\_\{BHE\}, we obtain

ρs​Cs​∂T~cI​L​S∂t​\(x,y\)−λx​∂T~cI​L​S∂x​\(x,y\)−\(λB\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~cI​L​S∂x2​\(x,y\)\\displaystyle\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial t\}\\left\(x,y\\right\)\-\\lambda\_\{x\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial x\}\\left\(x,y\\right\)\-\\left\(\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial x^\{2\}\}\\left\(x,y\\right\)−λy​∂T~cI​L​S∂y​\(x,y\)−\(λB\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~cI​L​S∂y2​\(x,y\)\\displaystyle\-\\lambda\_\{y\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial y\}\\left\(x,y\\right\)\-\\left\(\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial y^\{2\}\}\\left\(x,y\\right\)=\\displaystyle=λx​∂T~I​L​S∂x​\(xr​e​l,yr​e​l\)\+\(λB−λ0\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~I​L​S∂x2​\(xr​e​l,yr​e​l\)\\displaystyle\\lambda\_\{x\}\\frac\{\\partial\\widetilde\{T\}^\{ILS\}\}\{\\partial x\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+\\left\(\\lambda\_\{B\}\-\\lambda\_\{0\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}^\{ILS\}\}\{\\partial x^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+λy​∂T~I​L​S∂y​\(xr​e​l,yr​e​l\)\+\(λB−λ0\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~I​L​S∂y2​\(xr​e​l,yr​e​l\)\.\\displaystyle\+\\lambda\_\{y\}\\frac\{\\partial\\widetilde\{T\}^\{ILS\}\}\{\\partial y\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+\\left\(\\lambda\_\{B\}\-\\lambda\_\{0\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}^\{ILS\}\}\{\\partial y^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\.\(13\)Replacingxxbyxr​e​l\+xB​H​Ex\_\{rel\}\+x\_\{BHE\}andyybyyr​e​l\+yB​H​Ey\_\{rel\}\+y\_\{BHE\}in \([13](https://arxiv.org/html/2607.12271#S3.E13)\) and putting the horizontal coordinate origin at the central location of the borehole result in

ρs​Cs​∂T~cI​L​S∂t​\(xr​e​l,yr​e​l\)−λx​∂T~cI​L​S∂x​\(xr​e​l,yr​e​l\)−\(λB\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~cI​L​S∂x2​\(xr​e​l,yr​e​l\)\\displaystyle\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial t\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\-\\lambda\_\{x\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial x\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\-\\left\(\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial x^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)−λy​∂T~cI​L​S∂y​\(xr​e​l,yr​e​l\)−\(λB\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~cI​L​S∂y2​\(xr​e​l,yr​e​l\)\\displaystyle\-\\lambda\_\{y\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial y\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\-\\left\(\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}\_\{c\}^\{ILS\}\}\{\\partial y^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)=\\displaystyle=λx​∂T~I​L​S∂x​\(xr​e​l,yr​e​l\)\+\(λB−λ0\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~I​L​S∂x2​\(xr​e​l,yr​e​l\)\\displaystyle\\lambda\_\{x\}\\frac\{\\partial\\widetilde\{T\}^\{ILS\}\}\{\\partial x\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+\\left\(\\lambda\_\{B\}\-\\lambda\_\{0\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}^\{ILS\}\}\{\\partial x^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+λy​∂T~I​L​S∂y​\(xr​e​l,yr​e​l\)\+\(λB−λ0\+λx​xr​e​l\+λy​yr​e​l\)​∂2T~I​L​S∂y2​\(xr​e​l,yr​e​l\)\.\\displaystyle\+\\lambda\_\{y\}\\frac\{\\partial\\widetilde\{T\}^\{ILS\}\}\{\\partial y\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\+\\left\(\\lambda\_\{B\}\-\\lambda\_\{0\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}\\right\)\\frac\{\\partial^\{2\}\\widetilde\{T\}^\{ILS\}\}\{\\partial y^\{2\}\}\\left\(x\_\{rel\},y\_\{rel\}\\right\)\.\(14\)Thus, the solution of \([14](https://arxiv.org/html/2607.12271#S3.E14)\) represents the temperature correction distribution as a function of the relative position with respect to the borehole\.

Assuming FLS reference model uses the same source term as the heterogeneous problem and applying a similar procedure to obtaining the correction equation for FLS base solution arrives at

ρs​Cs​∂T~cF​L​S∂t−∇⋅\(λ​\(x,y\)​∇T~cF​L​S\)=∇⋅\(\(λ​\(x,y\)−λ0\)​∇T~F​L​S\)\.\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{FLS\}\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda\(x,y\)\\nabla\\widetilde\{T\}\_\{c\}^\{FLS\}\\right\)=\\nabla\\cdot\\left\(\\left\(\\lambda\(x,y\)\-\\lambda\_\{0\}\\right\)\\nabla\\widetilde\{T\}^\{FLS\}\\right\)\.\(15\)Here, the gradient and divergence operators are taken in three\-dimensional space\. The same idea leads to the MFLS\-based correction equation

ρ​C​∂T~cM​F​L​S∂t\+ρf​Cf​𝐯⋅∇T~cM​F​L​S−∇⋅\(λ​\(x,y\)​∇T~cM​F​L​S\)=∇⋅\(\(λ​\(x,y\)−λ0\)​∇T~M​F​L​S\)\.\\rho C\\frac\{\\partial\\widetilde\{T\}\_\{c\}^\{MFLS\}\}\{\\partial t\}\+\\rho\_\{f\}C\_\{f\}\\mathbf\{v\}\\cdot\\nabla\\widetilde\{T\}\_\{c\}^\{MFLS\}\-\\nabla\\cdot\\left\(\\lambda\(x,y\)\\nabla\\widetilde\{T\}\_\{c\}^\{MFLS\}\\right\)=\\nabla\\cdot\\left\(\\left\(\\lambda\(x,y\)\-\\lambda\_\{0\}\\right\)\\nabla\\widetilde\{T\}^\{MFLS\}\\right\)\.\(16\)Via formulatingλ​\(𝐱\)\\lambda\\left\(\\mathbf\{x\}\\right\)asλB\+λx​xr​e​l\+λy​yr​e​l\\lambda\_\{B\}\+\\lambda\_\{x\}x\_\{rel\}\+\\lambda\_\{y\}y\_\{rel\}in \([15](https://arxiv.org/html/2607.12271#S3.E15)\) and \([16](https://arxiv.org/html/2607.12271#S3.E16)\), the correction equation based on arbitrary analytical base model solves the radial distribution around the line source\.

As a result, the ILS, FLS and MFLS corrections embedded in the differential equations obey the same mechanism, where the left\-hand side of the equations transport the thermal correction under the target global physics, while the right\-hand side provides the residual forcing generated by evaluating a homogeneous analytical prior in a heterogeneous conductivity field\. Thus, the neural network does not need to directly deal with the singularity challenge arising from the original problems \([2](https://arxiv.org/html/2607.12271#S2.E2)\) and \([3](https://arxiv.org/html/2607.12271#S2.E3)\)\.

Since the temperature perturbation satisfies \([6](https://arxiv.org/html/2607.12271#S2.E6)\) and the hybrid decomposition givesT~cm=T~−T~m\\widetilde\{T\}\_\{c\}^\{m\}=\\widetilde\{T\}\-\\widetilde\{T\}^\{m\},m∈\{ILS,FLS,MFLS\}m\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}, the corresponding correction field satisfies the following initial condition

T~cm​\(𝐱,0\)=0,𝐱∈Ωs,m∈\{ILS,FLS,MFLS\}\.\\widetilde\{T\}\_\{c\}^\{m\}\(\\mathbf\{x\},0\)=0,\\quad\\mathbf\{x\}\\in\\Omega\_\{s\},\\quad\\mathrm\{m\}\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}\.On the other hand, the ILS and FLS based corrections mainly consist of the calibration of two physical mismatches: one induced by the imbalance between averaged and variable conductivities and one caused by difference between boundary conditions for infinite and finite domains\. Thus,

T~cm​\(𝐱,t\)=−T~m​\(𝐱,t\),𝐱∈∂Ωs,t∈I,m∈\{ILS,FLS\}\.\\displaystyle\\widetilde\{T\}\_\{c\}^\{m\}\(\\mathbf\{x\},t\)=\-\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t\),\\quad\\mathbf\{x\}\\in\\partial\\Omega\_\{s\},\\quad t\\in I,\\quad\\mathrm\{m\}\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\}\\\}\.

### 3\.2Boundary Condition Enforcement

As discussed above, the ILS\-based and FLS\-based temperature correction fields contain both the correction for mismatch from formation properties and the boundary conditions\. The former appears as an interior forcing term in the correction equations, while the latter is inherited from the difference between the idealized boundary assumptions of ILS/FLS models and the prescribed boundary conditions of the original problem\. For clarity, we consider a homogeneous finite\-domain solutionThT\_\{h\}with constant conductivityλ0\\lambda\_\{0\}; then the boundary correction can be defined form∈\{ILS,FLS\}\\mathrm\{m\}\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\}\\\}as

B​\(𝐱,t\)=Th​\(𝐱,t\)−\(T~m​\(𝐱,t\)\+Tu\)=T~h​\(𝐱,t\)−T~m​\(𝐱,t\),B\(\\mathbf\{x\},t\)=T\_\{h\}\(\\mathbf\{x\},t\)\-\\left\(\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t\)\+T\_\{u\}\\right\)=\\widetilde\{T\}\_\{h\}\(\\mathbf\{x\},t\)\-\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t\),which satisfies a homogeneous heat equation,

ρs​Cs​∂B∂t−λ0​Δ​B=0,𝐱∈Ωs,t∈I,\\rho\_\{s\}C\_\{s\}\\frac\{\\partial B\}\{\\partial t\}\-\\lambda\_\{0\}\\Delta B=0,\\qquad\\mathbf\{x\}\\in\\Omega\_\{s\},\\quad t\\in I,\(17\)with the inherited boundary condition

B​\(𝐱,t\)=−T~m​\(𝐱,t\),𝐱∈∂Ωs,t∈I\.B\(\\mathbf\{x\},t\)=\-\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t\),\\qquad\\mathbf\{x\}\\in\\partial\\Omega\_\{s\},\\quad t\\in I\.\(18\)Thus, the boundary correction is a well\-defined component governed by the inconsistency between the boundary conditions\.

In a vanilla PINN formulation, the Dirichlet boundary conditions are implemented using additional penalization terms as a contribution to the loss function\. However, this treatment is particularly inefficient for the present boundary correction problem \([17](https://arxiv.org/html/2607.12271#S3.E17)\) with \([18](https://arxiv.org/html/2607.12271#S3.E18)\)\. This difficulty originates from the fact that the line source response is strongly localized near the borehole at early times, while the thermal perturbation gradually propagates across the domain and interacts with the external boundary at later times\. This produces a pronounced imbalance in both spatial locations and temporal scales, making the approximation sensitive to the sampling of residual and boundary collocation points\. Moreover, statically assigned weights to different loss terms is generally inadequate to accommodate the continuously evolving distribution of the correction field, making the optimization sensitive to hyperparameter tuning and resulting in substantial computational cost\. Therefore, for the FLS\-based correction equation, we adopt a hard\-constrained strategy based on a distance function to match the specified correction boundary values\.

Inspired by\[[2](https://arxiv.org/html/2607.12271#bib.bib2)\], we add a non\-trainable layerNNat the end of the neural network to modify the previous outputT~θ,cF​L​S\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}by

N​T~θ,cF​L​S=−T~¯F​L​S\+ϕ​T~θ,cF​L​S\\displaystyle N\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}=\-\\bar\{\\widetilde\{T\}\}^\{FLS\}\+\\phi\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}to exactly satisfyN​T~θ,cF​L​S​\(𝐱,t\)=−T~F​L​S​\(𝐱,t\)N\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}\\left\(\\mathbf\{x\},t\\right\)=\-\\widetilde\{T\}^\{FLS\}\\left\(\\mathbf\{x\},t\\right\)on𝐱∈∂Ωs\\mathbf\{x\}\\in\\partial\\Omega\_\{s\}, whereT~¯F​L​S\\bar\{\\widetilde\{T\}\}^\{FLS\}is a continuous extension ofT~F​L​S\\widetilde\{T\}^\{FLS\}intoΩs\\Omega\_\{s\}, i\.e\.,T~¯F​L​S\|∂Ωs=T~F​L​S\\bar\{\\widetilde\{T\}\}^\{FLS\}\|\_\{\\partial\\Omega\_\{s\}\}=\\widetilde\{T\}^\{FLS\}andϕ\\phiis an approximate distance function to∂Ωs\\partial\\Omega\_\{s\}which satisfiesϕ​\(𝐱\)=0\\phi\\left\(\\mathbf\{x\}\\right\)=0only for𝐱∈∂Ωs\\mathbf\{x\}\\in\\partial\\Omega\_\{s\}andϕ​\(𝐱\)\>0\\phi\\left\(\\mathbf\{x\}\\right\)\>0for𝐱∈Ωs\\mathbf\{x\}\\in\\Omega\_\{s\}\. The training of the FLS\-based model is to minimize the loss corresponding toN​T~θ,cF​L​SN\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}\. For the sake of simplicity, we considerΩs=\(xmin,xmax\)×\(ymin,ymax\)×\(zmin,zmax\)\\Omega\_\{s\}=\\left\(x\_\{\\min\},x\_\{\\max\}\\right\)\\times\\left\(y\_\{\\min\},y\_\{\\max\}\\right\)\\times\\left\(z\_\{\\min\},z\_\{\\max\}\\right\)here, although approximate distance function can be generalized for more complex geometries in three\-dimensional space\. For the cuboidal domain,ϕ\\phiis assembled from the distances to the six boundary faces

d1=x−xmin,d2=xmax−x,d3=y−ymin,d4=ymax−y,d5=z−zmin,d6=zmax−z,d\_\{1\}=x\-x\_\{\\min\},\\quad d\_\{2\}=x\_\{\\max\}\-x,\\quad d\_\{3\}=y\-y\_\{\\min\},\\quad d\_\{4\}=y\_\{\\max\}\-y,\\quad d\_\{5\}=z\-z\_\{\\min\},\\quad d\_\{6\}=z\_\{\\max\}\-z,using the normalized form of orderpp

ϕ​\(𝐱\)=\(∑i=16di​\(𝐱\)−p\)−1/p,p≥1\.\\phi\(\\mathbf\{x\}\)=\\left\(\\sum\_\{i=1\}^\{6\}d\_\{i\}\(\\mathbf\{x\}\)^\{\-p\}\\right\)^\{\-1/p\},\\qquad p\\geq 1\.SinceT~F​L​S\\widetilde\{T\}^\{FLS\}is well\-defined and continuous insideΩs\\Omega\_\{s\}, we selectT~¯F​L​S​\(𝐱,t\)=T~F​L​S​\(𝐱,t\)\\bar\{\\widetilde\{T\}\}^\{FLS\}\\left\(\\mathbf\{x\},t\\right\)=\\widetilde\{T\}^\{FLS\}\\left\(\\mathbf\{x\},t\\right\),𝐱∈Ω¯s\\mathbf\{x\}\\in\\bar\{\\Omega\}\_\{s\}\. The transfinite interpolation can also be employed to construct the extended boundary function\[[2](https://arxiv.org/html/2607.12271#bib.bib2)\]\.

### 3\.3A Parametric PINN Method

To avoid retraining a separate network for each borehole location and each thermal conductivity realization, we formulate the neural corrector as a parametric PINN\. The model takes the local conductivity parameters and the spatiotemporal coordinates as inputs, and therefore learns a unified correction operator over a family of BHE configurations in linearly heterogeneous media\.

#### Location indicator augmentation\.

We note that explicitly informing the neural networks the distance of sampled data points and collocation points from the location of the borehole benefits the approximation of near\-source distribution\. For the ILS\-based model with no hard\-constrained boundary transformation used, we observe that the boundary condition was less accurately fitted in the direction orthogonal to the conductivity gradient\. This suggests that learning the solution under both absolute and relative coordinate systems using a single neural network structure requires additional efforts\. We therefore introduce normalized global coordinate indicators

x^=x−xminxmax−xmin,y^=y−yminymax−ymin\\displaystyle\\hat\{x\}=\\frac\{x\-x\_\{\\min\}\}\{x\_\{\\max\}\-x\_\{\\min\}\},\\qquad\\hat\{y\}=\\frac\{y\-y\_\{\\min\}\}\{y\_\{\\max\}\-y\_\{\\min\}\}to encode the position of each data or collocation point relative to the computational domain boundary\. These indicators are also used in the FLS correctors with

z^=z−zminzmax−zmin\\displaystyle\\hat\{z\}=\\frac\{z\-z\_\{\\min\}\}\{z\_\{\\max\}\-z\_\{\\min\}\}added to better represent the correction in the interior induced by the boundary condition and to better approach the prescribed boundary behavior near∂Ωs\\partial\\Omega\_\{s\}\. That is to say the relative positions with regard to the borehole which indicates the local geometry and to the boundary interfaces which reflects the global characteristics are both incorporated as feature inputs to the neural net model\.

#### Selection of training points\.

To achieve a higher accuracy of the universal corrector approximation, we use a small number of supervised reference solutions as sparse anchors to calibrate the amplitude and early\-time correction\. The data supervision is integrated with physics\-informed constraints which minimizes a combined residual loss over a set of sampled training points in space\-time\. To cope with the strong spatial and temporal scale imbalance exhibited in the correction field, both the supervised data points and the physics collocation points are selected using a source\-adaptive and time\-adaptive strategy\.

To evaluate the PDE residuals, the horizontal coordinates of each collocation point are sampled adaptively according to its distance from the borehole\. We aim to accurately resolve the steep local gradients near the source by increasing the distribution density of residual points in its vicinity\. Here the horizontal plane is divided into near\-source, intermediate and far\-field regions by

rcore<r≤r1,r1<r≤r2,r\>r2,r\_\{\\mathrm\{core\}\}<r\\leq r\_\{1\},\\qquad r\_\{1\}<r\\leq r\_\{2\},\\qquad r\>r\_\{2\},and fixed portions of collocation pointsα0\\alpha\_\{0\},α1\\alpha\_\{1\},α2\\alpha\_\{2\}are assigned to these regions\. Points inside the immediate line source core are excluded from the PDE residual evaluation to avoid the singularity\. One can alternatively adopt a polar sampling scheme\[[21](https://arxiv.org/html/2607.12271#bib.bib4)\]to control the sampling density by adaptively partitioning the radial and angular directions\.

For the FLS\-based and MFLS\-based models, the vertical coordinate is sampled independently with additional enrichment near the two endpoints of the finite line sourcez=Dz=Dandz=D\+Lz=D\+L\. A mixture sampler with fixed portionsβ0\\beta\_\{0\},β1\\beta\_\{1\}andβ2\\beta\_\{2\}allocated to the endpoints, interior of line source and uniformly sampled full vertical domain\.

To balance early transient behavior and long\-term evolution, the full simulation horizon is divided into several time intervals\{Ik\}k=1Nt\\\{I\_\{k\}\\\}\_\{k=1\}^\{N\_\{t\}\}, and each collocation time is drawn uniformly from these intervals with prescribed probabilities

t∼∑k=1Ntpk​𝒰​\(Ik\),∑k=1Ntpk=1\.t\\sim\\sum\_\{k=1\}^\{N\_\{t\}\}p\_\{k\}\\,\\mathcal\{U\}\(I\_\{k\}\),\\qquad\\sum\_\{k=1\}^\{N\_\{t\}\}p\_\{k\}=1\.
In addition to interior residual points, we generate collocation points to impose initial condition for three analytical\-based models and Dirichlet boundary condition especially for ILS\-based model\.

The supervised data points are sampled following a similar procedure, except that only the BHE coordinates are excluded in the near\-source region\. Logarithmic sampling is employed in time variable to provide a denser distribution at early times\.

Therefore, the spatiotemporal dimension of the total training point set for ILS\-based model is

𝒮s​t=𝒮p​d​e​⋃𝒮i​c​⋃𝒮b​c​⋃𝒮d​a​t​a,\\displaystyle\\mathcal\{S\}\_\{st\}=\\mathcal\{S\}\_\{pde\}\\bigcup\\mathcal\{S\}\_\{ic\}\\bigcup\\mathcal\{S\}\_\{bc\}\\bigcup\\mathcal\{S\}\_\{data\},where

𝒮p​d​e=\{\(xi,yi,ti\)\|\(xj,yj\)∼pr\(x,y\),ti∼∑k=1Ntpk𝒰\(Ik\)\}i=1Np​d​e,\\displaystyle\\mathcal\{S\}\_\{pde\}=\\bigl\\\{\\left\(x\_\{i\},y\_\{i\},t\_\{i\}\\right\)\|\\left\(x\_\{j\},y\_\{j\}\\right\)\\sim p\_\{r\}\\left\(x,y\\right\),t\_\{i\}\\sim\\sum\_\{k=1\}^\{N\_\{t\}\}p\_\{k\}\\,\\mathcal\{U\}\(I\_\{k\}\)\\bigl\\\}\_\{i=1\}^\{N\_\{pde\}\},with

pr​\(x,y\)=α0​𝒰​\(Ωs,rcore<r≤r1\)\+α1​𝒰​\(Ωs,r1<r≤r2\)\+α2​𝒰​\(Ωs,r\>r2\),\\displaystyle p\_\{r\}\\left\(x,y\\right\)=\\alpha\_\{0\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\_\{\\mathrm\{core\}\}<r\\leq r\_\{1\}\}\\right\)\+\\alpha\_\{1\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\_\{1\}<r\\leq r\_\{2\}\}\\right\)\+\\alpha\_\{2\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\>r\_\{2\}\}\\right\),𝒮i​c=\{\(xi,yi,0\)\|\(xi,yi\)∼𝒰​\(Ωs,r\>rc​o​r​e\)\}i=1Ni​c,\\displaystyle\\mathcal\{S\}\_\{ic\}=\\Bigl\\\{\\left\(x\_\{i\},y\_\{i\},0\\right\)\|\\left\(x\_\{i\},y\_\{i\}\\right\)\\sim\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\>r\_\{core\}\}\\right\)\\Bigr\\\}\_\{i=1\}^\{N\_\{ic\}\},𝒮b​c=\{\(xi,yi,0\)\|\(xi,yi\)∼𝒰​\(∂Ωs\)\}i=1Ni​c,\\displaystyle\\mathcal\{S\}\_\{bc\}=\\Bigl\\\{\\left\(x\_\{i\},y\_\{i\},0\\right\)\|\\left\(x\_\{i\},y\_\{i\}\\right\)\\sim\\mathcal\{U\}\\\!\\left\(\\partial\\Omega\_\{s\}\\right\)\\Bigr\\\}\_\{i=1\}^\{N\_\{ic\}\},𝒮d​a​t​a=\{\(xi,yi,ti\)\|\(xj,yj\)∼pr′\(x,y\),ti∼pl​o​g\(t\)\}i=1Nd​a​t​a,\\displaystyle\\mathcal\{S\}\_\{data\}=\\bigl\\\{\\left\(x\_\{i\},y\_\{i\},t\_\{i\}\\right\)\|\\left\(x\_\{j\},y\_\{j\}\\right\)\\sim p^\{\\prime\}\_\{r\}\\left\(x,y\\right\),t\_\{i\}\\sim p\_\{log\}\\left\(t\\right\)\\bigl\\\}\_\{i=1\}^\{N\_\{data\}\},with

pr′​\(x,y\)=α0​𝒰​\(Ωs,0<r≤r1\)\+α1​𝒰​\(Ωs,r1<r≤r2\)\+α2​𝒰​\(Ωs,r\>r2\),\\displaystyle p^\{\\prime\}\_\{r\}\\left\(x,y\\right\)=\\alpha\_\{0\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,0<r\\leq r\_\{1\}\}\\right\)\+\\alpha\_\{1\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\_\{1\}<r\\leq r\_\{2\}\}\\right\)\+\\alpha\_\{2\}\\mathcal\{U\}\\\!\\left\(\\Omega\_\{s,r\>r\_\{2\}\}\\right\),and

plog​\(t\)=1t​ln⁡\(T/tcore\),t∈\[tcore,T\]\.\\displaystyle p\_\{\\log\}\(t\)=\\frac\{1\}\{t\\ln\(T/t\_\{\\mathrm\{core\}\}\)\},\\qquad t\\in\[t\_\{\\mathrm\{core\}\},T\]\.
And the training point set for FLS\-based and MFLS\-based models is formulated as

𝒮s​t=𝒮p​d​e​⋃𝒮i​c​⋃𝒮d​a​t​a,\\displaystyle\\mathcal\{S\}\_\{st\}=\\mathcal\{S\}\_\{pde\}\\bigcup\\mathcal\{S\}\_\{ic\}\\bigcup\\mathcal\{S\}\_\{data\},where

𝒮p​d​e=\{\(xi,yi,zi,ti\)\|\(xj,yj\)∼pr\(x,y\),zi∼pz\(z\),ti∼∑k=1Ntpk𝒰\(Ik\)\}i=1Np​d​e,\\displaystyle\\mathcal\{S\}\_\{pde\}=\\bigl\\\{\\left\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\}\\right\)\|\\left\(x\_\{j\},y\_\{j\}\\right\)\\sim p\_\{r\}\\left\(x,y\\right\),z\_\{i\}\\sim p\_\{z\}\\left\(z\\right\),t\_\{i\}\\sim\\sum\_\{k=1\}^\{N\_\{t\}\}p\_\{k\}\\,\\mathcal\{U\}\(I\_\{k\}\)\\bigl\\\}\_\{i=1\}^\{N\_\{pde\}\},with

pz​\(z\)=β0​𝒰​\(\[D−rz,D\+rz\]∪\[D\+L−rz,D\+L\+rz\]\)\+β1​𝒰​\(\[D,D\+L\]\)\+β2​𝒰​\(\[zmin,zmax\]\),\\displaystyle p\_\{z\}\\left\(z\\right\)=\\beta\_\{0\}\\mathcal\{U\}\\\!\\left\(\[D\-r\_\{z\},D\+r\_\{z\}\]\\cup\[D\+L\-r\_\{z\},D\+L\+r\_\{z\}\]\\right\)\+\\beta\_\{1\}\\mathcal\{U\}\\\!\\left\(\[D,D\+L\]\\right\)\+\\beta\_\{2\}\\mathcal\{U\}\\\!\\left\(\[z\_\{\\min\},z\_\{\\max\}\]\\right\),
𝒮i​c=\{\(xi,yi,zi,0\)\|\(xi,yi\)∼pr​\(x,y\),zi∼pz​\(z\)\}i=1Ni​c,\\displaystyle\\mathcal\{S\}\_\{ic\}=\\bigl\\\{\\left\(x\_\{i\},y\_\{i\},z\_\{i\},0\\right\)\|\\left\(x\_\{i\},y\_\{i\}\\right\)\\sim p\_\{r\}\\left\(x,y\\right\),z\_\{i\}\\sim p\_\{z\}\\left\(z\\right\)\\Bigr\\\}\_\{i=1\}^\{N\_\{ic\}\},
𝒮d​a​t​a=\{\(xi,yi,zi,ti\)\|\(xj,yj\)∼pr′\(x,y\),zi∼pz\(z\),ti∼pl​o​g\(t\)\}i=1Nd​a​t​a\.\\displaystyle\\mathcal\{S\}\_\{data\}=\\bigl\\\{\\left\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\}\\right\)\|\\left\(x\_\{j\},y\_\{j\}\\right\)\\sim p^\{\\prime\}\_\{r\}\\left\(x,y\\right\),z\_\{i\}\\sim p\_\{z\}\\left\(z\\right\),t\_\{i\}\\sim p\_\{log\}\\left\(t\\right\)\\bigl\\\}\_\{i=1\}^\{N\_\{data\}\}\.
We denote the linear parametrization for the local conductivity as a parameter vector𝝁=\(λB,λx,λy\)\\boldsymbol\{\\mu\}=\\left\(\\lambda\_\{B\},\\lambda\_\{x\},\\lambda\_\{y\}\\right\), and sample the training vectors in parametric space𝒫=\{𝝁:λB∈\[λBmin,λBmax\],λx∈\[λxmin,λxmax\],λy∈\[λymin,λymax\]\}\\mathcal\{P\}=\\bigl\\\{\\boldsymbol\{\\mu\}:\\lambda\_\{B\}\\in\\left\[\\lambda\_\{B\}^\{\\min\},\\lambda\_\{B\}^\{\\max\}\\right\],\\lambda\_\{x\}\\in\\left\[\\lambda\_\{x\}^\{\\min\},\\lambda\_\{x\}^\{\\max\}\\right\],\\lambda\_\{y\}\\in\\left\[\\lambda\_\{y\}^\{\\min\},\\lambda\_\{y\}^\{\\max\}\\right\]\\bigl\\\}to evaluate the solutions corresponding to sampled thermal conditions on training points\. For the supervised data points, we generate high\-fidelity reference solutions using the finite difference method for randomly selected conductivity parameters𝝁\(j\)∈𝒫\\boldsymbol\{\\mu\}^\{\(j\)\}\\in\\mathcal\{P\}\. For each sampled parameter realization, the supervised label is computed asT~d​a​t​a,cm​\(𝐱,t;𝝁\(j\)\)=Th​\(𝐱,t;𝝁\(j\)\)−T~m​\(𝐱,t;𝝁\(j\)\)\\widetilde\{T\}\_\{data,c\}^\{m\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)=T\_\{h\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\-\\widetilde\{T\}^\{m\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\),m∈\{ILS,FLS,MFLS\}m\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}, whereTh​\(𝐱,t;𝝁\(j\)\)T\_\{h\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)is the FDM solution with the borehole at the horizontal origin\. For the MFLS model, the advection velocity𝐯=\(vx,vy,0\)\\mathbf\{v\}=\(v\_\{x\},v\_\{y\},0\)is also specified in the analytical baseline\.

Thus, the parametric neural network is trained to approximate a family of correction fields over the parameter space𝒫\\mathcal\{P\}\. Combining the training points and their location indicator augmentation, we define the loss function for ILS\-based model as

ℒILS=wdataILS​ℒdataILS\+wpdeILS​ℒpdeILS\+wicILS​ℒicILS\+wbcILS​ℒbcILS,\\mathcal\{L\}^\{\\mathrm\{ILS\}\}=w\_\{\\mathrm\{data\}\}^\{\\mathrm\{ILS\}\}\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{ILS\}\}\+w\_\{\\mathrm\{pde\}\}^\{\\mathrm\{ILS\}\}\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{ILS\}\}\+w\_\{\\mathrm\{ic\}\}^\{\\mathrm\{ILS\}\}\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{ILS\}\}\+w\_\{\\mathrm\{bc\}\}^\{\\mathrm\{ILS\}\}\\mathcal\{L\}\_\{\\mathrm\{bc\}\}^\{\\mathrm\{ILS\}\},where

ℒdataILS=1Np​Nd​a​t​a​∑j=1Np∑i=1Nd​a​t​a\|T~θ,cI​L​S​\(xi,yi,ti,ri,x^i,y^i;𝝁\(j\)\)−T~d​a​t​a,cI​L​S​\(xi,yi,ti;𝝁\(j\)\)Cscale\|2,\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{ILS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{data\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{data\}\}\\left\|\\frac\{\\widetilde\{T\}^\{ILS\}\_\{\\theta,c\}\(x\_\{i\},y\_\{i\},t\_\{i\},r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\)\-\\widetilde\{T\}^\{ILS\}\_\{data,c\}\\left\(x\_\{i\},y\_\{i\},t\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},ℒpdeILS=1Np​Np​d​e​∑j=1Np∑i=1Np​d​e\|\[ρs​Cs​∂T~θ,cI​L​S∂t−∇⋅\(λ\(j\)​∇T~θ,cI​L​S\)\]​\(xi,yi,ti,ri,x^i,y^i;𝝁\(j\)\)−δ​\(xi\)​δ​\(yi\)Rscale\|2,\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{ILS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{pde\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{pde\}\}\\left\|\\frac\{\\left\[\\rho\_\{s\}C\_\{s\}\\frac\{\\partial\\widetilde\{T\}^\{ILS\}\_\{\\theta,c\}\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda^\{\\left\(j\\right\)\}\\nabla\\widetilde\{T\}^\{ILS\}\_\{\\theta,c\}\\right\)\\right\]\\left\(x\_\{i\},y\_\{i\},t\_\{i\},r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\-\\delta\\left\(x\_\{i\}\\right\)\\delta\\left\(y\_\{i\}\\right\)\}\{R\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},ℒicILS=1Np​Ni​c​∑j=1Np∑i=1Ni​c\|T~θ,cI​L​S​\(xi,yi,0,ri,x^i,y^i;𝝁\(j\)\)Cscale\|2,\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{ILS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{ic\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{ic\}\}\\left\|\\frac\{\\widetilde\{T\}\_\{\\theta,c\}^\{ILS\}\\left\(x\_\{i\},y\_\{i\},0,r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},and

ℒbcILS=1Np​Nb​c​∑j=1Np∑i=1Nb​c\|T~θ,cI​L​S​\(xi,yi,ti,ri,x^i,y^i;𝝁\(j\)\)\+T~I​L​S​\(xi,yi,tj;𝝁\(j\)\)Cscale\|2,\\mathcal\{L\}\_\{\\mathrm\{bc\}\}^\{\\mathrm\{ILS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{bc\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{bc\}\}\\left\|\\frac\{\\widetilde\{T\}^\{ILS\}\_\{\\theta,c\}\(x\_\{i\},y\_\{i\},t\_\{i\},r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\)\+\\widetilde\{T\}^\{ILS\}\\left\(x\_\{i\},y\_\{i\},t\_\{j\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},with predefined normalization scalesCs​c​a​l​eC\_\{scale\},Rs​c​a​l​eR\_\{scale\},λ\(j\)=λB\(j\)\+λx\(j\)​x\+λy\(j\)​y\\lambda^\{\\left\(j\\right\)\}=\\lambda\_\{B\}^\{\\left\(j\\right\)\}\+\\lambda\_\{x\}^\{\\left\(j\\right\)\}x\+\\lambda\_\{y\}^\{\\left\(j\\right\)\}y,ri=xi2\+yi2r\_\{i\}=\\sqrt\{x\_\{i\}^\{2\}\+y\_\{i\}^\{2\}\}, andx^i=xi−xminxmax−xmin\\hat\{x\}\_\{i\}=\\frac\{x\_\{i\}\-x\_\{\\min\}\}\{x\_\{\\max\}\-x\_\{\\min\}\},y^i=yi−yminymax−ymin\\hat\{y\}\_\{i\}=\\frac\{y\_\{i\}\-y\_\{\\min\}\}\{y\_\{\\max\}\-y\_\{\\min\}\}\. And the loss function for FLS\-based model is formulated as

ℒFLS=wdataFLS​ℒdataFLS\+wpdeFLS​ℒpdeFLS\+wicFLS​ℒicFLS,\\mathcal\{L\}^\{\\mathrm\{FLS\}\}=w\_\{\\mathrm\{data\}\}^\{\\mathrm\{FLS\}\}\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{FLS\}\}\+w\_\{\\mathrm\{pde\}\}^\{\\mathrm\{FLS\}\}\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{FLS\}\}\+w\_\{\\mathrm\{ic\}\}^\{\\mathrm\{FLS\}\}\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{FLS\}\},where

ℒdataFLS=1Np​Nd​a​t​a​∑j=1Np∑i=1Nd​a​t​a\|T~θ,cF​L​S​\(xi,yi,zi,ti,ri,x^i,y^i,z^i;𝝁\(j\)\)−T~d​a​t​a,cF​L​S​\(xi,yi,zi,ti;𝝁\(j\)\)Cscale\|2,\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{FLS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{data\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{data\}\}\\left\|\\frac\{\\widetilde\{T\}^\{FLS\}\_\{\\theta,c\}\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\},r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\},\\hat\{z\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\)\-\\widetilde\{T\}^\{FLS\}\_\{data,c\}\\left\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},ℒpdeFLS=\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{FLS\}\}=1Np​Np​d​e​∑j=1Np∑i=1Np​d​e\\displaystyle\\frac\{1\}\{N\_\{p\}N\_\{pde\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{pde\}\}\|\[ρs​Cs​∂N​T~θ,cF​L​S∂t−∇⋅\(λ\(j\)​∇N​T~θ,cF​L​S\)\]​\(xi,yi,zi,ti,ri,x^i,y^i,z^i;𝝁\(j\)\)−δ​\(xi\)​δ​\(yi\)​χ\[D,D\+L\]​\(zi\)Rscale\|2,\\displaystyle\\left\|\\frac\{\\left\[\\rho\_\{s\}C\_\{s\}\\frac\{\\partial N\\widetilde\{T\}^\{FLS\}\_\{\\theta,c\}\}\{\\partial t\}\-\\nabla\\cdot\\left\(\\lambda^\{\\left\(j\\right\)\}\\nabla N\\widetilde\{T\}^\{FLS\}\_\{\\theta,c\}\\right\)\\right\]\\left\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\},r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\},\\hat\{z\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\-\\delta\\left\(x\_\{i\}\\right\)\\delta\\left\(y\_\{i\}\\right\)\\chi\_\{\[D,D\+L\]\}\(z\_\{i\}\)\}\{R\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},and

ℒicFLS=1Np​Ni​c​∑j=1Np∑i=1Ni​c\|N​T~θ,cF​L​S​\(xi,yi,zi,0,ri,x^i,y^i,z^i;𝝁\(j\)\)Cscale\|2,\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{FLS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{ic\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{ic\}\}\\left\|\\frac\{N\\widetilde\{T\}\_\{\\theta,c\}^\{FLS\}\\left\(x\_\{i\},y\_\{i\},z\_\{i\},0,r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\},\\hat\{z\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},withri=xi2\+yi2\+zi2r\_\{i\}=\\sqrt\{x\_\{i\}^\{2\}\+y\_\{i\}^\{2\}\+z\_\{i\}^\{2\}\},z^i=zi−zminzmax−zmin\\hat\{z\}\_\{i\}=\\frac\{z\_\{i\}\-z\_\{\\min\}\}\{z\_\{\\max\}\-z\_\{\\min\}\}\. Similarly for MFLS\-based model, we have

ℒMFLS=wdataMFLS​ℒdataMFLS\+wpdeMFLS​ℒpdeMFLS\+wicMFLS​ℒicMFLS,\\mathcal\{L\}^\{\\mathrm\{MFLS\}\}=w\_\{\\mathrm\{data\}\}^\{\\mathrm\{MFLS\}\}\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{MFLS\}\}\+w\_\{\\mathrm\{pde\}\}^\{\\mathrm\{MFLS\}\}\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{MFLS\}\}\+w\_\{\\mathrm\{ic\}\}^\{\\mathrm\{MFLS\}\}\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{MFLS\}\},whereℒdataMFLS\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{MFLS\}\}follows the same definition asℒdataFLS\\mathcal\{L\}\_\{\\mathrm\{data\}\}^\{\\mathrm\{FLS\}\},

ℒpdeMFLS=1Np​Np​d​e​∑j=1Np∑i=1Np​d​e\|ℛ​\(xi,yi,zi,ti,ri;𝝁\(j\)\)−δ​\(xi\)​δ​\(yi\)​χ\[D,D\+L\]​\(zi\)Rscale\|2,\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{pde\}\}^\{\\mathrm\{MFLS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{pde\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{pde\}\}\\left\|\\frac\{\\mathcal\{R\}\\left\(x\_\{i\},y\_\{i\},z\_\{i\},t\_\{i\},r\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\-\\delta\\left\(x\_\{i\}\\right\)\\delta\\left\(y\_\{i\}\\right\)\\chi\_\{\[D,D\+L\]\}\(z\_\{i\}\)\}\{R\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\},with

ℛ=ρ​C​∂T~θ,cM​F​L​S∂t\+ρf​Cf​𝐯⋅∇T~θ,cM​F​L​S−∇⋅\(λ\(j\)​∇T~θ,cM​F​L​S\),\\displaystyle\\mathcal\{R\}=\\rho C\\frac\{\\partial\\widetilde\{T\}^\{MFLS\}\_\{\\theta,c\}\}\{\\partial t\}\+\\rho\_\{f\}C\_\{f\}\\mathrm\{\\mathbf\{v\}\}\\cdot\\nabla\\widetilde\{T\}^\{MFLS\}\_\{\\theta,c\}\-\\nabla\\cdot\\left\(\\lambda^\{\\left\(j\\right\)\}\\nabla\\widetilde\{T\}^\{MFLS\}\_\{\\theta,c\}\\right\),and

ℒicMFLS=1Np​Ni​c​∑j=1Np∑i=1Ni​c\|T~θ,cM​F​L​S​\(xi,yi,zi,0,ri,x^i,y^i,z^i;𝝁\(j\)\)Cscale\|2\.\\mathcal\{L\}\_\{\\mathrm\{ic\}\}^\{\\mathrm\{MFLS\}\}=\\frac\{1\}\{N\_\{p\}N\_\{ic\}\}\\sum\_\{j=1\}^\{N\_\{p\}\}\\sum\_\{i=1\}^\{N\_\{ic\}\}\\left\|\\frac\{\\widetilde\{T\}\_\{\\theta,c\}^\{MFLS\}\\left\(x\_\{i\},y\_\{i\},z\_\{i\},0,r\_\{i\},\\hat\{x\}\_\{i\},\\hat\{y\}\_\{i\},\\hat\{z\}\_\{i\};\\boldsymbol\{\\mu\}^\{\\left\(j\\right\)\}\\right\)\}\{C\_\{\\mathrm\{scale\}\}\}\\right\|^\{2\}\.

### 3\.4Spatial and Temporal Superposition

As discussed above, the trained parametric corrector provides the thermal response for a single BHE with unit line load in a linearly heterogeneous computational domain\. For an arbitrary borehole located at\(xB​H​Ek,yB​H​Ek\)\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\), the unit load hybrid response is written as

Tθ∗,km​\(𝐱,t;𝝁k\)=Tu\+T~m​\(𝐱,t;λ0\)\+T~θ∗,cm​\(x−xB​H​Ek,y−yB​H​Ek,z,t,r,x^,y^,z^;𝝁k\),T\_\{\\theta^\{\*\},k\}^\{m\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}\_\{k\}\\right\)=T\_\{u\}\+\\widetilde\{T\}^\{m\}\(\\mathbf\{x\},t;\\lambda\_\{0\}\)\+\\widetilde\{T\}^\{m\}\_\{\\theta^\{\*\},c\}\(x\-x\_\{BHE\}^\{k\},y\-y\_\{BHE\}^\{k\},z,t,r,\\hat\{x\},\\hat\{y\},\\hat\{z\};\\boldsymbol\{\\mu\}\_\{k\}\),where𝝁k=\(λ​\(xB​H​Ek,yB​H​Ek\),λx​\(xB​H​Ek,yB​H​Ek\),λy​\(xB​H​Ek,yB​H​Ek\)\)\\boldsymbol\{\\mu\}\_\{k\}=\\left\(\\lambda\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\),\\lambda\_\{x\}\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\),\\lambda\_\{y\}\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\)\\right\),θ∗\\theta^\{\*\}are the updated trainable parameters of the network andTkm​\(𝐱,t;𝝁k\)T\_\{k\}^\{m\}\\left\(\\mathbf\{x\},t;\\boldsymbol\{\\mu\}\_\{k\}\\right\)denotes the predicted soil temperature with regard to analytical modelm∈\{ILS,FLS,MFLS\}m\\in\\\{\\mathrm\{ILS\},\\mathrm\{FLS\},\\mathrm\{MFLS\}\\\}\. For the two\-dimensional ILS model, thezz\-dependent terms are omitted\.

With the assumption that the soil thermal properties do not depend on the temperature distribution, the governing heat diffusion equation, advection\-diffusion equation and the analytical solutions are linear with respect to the source load for a prescribed conductivity field\. Due to the additivity of the thermal energy, the response to multiple BHEs and time\-varying loads can be reconstructed by spatial and temporal superposition\[[5](https://arxiv.org/html/2607.12271#bib.bib3)\]\. Let the load history be discretized intommpiecewise\-constant time intervals, with

tl=l​Δ​t,l=0,…,m,t\_\{l\}=l\\Delta t,\\qquad l=0,\\ldots,m,and letqk,lq\_\{k,l\}denote the constant load of thekk\-th borehole during\[tl−1,tl\)\[t\_\{l\-1\},t\_\{l\}\)\. The predicted temperature response at time steptm−1<t≤tmt\_\{m\-1\}<t\\leq t\_\{m\}concerningKKBHEs in the computation domain with a certain conductivityλ​\(x,y\)\\lambda\\left\(x,y\\right\)is calculated as

Tθ∗m​\(𝐱,t;λ​\(x,t\)\)=Tu\+∑k=1K∑l=0m−1\(T~k,lm​\(𝐱,t;λ0\)\+\(ql−ql−1\)​T~θ∗,cm​\(x−xB​H​Ek,y−yB​H​Ek,z,t−tl,r,x^,y^,z^;𝝁k\)\)\\displaystyle T^\{m\}\_\{\\theta^\{\*\}\}\\left\(\\mathbf\{x\},t;\\lambda\\left\(x,t\\right\)\\right\)=T\_\{u\}\+\\sum\_\{k=1\}^\{K\}\\sum\_\{l=0\}^\{m\-1\}\\left\(\\widetilde\{T\}^\{m\}\_\{k,l\}\\left\(\\mathbf\{x\},t;\\lambda\_\{0\}\\right\)\+\\left\(q\_\{l\}\-q\_\{l\-1\}\\right\)\\widetilde\{T\}^\{m\}\_\{\\theta^\{\*\},c\}\\left\(x\-x\_\{BHE\}^\{k\},y\-y\_\{BHE\}^\{k\},z,t\-t\_\{l\},r,\\hat\{x\},\\hat\{y\},\\hat\{z\};\\boldsymbol\{\\mu\}\_\{k\}\\right\)\\right\)Here,T~k,lm\\widetilde\{T\}^\{m\}\_\{k,l\}denotes the single step analytical approximation concerningqk,l−qk,l−14​π​λ0\\frac\{q\_\{k,l\}\-q\_\{k,l\-1\}\}\{4\\pi\\lambda\_\{0\}\}in \([8](https://arxiv.org/html/2607.12271#S3.E8)\), \([10](https://arxiv.org/html/2607.12271#S3.E10)\) and \([11](https://arxiv.org/html/2607.12271#S3.E11)\)\. The superposition principle allows the trained parametric universal correctorT~θ∗,cm\\widetilde\{T\}^\{m\}\_\{\\theta^\{\*\},c\}to be used for multi\-BHE field simulations without retraining the neural network for various borehole locations and thermal heterogeneities\.

## 4Case Study

In this section, to illustrate the effectiveness of the proposed parametric hybrid\-PINN corrector method, we present several numerical results in two and three spatial dimensions for BHE\-induced global thermal problem with line sources\. After the training process is finished, we examine the accuracy of predicted hybrid solution by computing the relativeL2L^\{2\}andL∞L^\{\\infty\}errors between the learned solution and reference solution provided by finite difference method\. To be specific, we use the grid points for visualization as the test points to evaluate the errors of a certain analytical model based corrections at a selected time stepttby

ℰ∞,hyb,mt=‖Tθ∗m−Th‖∞‖Th‖∞,ℰ2,hyb,mt=‖Tθ∗m−Th‖2‖Th‖2,\\displaystyle\\mathcal\{E\}\_\{\\infty,\\mathrm\{hyb\},m\}^\{t\}=\\frac\{\\\|T^\{m\}\_\{\\theta^\{\*\}\}\-T\_\{h\}\\\|\_\{\\infty\}\}\{\\\|T\_\{h\}\\\|\_\{\\infty\}\},\\qquad\\mathcal\{E\}\_\{2,\\mathrm\{hyb\},m\}^\{t\}=\\frac\{\\\|T^\{m\}\_\{\\theta^\{\*\}\}\-T\_\{h\}\\\|\_\{2\}\}\{\\\|T\_\{h\}\\\|\_\{2\}\},respectively, where

‖Th‖∞=max1≤i≤Nt​e​s​t⁡\|Th​\(𝐱i,t\)\|,‖Th‖2=1Nt​e​s​t​∑i=1Nt​e​s​tTh2​\(𝐱i,t\)\.\\\|T\_\{h\}\\\|\_\{\\infty\}=\\max\\limits\_\{1\\leq i\\leq N\_\{test\}\}\|T\_\{h\}\\left\(\\mathbf\{x\}\_\{i\},t\\right\)\|,\\qquad\\\|T\_\{h\}\\\|\_\{2\}=\\sqrt\{\\frac\{1\}\{N\_\{test\}\}\\sum\_\{i=1\}^\{N\_\{test\}\}T\_\{h\}^\{2\}\\left\(\\mathbf\{x\}\_\{i\},t\\right\)\}\.Similarly, we demonstrate the improved effect of the hybrid corrector method compared to the analytical solutions using a bulk conductivity or local conductivities

T~m​\(𝐱,t\)=∑k=1K∑l=0m−1T~k,lm​\(𝐱,t\),T~l​om​\(𝐱,t\)=∑k=1K∑l=0m−1T~k,l,l​om​\(𝐱,t\),\\displaystyle\\widetilde\{T\}^\{m\}\\left\(\\mathbf\{x\},t\\right\)=\\sum\_\{k=1\}^\{K\}\\sum\_\{l=0\}^\{m\-1\}\\widetilde\{T\}^\{m\}\_\{k,l\}\\left\(\\mathbf\{x\},t\\right\),\\qquad\\widetilde\{T\}^\{m\}\_\{lo\}\\left\(\\mathbf\{x\},t\\right\)=\\sum\_\{k=1\}^\{K\}\\sum\_\{l=0\}^\{m\-1\}\\widetilde\{T\}^\{m\}\_\{k,l,lo\}\\left\(\\mathbf\{x\},t\\right\),whereT~k,lm\\widetilde\{T\}^\{m\}\_\{k,l\}represents the single step response regarding modelmmmodulated byqk,l−qk,l−14​π​λ0\\frac\{q\_\{k,l\}\-q\_\{k,l\-1\}\}\{4\\pi\\lambda\_\{0\}\}, whileT~k,l,l​om\\widetilde\{T\}^\{m\}\_\{k,l,lo\}the response modulated byqk,l−qk,l−14​π​λ​\(xB​H​Ek,yB​H​Ek\)\\frac\{q\_\{k,l\}\-q\_\{k,l\-1\}\}\{4\\pi\\lambda\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\)\}andαk=λ​\(xB​H​Ek,yB​H​Ek\)/\(ρs​Cs\)\\alpha\_\{k\}=\\lambda\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\)/\\left\(\\rho\_\{s\}C\_\{s\}\\right\)orαk=λ​\(xB​H​Ek,yB​H​Ek\)/\(ρ​C\)\\alpha\_\{k\}=\\lambda\\left\(x\_\{BHE\}^\{k\},y\_\{BHE\}^\{k\}\\right\)/\\left\(\\rho C\\right\)\. The corresponding errors are computed by

ℰ∞,ana,mt=‖Tu\+T~m−Th‖∞‖Th‖∞,ℰ2,ana,mt=‖Tu\+T~m−Th‖2‖Th‖2\.\\displaystyle\\mathcal\{E\}\_\{\\infty,\\mathrm\{ana\},m\}^\{t\}=\\frac\{\\\|T\_\{u\}\+\\widetilde\{T\}^\{m\}\-T\_\{h\}\\\|\_\{\\infty\}\}\{\\\|T\_\{h\}\\\|\_\{\\infty\}\},\\qquad\\mathcal\{E\}\_\{2,\\mathrm\{ana\},m\}^\{t\}=\\frac\{\\\|T\_\{u\}\+\\widetilde\{T\}^\{m\}\-T\_\{h\}\\\|\_\{2\}\}\{\\\|T\_\{h\}\\\|\_\{2\}\}\.SubstitutingT~m\\widetilde\{T\}^\{m\}byT~l​om\\widetilde\{T\}\_\{lo\}^\{m\}gives the errors for locally estimated analytical solutions\. For each selected time step, the FDM reference solution is interpolated onto the same evaluation grid as the analytical and hybrid predictions\. In general, withΔ​x=Δ​y=Δ​z=1​m\\Delta x=\\Delta y=\\Delta z=1\\,\\mathrm\{m\}, we setNt​e​s​t=81×121=9801N\_\{test\}=81\\times 121=9801orNt​e​s​t=114×121=13794N\_\{test\}=114\\times 121=13794points for vertical slice visualization depending on the scenario considered, andNt​e​s​t=81×81=6561N\_\{test\}=81\\times 81=6561for each horizontal slice\.

Throughout all numerical experiments for all three investigated base models, the neural corrector is implemented as a fully connected network with five hidden layers, 128 neurons per layer andtanh\\tanhactivations\. The linear output layer returns the correction fieldT~θ∗,cI​L​S\\widetilde\{T\}^\{ILS\}\_\{\\theta^\{\*\},c\},T~θ∗,cM​F​L​S\\widetilde\{T\}^\{MFLS\}\_\{\\theta^\{\*\},c\}orT~θ∗,cF​L​S\\widetilde\{T\}^\{FLS\}\_\{\\theta^\{\*\},c\}to generate hard\-constrainedN​T~θ∗,cF​L​SN\\widetilde\{T\}^\{FLS\}\_\{\\theta^\{\*\},c\}\. The network model is trained using AdamW optimizer with a learning rate of10−410^\{\-4\}, weight decay10−510^\{\-5\}and gradient clipping with threshold0\.50\.5\.

In the following examples, we consider two thermal conductivity scenarios: one varying linearly along thexxdirection and the other along both thexxandyydirections, as shown in Fig\.[1](https://arxiv.org/html/2607.12271#S4.F1)\.

![Refer to caption](https://arxiv.org/html/2607.12271v1/conductivity_x.png)\(a\)Scenario 1
![Refer to caption](https://arxiv.org/html/2607.12271v1/conductivity_xy.png)\(b\)Scenario 2

Figure 1:Thermal conductivity fields considered for testing ILS, FLS and MFLS based models\.The remaining parameters featuring the soil and fluid for the ILS/FLS\-based and MFLS\-based models are listed in Tables[1](https://arxiv.org/html/2607.12271#S4.T1)and[2](https://arxiv.org/html/2607.12271#S4.T2)\.

Table 1:Physical property settings for demonstration examples used for training and inference procedure of parametric hybrid\-PINN corrector method based on ILS and FLS models\.Table 2:Physical property settings for demonstration examples used for training and inference procedure of parametric hybrid\-PINN corrector method based on MFLS model\.#### Example 1\.\(ILS\-based Model\)

We choose a square domainΩs=\[−10,70\]​m×\[−10,70\]​m\\Omega\_\{s\}=\\left\[\-10,70\\right\]\\mathrm\{m\}\\times\\left\[\-10,70\\right\]\\mathrm\{m\}and a simulation time span of3030years, and we consider a case where2525BHEs are arranged as a regular5×55\\times 5array with a spacing1010m\. In terms of model training, we derive the bulk thermal conductivityλ0\\lambda\_\{0\}as the average of the conductivities at all locations of the boreholes\. To test the proposed model on the two scenarios, we introduce an interpolation parameterγ∼𝒰​\(0,1\)\\gamma\\sim\\mathcal\{U\}\(0,1\)to be sampled\. Letλmintrain=0\.5\\lambda\_\{\\min\}^\{\\mathrm\{train\}\}=0\.5,λmaxtrain=10\\lambda\_\{\\max\}^\{\\mathrm\{train\}\}=10and the conductivity gradient parameters are defined as

λx​\(γ\)=\(1−γ\)​λmaxtrain−λmintrainxmax−xmin\+γ​λmaxtrain−λmintrain2​\(xmax−xmin\),λy​\(γ\)=γ​λmaxtrain−λmintrain2​\(ymax−ymin\)\.\\lambda\_\{x\}\(\\gamma\)=\(1\-\\gamma\)\\frac\{\\lambda\_\{\\max\}^\{\\mathrm\{train\}\}\-\\lambda\_\{\\min\}^\{\\mathrm\{train\}\}\}\{x\_\{\\max\}\-x\_\{\\min\}\}\+\\gamma\\frac\{\\lambda\_\{\\max\}^\{\\mathrm\{train\}\}\-\\lambda\_\{\\min\}^\{\\mathrm\{train\}\}\}\{2\(x\_\{\\max\}\-x\_\{\\min\}\)\},\\qquad\\lambda\_\{y\}\(\\gamma\)=\\gamma\\frac\{\\lambda\_\{\\max\}^\{\\mathrm\{train\}\}\-\\lambda\_\{\\min\}^\{\\mathrm\{train\}\}\}\{2\(y\_\{\\max\}\-y\_\{\\min\}\)\}\.The selected conductivity parameter is then derived by combining

\(λB​\(γ\),λx​\(γ\),λy​\(γ\)\),\\displaystyle\\left\(\\lambda\_\{B\}\\left\(\\gamma\\right\),\\lambda\_\{x\}\\left\(\\gamma\\right\),\\lambda\_\{y\}\\left\(\\gamma\\right\)\\right\),where

λB​\(γ\)=λ​\(xB​H​E,yB​H​E;γ\)=λmintrain\+λx​\(γ\)​\(xB​H​E−xmin\)\+λy​\(γ\)​\(yB​H​E−ymin\)\.\\displaystyle\\lambda\_\{B\}\\left\(\\gamma\\right\)=\\lambda\(x\_\{BHE\},y\_\{BHE\};\\gamma\)=\\lambda\_\{\\min\}^\{\\mathrm\{train\}\}\+\\lambda\_\{x\}\(\\gamma\)\(x\_\{BHE\}\-x\_\{\\min\}\)\+\\lambda\_\{y\}\(\\gamma\)\(y\_\{BHE\}\-y\_\{\\min\}\)\.The training samples are generated in the predefined computational domain under the sampled conductivity conditions with the borehole fixed at the horizontal origin\. We generate reference solutions under3232and6464conductivity conditions for training ILS and FLS/MFLS based models, respectively\. The spatial coordinates of physics collocation points and the interior data points are both chosen using a three\-layer adaptive strategy around the BHE\. We use

rcore=2​m,r1=6​m,r2=20​m,\\displaystyle r\_\{\\mathrm\{core\}\}=2~\\mathrm\{m\},\\qquad r\_\{1\}=6~\\mathrm\{m\},\\qquad r\_\{2\}=20~\\mathrm\{m\},with portions

α0=0\.50,α1=0\.35,α2=0\.15\.\\displaystyle\\alpha\_\{0\}=0\.50,\\qquad\\alpha\_\{1\}=0\.35,\\qquad\\alpha\_\{2\}=0\.15\.The PDE collocation times are sampled from four stratified intervals\[1,24\]\[1,24\],\[24,72\]\[24,72\],\[72,180\]\[72,180\]and\[180,360\]\[180,360\]months, with equal probability0\.250\.25for each interval\. During the inference, we adopt the source load profile from\[[5](https://arxiv.org/html/2607.12271#bib.bib3)\]where the energy extraction distribution for 1 year is given and we apply the monthly adapted pattern as shown in Fig\.[2](https://arxiv.org/html/2607.12271#S4.F2)repeatedly over the 30\-year period\. And a uniform heat rate per unit length is assumed along each active borehole segment\. The trained universal corrector is then evaluated for each BHE using its own local parameter vector𝝁k​\(k=1,…,25\)\\boldsymbol\{\\mu\}\_\{k\}\\left\(k=1,\.\.\.,25\\right\)regrading two scenarios\. The multi\-BHE field response is then obtained by spatial superposition over the2525BHEs and temporal superposition over the piecewise\-constant monthly load history\.

![Refer to caption](https://arxiv.org/html/2607.12271v1/extracted_energy_profile.png)\(a\)Extracted energy per month\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/heat_transfer_rate_profile.png)\(b\)Heat transfer rate per month of the considered system consisting of 25 BHEs with identical lengthsL=100L=100m and heat transfer rates\.

Figure 2:Annual source load profile repeatedly used in the simulations\.The compared distributions by numerical methods, analytical models and hybrid correctors are shown in Fig\.[3](https://arxiv.org/html/2607.12271#S4.F3)and Fig\.[4](https://arxiv.org/html/2607.12271#S4.F4), while the error metrics are summarized in Tables[3](https://arxiv.org/html/2607.12271#S4.T3)and[4](https://arxiv.org/html/2607.12271#S4.T4)\. For brevity, we only present the results for December of the 15th year and provide those for other selected time steps in the Appendix\. As is shown, the ILS error exhibits a structured asymmetry along the conductivity gradient direction and significant mismatch near the boundary, and the corrected solution largely removes the bias and recovers the residual introduced by the soil heterogeneity\.

![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_x_15.png)Figure 3:The first row: temperature distributions in December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), ILS analytical formula \(ILS solution\) and hybrid\-PINN corrector based on ILS solution \(Predicted solution\); the second row: the difference distributions of ILS solution and predicted solution with respect to the ground truth\.Table 3:L2L^\{2\}andL∞L^\{\\infty\}relative errors of ILS analytical solutions using the averaged thermal conductivity over the domain, local conductivities at each borehole location and the hybrid\-PINN predicted solution in Decembers of the55th,1515th and3030th years for Scenario 1\.Table 4:L2L^\{2\}andL∞L^\{\\infty\}relative errors of ILS analytical solutions using the averaged thermal conductivity over the domain, local conductivities at each borehole location and the hybrid\-PINN predicted solution in Decembers of the55th,1515th and3030th years for Scenario 2\.![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_xy_15.png)Figure 4:The first row: temperature distributions in December of the1515th year for Scenario 2, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), ILS analytical formula \(ILS solution\) and hybrid\-PINN corrector based on ILS solution \(Predicted solution\); the second row: the difference distributions of ILS solution and predicted solution with respect to the ground truth\.
#### Example 2\.\(FLS\-based Model\)

In this example, the same computational domain, simulation period, BHE arrangement and input sampling are applied as those of ILS\-based model\. Here, FLS models each BHE as a vertical line source with zero radius starting at depthD=10D=10m and extending toD\+L=110D\+L=110m\. The vertical coordinates of training points are sampled from the source endpoint regions \(z∈\[4,16\]⋃\[104,116\]\)z\\in\[4,16\]\\bigcup\[104,116\]\)m, the active source segment \(z∈\[10,110\]z\\in\[10,110\]m\) and the complete vertical plane with probabilitiesβ0=0\.45\\beta\_\{0\}=0\.45,β1=0\.35\\beta\_\{1\}=0\.35andβ2=0\.2\\beta\_\{2\}=0\.2, respectively\.

The results tested in December of the 15th year are presented in Figs\.[5](https://arxiv.org/html/2607.12271#S4.F5)and[6](https://arxiv.org/html/2607.12271#S4.F6), where the temperature fields are visualized on they=30y=30andy=xy=xcross sections for Scenario 1 and Scenario 2, respectively\. And the relative error for all three tested time steps are reported in Tables[5](https://arxiv.org/html/2607.12271#S4.T5)and[6](https://arxiv.org/html/2607.12271#S4.T6)\.

Table 5:L2L^\{2\}andL∞L^\{\\infty\}relative error comparison among the global FLS solution, local FLS solution and hybrid\-PINN predictions in Decembers of the55th,1515th and3030th years for Scenario 1, computed over the vertical cross section aty=30y=30\.Table 6:L2L^\{2\}andL∞L^\{\\infty\}relative error comparison among the global FLS solution, local FLS solution and hybrid\-PINN predictions in Decembers of the55th,1515th and3030th years for Scenario 2, computed over the vertical cross section aty=xy=x\.![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_x_15.png)Figure 5:The first row: temperature distributions on the vertical cross section aty=30y=30in December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), FLS analytical formula \(FLS solution\) and hybrid\-PINN corrector based on FLS solution \(Predicted solution\); the second row: the difference distributions of FLS solution and predicted solution with respect to the ground truth on the same cross section\.![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_xy_15.png)Figure 6:The first row: temperature distributions on the vertical cross section aty=xy=xin December of the1515th year for Scenario 2, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), FLS analytical formula \(FLS solution\) and hybrid\-PINN corrector based on FLS solution \(Predicted solution\); the second row: the difference distributions of FLS solution and predicted solution with respect to the ground truth on the same cross section\. The horizontal coordinate is labeled by the projection onto thexxaxis\.
#### Example 3\.\(MFLS\-based Model\)

We employ the same computational setup and training configuration as FLS\-based model to this case\. We choose𝐯=\(0,2,0\)\\mathrm\{\\mathbf\{v\}\}=\\left\(0,2,0\\right\)m/year for training and inference in this test\. To examine the effects of conductivity and groundwater advection, we evaluate the predicted temperature fields on two representative vertical sections which are aligned with the conductivity gradient and groundwater flow direction respectively for both scenarios\. The simulation results for Scenario 1 are exhibited in Fig\.[7](https://arxiv.org/html/2607.12271#S4.F7)and Fig\.[8](https://arxiv.org/html/2607.12271#S4.F8)and Scenario 2 in Fig\.[9](https://arxiv.org/html/2607.12271#S4.F9)and Fig\.[10](https://arxiv.org/html/2607.12271#S4.F10)\. The corresponding horizontal slices atz=60z=60are provided in Appendix\.

![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_15_ymid.png)Figure 7:The first row: temperature distributions on the vertical cross section aty=30y=30in December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector based on MFLS solution \(Predicted solution\); the second row: the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_15_xmid.png)Figure 8:The first row: temperature distributions on the vertical cross section atx=30x=30in December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector based on MFLS solution \(Predicted solution\); the second row: the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.Tables[7](https://arxiv.org/html/2607.12271#S4.T7)and[8](https://arxiv.org/html/2607.12271#S4.T8)show the comparison between the accuracy performance of global MFLS solution, local MFLS solution and predicted approximation on the cross sectionsy=30y=30,y=xy=xfor Scenarios 1 and 2, respectively\.

![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_15_ymid.png)Figure 9:The first row: temperature distributions on the vertical cross section aty=xy=xin December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector based on MFLS solution \(Predicted solution\); the second row: the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\. The horizontal coordinate is labeled by the projection onto thexxaxis\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_15_xmid.png)Figure 10:The first row: temperature distributions on the vertical cross section atx=30x=30in December of the1515th year for Scenario 1, derived from the finite difference solver of mesh sizeh=1h=1m \(Ground truth\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector based on MFLS solution \(Predicted solution\); the second row: the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.Table 7:L2L^\{2\}andL∞L^\{\\infty\}relative error comparison among the global MFLS solution, local MFLS solution and hybrid\-PINN predictions in Decembers of the55th,1515th and3030th years for Scenario 1, computed over the vertical cross section aty=30y=30\.Table 8:L2L^\{2\}andL∞L^\{\\infty\}relative error comparison among the global MFLS solution, local MFLS solution and hybrid\-PINN predictions in Decembers of the55th,1515th and3030th years for Scenario 2, computed over the vertical cross section aty=xy=x\.

## 5Conclusion

In this paper, we develop a parametric hybrid\-PINN method, based on ILS, FLS and MFLS analytical models, for thermal simulation of borehole heat exchangers, i\.e\. singular line sources, in heterogeneous subsurface\. By decomposing the global soil temperature distribution into an analytical baseline and a learned correction, we remove the challenge for the network to accurately capture the singularity and avoid the need to approach the delta function using traditional regularization tools\. For the ILS and FLS based models, the goal of the correction component is to characterize both mismatches induced by formation properties and boundary conditions on account of the idealized assumptions of the analytical models on the soil homogeneity and domain boundary\. The MFLS\-based model mainly accounts for the former under an open boundary setting\.

We propose a parametric physics\-informed neural network to approximate the temperature correction governed by the correction equation and initial condition, with Dirichlet boundary treatment enforced for conduction\-dominated cases\. To inform the network the relative position of any input coordinate with respect to the borehole source and boundary interfaces, we include location indicator functions as additional feature inputs to the network, endowing it with capability of learning local and global behaviors simultaneously\. In addition, we adopt a source\-centered sampling strategy for generating physics collocation and data points to handle the correction distribution imbalance over the space and time scales\. The physics\-informed and data\-supervised loss function is constructed via evaluating all the sampled gradient conductivity conditions on these adaptively placed training points\. The inference efficiency is attained through applying the superposition principles to the trained universal corrector\. We perform numerical experiments to show the effectiveness of the proposed method on top of all ILS, FLS and MFLS models in a multi\-BHE domain over a long simulation period\. As demonstrated in presented results, the tested examples can be solved with acceptable accuracy by the parametric hybrid\-PINN universal corrector with moderate number of supervised data samples\.

Future research will extend to providing more comprehensive models for BHE simulation by coupling the heat transport inside the borehole to soil temperature response\. The fast and differentiable nature of the learned corrector also makes it promising for PDE\-constrained optimization, including borehole layout design, installation length, load profile planning and long\-term operational scheduling\. Another future development is to explore frameworks which can solve more complex soil property fields and include a broader set of design and operation parameters\. As the parametric space becomes higher\-dimensional, more expressive learning architectures including neural operators, encoder\-decoder structures and generative deep learning models, may be required to construct scalable surrogate models\. Incorporating these directions in future work will further improve the proposed methodology, enabling it to serve as a more extensive tool for the efficient simulation and optimization of geothermal energy systems\.

## Conflict of Interest

The authors declare no conflicts of interest\.

## Acknowledgement

The work presented in this paper was supported by the German Federal Ministry of Economics and Climate Protection \(BMWK\) within the scope of the research project OptiGeoS \(01256700/1\)\.

## References

- \[1\]R\. Al\-Khoury, T\. Kölbel, and R\. Schramedei\(2010\)Efficient numerical modeling of borehole heat exchangers\.Computers & Geosciences36\(10\),pp\. 1301–1315\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[2\]S\. Berrone, C\. Canuto, M\. Pintore, and N\. Sukumar\(2023\)Enforcing Dirichlet boundary conditions in physics\-informed neural networks and variational physics\-informed neural networks\.Heliyon9\(8\),pp\. e18820\.Cited by:[§3\.2](https://arxiv.org/html/2607.12271#S3.SS2.p3.2),[§3\.2](https://arxiv.org/html/2607.12271#S3.SS2.p3.22)\.
- \[3\]L\. Cassina, L\. Laloui, and A\. F\. Rotta Loria\(2022\)Thermal interactions among vertical geothermal borehole fields\.Renewable Energy194,pp\. 1204–1220\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[4\]M\. de Paly, J\. Hecht\-Méndez, M\. Beck, P\. Blum, A\. Zell, and P\. Bayer\(2012\)Optimization of energy extraction for closed shallow geothermal systems using linear programming\.Geothermics43,pp\. 57–65\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[5\]M\. de Paly, J\. Hecht\-Méndez, M\. Beck, P\. Blum, A\. Zell, and P\. Bayer\(2012\)Optimization of energy extraction for closed shallow geothermal systems using linear programming\.Geothermics43,pp\. 57–65\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1),[§3\.1](https://arxiv.org/html/2607.12271#S3.SS1.SSS0.Px1.p1.11),[§3\.4](https://arxiv.org/html/2607.12271#S3.SS4.p2.1),[§4](https://arxiv.org/html/2607.12271#S4.SS0.SSS0.Px1.p1.18)\.
- \[6\]H\.\-J\. G\. Diersch, D\. Bauer, W\. Heidemann, W\. Rühaak, and P\. Schätzl\(2011\)Borehole heat exchangers: a finite element formulation and the integration of transient heat conduction and advection processes in the 3D porous medium\. part I: steady\-state description\.Computers & Geosciences37\(1\),pp\. 1–13\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[7\]J\. S\. Figueira, A\. García\-Gil, A\. Vieira, A\. K\. Michopoulos, D\. P\. Boon, F\. Loveridge, F\. Cecinato, G\. Götzl, J\. Epting, K\. Zosseder, M\. Bloemendal, M\. Woods, P\. Christodoulides, P\. J\. Vardon, S\. P\. Borg, S\. E\. Poulsen, and T\. R\. Andersen\(2024\)Shallow geothermal energy systems for district heating and cooling networks: review and technological progression through case studies\.Renewable Energy236,pp\. 121436\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[8\]Z\. Gao, T\. Tang, L\. Yan, and T\. Zhou\(2024\)Failure\-informed adaptive sampling for PINNs, part II: combining with re\-sampling and subset simulation\.Communications on Applied Mathematics and Computation6\(3\),pp\. 1720–1741\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[9\]Z\. Gao, L\. Yan, and T\. Zhou\(2023\)Failure\-informed adaptive sampling for PINNs\.SIAM Journal on Scientific Computing45\(4\),pp\. A1971–A1994\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[10\]S\. Halilovic, F\. Böttcher, S\. C\. Kramer, M\. D\. Piggott, K\. Zosseder, and T\. Hamacher\(2022\)Well layout optimization for groundwater heat pump systems using the adjoint approach\.Energy Conversion and Management268,pp\. 116033\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[11\]S\. Halilovic, F\. Böttcher, K\. Zosseder, and T\. Hamacher\(2023\)Optimization approaches for the design and operation of open\-loop shallow geothermal systems\.Advances in Geosciences62,pp\. 57–66\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[12\]S\. Halilovic, F\. Böttcher, K\. Zosseder, and T\. Hamacher\(2023\)Optimizing the spatial arrangement of groundwater heat pumps and their well locations\.Renewable Energy217,pp\. 119148\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[13\]S\. Halilovic, F\. Böttcher, K\. Zosseder, and T\. Hamacher\(2024\)Spatial analysis of thermal groundwater use based on optimal sizing and placement of well doublets\.Energy304,pp\. 132058\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[14\]S\. Halilovic, L\. Odersky, and T\. Hamacher\(2022\)Integration of groundwater heat pumps into energy system optimization models\.Energy238,pp\. 121607\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[15\]S\. Halilovic\(2026\)Deep learning methods for the simulation and optimization of shallow geothermal energy systems\.EarthArXiv eprints\.External Links:[Document](https://dx.doi.org/10.31223/X58N3K)Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p3.1)\.
- \[16\]O\. Heinzel, S\. Halilovic, T\. Hamacher, and M\. Ulbrich\(2026\)Optimization of closed\-loop shallow geothermal systems using analytical models\.arXiv preprint arXiv:2603\.24957\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[17\]G\. Hellström\(1989\)Duct ground heat storage model: manual for computer code\.Technical reportDepartment of Mathematical Physics, University of Lund,Lund, Sweden\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[18\]T\. Hu, B\. Jin, and Z\. Zhou\(2023\)Solving elliptic problems with singular sources using singularity splitting deep Ritz method\.SIAM Journal on Scientific Computing45\(4\),pp\. A2043–A2074\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[19\]T\. Hu, B\. Jin, and Z\. Zhou\(2024\)Solving Poisson problems in polygonal domains with singularity enriched physics\-informed neural networks\.SIAM Journal on Scientific Computing46\(4\),pp\. C369–C398\.External Links:[Document](https://dx.doi.org/10.1137/23M1601195)Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[20\]X\. Huang, H\. Liu, B\. Shi, Z\. Wang, K\. Yang, Y\. Li, M\. Wang, H\. Chu, J\. Zhou, F\. Yu, B\. Hua, B\. Dong, and L\. Chen\(2022\)A universal PINNs method for solving partial differential equations with a point source\.InProceedings of the Thirty\-First International Joint Conference on Artificial Intelligence,IJCAI\-22,pp\. 3839–3846\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[21\]E\. Kara and P\. Stinis\(2025\)Physics\-informed DeepONet coupled with FEM for convective transport in porous media with sharp Gaussian sources\.Note:arXiv:2508\.19847Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1),[§3\.3](https://arxiv.org/html/2607.12271#S3.SS3.SSS0.Px2.p2.3)\.
- \[22\]M\.\-C\. Lai, C\.\-C\. Chang, W\.\-S\. Lin, W\.\-F\. Hu, and T\.\-S\. Lin\(2022\)A shallow Ritz method for elliptic problems with singular sources\.Journal of Computational Physics469,pp\. 111547\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p4.1)\.
- \[23\]M\. Li and A\. C\. K\. Lai\(2015\)Review of analytical models for heat transfer by vertical ground heat exchangers \(GHEs\): a perspective of time and space scales\.Applied Energy151,pp\. 178–191\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[24\]N\. Molina\-Giraldo, P\. Blum, K\. Zhu, P\. Bayer, and Z\. Fang\(2011\)A moving finite line source model to simulate borehole heat exchangers with groundwater advection\.International Journal of Thermal Sciences50\(12\),pp\. 2506–2513\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1),[§3\.1](https://arxiv.org/html/2607.12271#S3.SS1.SSS0.Px3.p1.3)\.
- \[25\]M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis\(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p3.1)\.
- \[26\]J\. Randow, S\. Chen, K\. Lubashevsky, S\. Thiel, T\. Reinhardt, K\. Rink, R\. Grimm, A\. Bucher, O\. Kolditz, and H\. Shao\(2022\)Modeling neighborhood\-scale shallow geothermal energy utilization: a case study in berlin\.Geothermal Energy10,pp\. 1\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[27\]D\. Romanov and B\. Leiss\(2022\)Geothermal energy at different depths for district heating and cooling of existing and future building stock\.Renewable and Sustainable Energy Reviews167,pp\. 112727\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[28\]H\. Sadeghi, R\. Jalali, and R\. M\. Singh\(2024\)A review of borehole thermal energy storage and its integration into district heating systems\.Renewable and Sustainable Energy Reviews192,pp\. 114236\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p1.1)\.
- \[29\]H\. Shao, P\. Hein, D\. Bertermann, O\. Kolditz, and W\. Rühaak\(2016\)Numerical modeling of borehole heat exchangers in OpenGeoSys and its validation using field experiment data\.Environmental Earth Sciences75\(15\),pp\. 1147\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.
- \[30\]H\. Y\. Zeng, N\. R\. Diao, and Z\. H\. Fang\(2002\)A finite line\-source model for boreholes in geothermal heat exchangers\.Heat Transfer–Asian Research31\(7\),pp\. 558–567\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1),[§3\.1](https://arxiv.org/html/2607.12271#S3.SS1.SSS0.Px2.p1.2)\.
- \[31\]C\. Zhang, Y\. Wang, Y\. Liu, X\. Kong, and Q\. Wang\(2018\)Computational methods for ground thermal response of multiple borehole heat exchangers: a review\.Renewable Energy127,pp\. 461–473\.Cited by:[§1](https://arxiv.org/html/2607.12271#S1.p2.1)\.

## Appendix

![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_x_5.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_x_30.png)\(b\)December of the3030th simulated year\.

Figure 11:Temperature distributions for Scenario 1, derived from the finite difference solver \(Ground truth,h=1h=1m\), ILS analytical formula \(ILS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of ILS solution and predicted solution with respect to the ground truth\.![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_xy_5.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/ils_xy_30.png)\(b\)December of the3030th simulated year\.

Figure 12:Temperature distributions for Scenario 2, derived from the finite difference solver \(Ground truth,h=1h=1m\), ILS analytical formula \(ILS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of ILS solution and predicted solution with respect to the ground truth\.![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_x_5.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_x_30.png)\(b\)December of the3030th simulated year\.

Figure 13:Temperature distributions on the vertical cross sectiony=30y=30for Scenario 1, derived from the finite difference solver \(Ground truth,h=1h=1m\), FLS analytical formula \(FLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of FLS solution and predicted solution with respect to the ground truth on the same cross section\.![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_xy_5.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/fls_xy_30.png)\(b\)December of the3030th simulated year\.

Figure 14:Temperature distributions on the vertical sectiony=xy=xfor Scenario 2, derived from the finite difference solver \(Ground truth,h=1h=1m\), FLS analytical formula \(FLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of FLS solution and predicted solution with respect to the ground truth on the same cross section\. The horizontal coordinate is labeled by the projection onto thexxaxis\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_5_ymid.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_30_ymid.png)\(b\)December of the3030th simulated year\.

Figure 15:Temperature distributions on the vertical cross sectiony=30y=30for Scenario 1, derived from the finite difference solver \(Ground truth,h=1h=1m\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_5_ymid.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_30_ymid.png)\(b\)December of the3030th simulated year\.

Figure 16:Temperature distributions on the vertical sectiony=xy=xfor Scenario 2, derived from the finite difference solver \(Ground truth,h=1h=1m\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\. The horizontal coordinate is labeled by the projection onto thexxaxis\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_5_zmid.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_x_30_zmid.png)\(b\)December of the3030th simulated year\.

Figure 17:Temperature distributions on the horizontal sectionz=60z=60for Scenario 1, derived from the finite difference solver \(Ground truth,h=1h=1m\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_5_zmid.png)\(a\)December of the55th simulated year\.
![Refer to caption](https://arxiv.org/html/2607.12271v1/mfls_xy_30_zmid.png)\(b\)December of the3030th simulated year\.

Figure 18:Temperature distributions on the horizontal sectionz=60z=60for Scenario 2, derived from the finite difference solver \(Ground truth,h=1h=1m\), MFLS analytical formula \(MFLS solution\) and hybrid\-PINN corrector \(Predicted solution\) and the difference distributions of MFLS solution and predicted solution with respect to the ground truth on the same cross section\.