Estimation and prediction of shallow ground source heat resources subjected to complex soil and atmospheric boundary conditions

5th generation district heating and cooling networks operating at near ground temperature offer a low-cost, zero- carbon energy solution. Detailed understanding and accurate estimation of ground behaviour for its heat storage and recharge potential are of paramount importance for the success of such networks. In this paper, an advanced modelling tool, based on a coupled Thermal-Hydraulic (TH) modelling framework, is presented to calculate and predict temperature and soil-moisture behaviour of a shallow ground under complex atmospheric, temperature and hydraulic boundary conditions. Atmospheric data e.g., solar radiation, rainfall, humidity, air temperature, wind velocity is considered together with subsurface soil data to investigate thermal and hydraulic responses of the ground, and its individual soil layers. Furthermore, a transient method for estimating shallow ground source heat (SGSH) resources is proposed based on the simulated temperature and saturation distributions of the ground. The model is applied to predict the long-term ground temperature and saturation level of a test site located in Warwickshire County, UK. The total heat content per unit area and the annual/seasonal/monthly net heat content per unit area of the site are predicted for a five-year period. The total heat content of the sandy clay layer varied between 2.32 and 11.6 MJ/m 2 , silty clay from 34.0 to 50.5 MJ/m 2 , and mudstone from 50.7 to 55.0 MJ/m 2 . A parametric sensitivity study is also conducted to investigate the effects of soil types and hydraulic drainage conditions on the ground heat supply potential, and it revealed that the spatial and temporal distri- butions of ground heat is significantly affected by the underlying soils. This study highlights the influences of atmospheric conditions and coupled ground processes, and the parameters that should be considered for designing a 5th generation low-temperature heat network.


Introduction
Great ambitions on Climate Change, such as the Glasgow Climate Pact and the Paris Agreement to limit the global temperature rise to 1.5 • C above pre-industrial levels [1], and the European Green Deal to achieve no net emissions of greenhouse gases (GHG) by 2050 [2], require significant reduction of fossil fuel usage and promotion of zero-carbon, renewable energy solutions. Energy sector is one of the largest emitters of GHG, dominated by space heating and cooling demand. The International Energy Agency (IEA) reported that half of the global energy consumption is associated with heating homes, industries, and other applications [3]. Heating and cooling accounted for 40% of the overall energy use in Europe [4] and heating alone accounted for 44% of the total energy consumption (in the year 2017) in the United Kingdom [5]. Meanwhile, the warmer summer periods have accelerated the cooling demand in recent years. Therefore, decarbonisation of the heating and cooling systems is essential to reduce GHG emission and to achieve a net-zero energy sector in Europe and across the globe.
Numerous efforts are ongoing to decarbonise heating, including exploitation of geothermal energy, solar energy, waste heat, combined heat and power (CHP) plants, etc. [6][7][8][9]. These low carbon heating options can be integrated, individually or collectively, to form heat networks, i.e., district heating and cooling networks (DHC). The DHCs are recognised for their effectiveness in reducing fossil fuel or primary energy consumption and associated emissions [8]. The conventional heating networks, based on high temperature CHP systems, yet requires significant input from fossil fuels, but the 4th and 5th generation district heating and cooling systems operate at the low temperature and even the near-ambient/ground temperature, and offer a low-cost, zero carbon supply of heating and cooling with flexibility of integrating other renewable energy sources [10].
Geothermal energy, such as shallow ground source heat (SGSH) energy [11,12], is becoming increasingly popular for heating and cooling of residential, commercial, and public buildings [13], and it is aligned with the development of the 5th generation heating and cooling networks.
To support design and installation of ground source heat pump (GSHP) for exploiting SGSH for a 5th generation network, the amount of thermal energy available in a ground should be thoroughly investigated. If the reserves of SGSH cannot meet the demand, GSHP can be combined with other systems, such as solar thermal panels, to create hybrid systems for greater performance [13][14][15]. To estimate ground temperatures accurately, it is essential to understand dynamics of ground thermal behaviour which varies seasonally, with soil types and covers, and with atmospheric conditions, such as precipitation, condensation, evaporation etc. [11,12,[16][17][18].
Reserves of SGSH are closely related to ground temperature distribution and heat fluxes through the soil-atmosphere interface [11,12]. Research has studied temperature profiles and heat fluxes in shallow grounds. Oliver et al. [16] measured soil temperature and soil heat fluxes from several sites in the U.K. and Syria, and they reported that daily soil heat fluxes could be estimated reasonably from soil temperatures. Mihalakakou et al. [19] proposed a mathematical model based on the transient heat conduction equation and the energy balance equation, which was used to predict ground surface temperature of bare and short-grass covered soils. They validated their model against experimental data obtained from Athens and Dublin and found that the surface temperatures were sensitive to wind speed, soil absorptivity, evaporation, and relative humidity. Staniec and Nowak [20] predicted annual soil temperature distributions via a mathematical model, which was also based on the transient heat conduction equation with energy balance at the soil surface boundary. From sensitivity studies they revealed that the model was mostly sensitive to changes of air temperature. Bryś et al. [11] gathered nine-years of temperature and radiation data of bare soil and grassy soil in Poland. Annual and monthly variabilities of positive heat flux density at various locations of a subsurface shallow depth soil layer (SSDSL) were analysed based on the energy balance equation. In a follow-up study, Bryś et al. [12] investigated the dynamics of diurnal, monthly, seasonal, and annual positive and negative heat fluxes in SSDSL based on ten-year measured data of bare soil and grassy soil in Poland. Factors affecting the energy balance equation and thermal reserves of SSDSL were discussed.
Existing models predicting soil temperature are often confronted with limitations. For example, they are applicable only to certain meteorological conditions [21], simplified in nature, and often neglect coupled thermo-hydraulic behaviours of a ground [19][20][21]. In the aforementioned studies, soil temperature and heat fluxes were analysed primarily from the perspective of energy balance. However, soil is a three-phase system consisting of solid aggregates, pore-water, and pore-air. Heat and moisture exchange occur constantly through the soil-atmosphere interface and soil thermal behaviours depend on responses of these constituents, both individually and collectively. For example, thermal properties of soils are influenced by its water content [17,18]. Saturated soils with water-filled pores facilitate higher thermal conductivity than partially or unsaturated soils, where the pore spaces are occupied by both water and air. Akrouch et al. [22] conducted experimental investigation into sandy soil and observed that the thermal conductivity reduced from 2.65 W/m/K in saturated condition to 0.9 W/m/K in dry soil. To obtain reasonable temperature and heat flux, as well as to calculate the key GSHP parameters, e.g., loop length, heat exchange efficiency, and heat loss, it is essential to investigate the coupled thermo-hydraulic (TH) behaviour of the ground including the complexities of variable climatic conditions [23]. This is a pre-requisite for accurate design and estimation of an effective and efficient low-temperature heat network.
In this paper, dynamic thermal and hydraulic behaviours of a ground that is subjected to complex atmospheric conditions are investigated for shallow ground source heat (SGSH) utilization. A previously developed thermo-hydraulic model [24][25][26], to study heat and moisture transport and their coupled interactions in soil, is further improved within the scope of this study. The advanced model is applied to estimate theoretical heat resources of a potential site intended to support a shallow or 5th generation heating and cooling network. To the authors' knowledge, such an advanced and on-purpose modelling tool is rare in the literature. The theoretical heat capacity is defined as the seasonal/annual heat reserve of an intact ground in response to the local climatic conditions. The model represents the soil-atmosphere interface as a surface boundary which allows complex climatic data to be incorporated in a robust manner. For example, seasonal variations of solar radiation, rainfall, humidity, air temperature, wind velocity are included in the model to predict their influence on thermal behaviour, e.g., ground heat content and temperature, and moisture behaviour, e.g., saturation level. It considers heat and moisture dissipation not only through the soil-atmosphere surface boundary but also through the surrounding soil boundaries, and therefore represents realistic field conditions. Based on the simulated temperature and saturation distributions, a transient method to estimate the capacity of SGSH is proposed. Developments of the 5th generation networks are at their early stages, and the detailed understanding of the ground behaviour in terms of its heat storage and recharge potential are of paramount importance for the success of this network [4]. This is also the focus of this paper.
The structure of this paper is as follows. Firstly, the theoretical background of the coupled model is presented. Secondly, a validation exercise is carried out to demonstrate the accuracy of the model for predicting temperature and moisture distributions of a site selected from literatures. After that, the model is applied to predict and assess the longterm ground temperature, saturation level, and heat content of a test site located in Warwickshire County, UK. Finally, a parametric sensitivity study is conducted to highlight the effects of soil types and hydraulic drainage conditions on the heat storage capacities of the ground.

Theoretical and numerical formulations
The formulations presented here are for the coupled TH behaviour of an unsaturated rigid porous soil in shallow ground, from the ground surface to around 10 m depth in Earth. The flow of heat and moisture in the rigid soil constitutes the two main processes considered in the model. The heat transport formulation is based on the conservation of energy, whereas the moisture flow is considered via mass conservation. Two primary variables, pore-water pressure u l and temperature T, are used to formulate the governing equations. The formulated ground surface boundary considers the soil-atmosphere interface in terms of heat and moisture flows.

Moisture transfer
Moisture flow within unsaturated soils can be descried as a twophase process, compromising of both liquid and vapour flows. The general equation for the moisture flow can be described as: where θ l is the volumetric liquid content, θ a is the volumetric air content, ∂V is the incremental volume, t is the time, ∇ is the gradient operator, ρ l is the density of the liquid water, ρ v is the density of water vapour, v l is the velocity of the liquid water, and v v is the equivalent velocity of vapour. The volumetric liquid and air contents θ l and θ a can be expressed in terms of the porosity ϕ and the degree of saturation of pore water S l : where S a is the degree of saturation of pore air.
The incremental volume ∂V is defined as the summation of the incremental solid volume ∂V s and void ratio e: Pressure head and elevation head are considered as the main mechanisms to contribute to liquid water flow, and Darcy's law is applied: where γ l is the unit weight of the liquid, y is the elevation, and K l is the unsaturated hydraulic conductivity, which is modelled using the Brooks and Corey [27] Model herein: in which K ls is the saturated hydraulic conductivity, θ ls is the saturated water content, and η is the shape parameter.
The Van Genuchten [28] Model is used to characterize the Soil Water Characteristic Curve (SWCC) of soils, and there is: where h is the pressure head, h g is the scale parameter, and n is the material coefficient.
Philip and de Vries [29] proposed a vapour flow law in unsaturated soil, where the velocity of vapour is described by: ∇ρ v (8) in which D atms is the molecular diffusivity of vapour through air, v v is a mass flow factor, τ v is a tortuosity factor, and ∇ρ v is the spatial vapour density gradient. The following expression proposed by Krischer and Rohnalter [30] is used to obtain the molecular diffusivity of vapour through air: in which u a is the pore-air pressure that was fixed as 1atm herein. An expression proposed by Partington [31] is adopted for the mass flow factor: where u v is the partial pressure of vapour, which can be calculated using the following thermodynamic relationship under assuming 'ideal gas' behaviour [32]: where R v is the specific gas constant for water vapour (461.5 J/kg/K). A thermodynamic relationship was proposed by Edlefsen and Anderson [33], known as the psychrometric law, showing the density of water vapour can be given as: where ρ 0 is the density of saturated water vapour, g is the gravitational constant, and φ is the capillary potential and can be defined as: In addition to that, the density of saturated water vapour can be obtained by the relationship proposed by Ewen and Thomas [34] as follows:

Heat transfer
The primary mechanisms of heat transfer, including conduction and convection are considered in the coupled TH model. Conservation of energy dictates the temporal derivative of the heat content Ω is equal to the spatial derivative of the heat flux Q expressed as: For an unsaturated soil, its heat content per unit volume is the sum of the product of the heat storage capacity with temperature change and the contribution of enthalpy characterised by the latent heat of vaporisation [34]. There is: (16) where L is the latent heat of vaporisation. The following equation was adopted by Ewen and Thomas [34] to calculate the heat capacity of an unsaturated soil H c at a reference temperature T r : where C ps , C pl , and C pv are the specific heat capacities of the solid, liquid water, and vapour, respectively, and ρ s is the density of the solid.
Following Thomas and He [24], the heat flux Q can be defined by considering three components of heat transportation as shown below: in which the first term describes the thermal conduction in accordance with Fourier's law, the second term describes the latent heat flow associated with vapour movement, and the third term describes the heat convection in terms of the liquid phase movement and vapour phase associated with a vapour pressure gradient.
Following the approach of Thomas and Rees [25] the coefficient of thermal conductivity for an unsaturated soil λ T can be defined as follows: where λ s , λ w , and λ a is the thermal conductivity coefficient corresponding to the solid, water, and air component of the soil, respectively, and χ s , χ w , and χ a is the volume fraction corresponding to the solid, water, and air component of the soil, respectively, which can be presented by the following expressions:

Heat flow at the soil-atmosphere interface
In the coupled TH model, four major mechanisms of heat exchange between the ground surface and atmosphere are considered, that is, shortwave heat radiation, longwave heat radiation, sensible heat radiation, and latent heat radiation. The effects of organisms, such as vegetation, are not considered. Therefore, the energy balance equation at the soil-atmosphere interface can be presented as follows [11,19,35]: in which H is the total heat radiation flux absorbed or emitted at the soil surface, H Absorbed SW is the absorbed shortwave radiation flux, H Absorbed LW is the longwave radiation flux absorbed at the ground surface, H Emitted LW is the longwave radiation flux emitted at the ground surface, H SEN is the sensible heat radiation flux, and H LE the latent heat radiation flux.
The absorbed shortwave radiation flux can be expressed as follows [36]: where ε SW is the shortwave reflection factor associated with the ground surface type, and H SW is the shortwave radiation flux striking the ground surface (shortwave solar radiation).
The longwave radiation flux absorbed at the ground surface from the atmosphere is given as [37,38]: where σ is the Stefan-Boltzmann constant (5.67 × 10 − 8 W/m 2 /K 4 ), C cloud is the fractional cloud cover coefficient (0 for clear sky and 1 for totally overcast), T air is the absolute temperature of the air adjacent to the ground surface (ambient air temperature), and ε A LW is the long-wave emissivity of the air at ground level (non-dimensional) which can be obtained by: The longwave radiation flux emitted at the ground surface can be calculated by the Stefan-Boltzmann's law [39]: where ε LW is the long-wave emissivity of the ground (non-dimensional).
The sensible heat radiation flux is given as the following equation [36]: where ρ a is the air density, C pa is the specific heat capacity of air, u y is the wind speed at an elevation y (wind speed), and C y is a drag coefficient (non-dimensional).
The latent heat flux due to evaporation is calculated as follows [36]: where L is the latent heat of vaporisation, and E is the evaporation flux (saturated state), and it can be obtained by the following equation [36,40]: in which q is the specific humidity of the soil at the ground surface, and q air is the specific humidity of air.
The value obtained by Equation (30) is equal to the so-called potential evaporation which can be described as the upper limit or the maximum rate of evaporation from a surface when it is saturated [41]. As the upper surface of soil would enter the unsaturated state, the following modification is made to E and there is [42]: where E a is the actual evaporation flux (unsaturated state), h ground is the relative humidity of the ground surface, and h air the relative humidity of air at the ground surface (air relative humidity).

Moisture flow at the soil-atmosphere interface
Assuming the net moisture flux at the surface is a function of hydrological process only, the ground moisture flux can be presented as follows [35,43]: where Q m is the net moisture mass flux at the ground surface, P is the precipitation mass flux (rainfall), and R is the run-off. Therefore, based on Equations 23-32, there are five climatic variables needed in the coupled TH model to model the ground surface boundary in terms of heat and moisture exchanges between the ground surface and the atmosphere, including ambient air temperature T air , shortwave solar radiation H SW , air relative humidity h air , wind speed u y , and rainfall P.

Numerical formulation
The coupled TH model presented in this paper has been developed within the framework of a bespoke thermo-hydraulic-chemicalmechanical (THCM) modelling code, namely, COMPASS (COde of Modelling PArtially Saturated Soils) [24].
In the model, the governing transport equations are expressed in terms of the primary variables, i.e., pore-water pressure u l and temperature T. Equation (1) and Equation (15) can be expressed as follows: The C and K terms represent storage and flux, respectively, and detailed in the Supplementary Information (SI). Within the framework of COMPASS, the Galerkin finite-element method [44] is adopted to spatially discretize the model equations, e.g., Equations (33) and (34), and an implicit mid-interval backward difference time-stepping algorithm is employed for temporal discretisation. The discretized system of linear equations is solved iteratively using a predictor-corrector algorithm [45]. Furthermore, regarding the ground surface boundary, which is prescribed by two fluxes, one representing thermal interactions and the other representing hydraulic surface interactions, a sequential algorithm was developed and implemented within COMPASS code [26].

Heat content and heat storage of shallow ground
As presented earlier, the heat content for an unsaturated soil per unit volume Ω (J/m 3 ) in Equation (15) denotes the intrinsic energy for the ground in a certain status: Assuming the ground consists of multi-layers of soils whose heat content can be distinct but only varies along with the depth, the layerwise summation method is used to approximately calculate the total heat content per unit area for the ground, which is denoted by Λ total (J/ m 2 ).
As shown in Fig. 1, there are N inspected points (IP) evenly taken along the domain depth. For the ith inspected point, it is assumed that the soil's density, specific heat capacity, saturation, and temperature are constant, namely, the heat content per unit volume Ω i (J/m 3 ) is constant based on Equation (35) in the thickness of Δy (m) (for IP 1 on the top of the domain and IP N on the bottom of the domain, the corresponding thickness is Δy/2). Therefore, the total heat content per unit area Λ total (J/m 2 ) for the whole domain can be estimated by the following equation: where: in which Y is the total depth of the ground domain (m). It should be noted that when the number of inspected points is great enough, the calculated total heat content per unit area Λ total can be close to the actual value of the ground. A similar method can be adopted to assess the total heat content per unit area for each soil layer.
Under the influences of the varied climatic conditions, the total heat content per unit area for the whole domain changes with time, which can be hypothetically illustrated by Fig. 2. In a calendar year with four seasons, the value of Λ total would decrease with the decline of temperature and solar radiation in Autumn and Winter, while increase with the rise of temperature and solar radiation in Spring and Summer. Correspondingly, in contrast to the starting time t = 0 in a year with the initial total heat content per unit area Λ 0 total , a +ve value of Λ total in Summer and Autumn means heat accumulation, while a -ve value of Λ total indicates heat is loss in the rest of the year.
Due to heat being accumulated and lost in a year, the annual net heat content per unit area N HC (J/m 2 ) for the whole domain can also be estimated by the summation method as follows: in which M is the number of Λ j total in a year and M = 365, and Λ j total is the heat content per unit area for the whole domain in the jth day (J/ m 2 ). A similar method can be employed to obtain the annual net heat content for each soil layer and seasonal/monthly net heat content per unit area for the whole domain and each soil layer.
Assuming the area of the ground model domain in known as A (m 2 ), then the annual net heat storage N HS (J) for the ground domain can be obtained as follows:

Model validation
Based on the literature review and sourced data, the validation of the coupled TH model with the ground surface boundary was carried out by comparing the simulated evaporation at the soil atmospheric interface against the measurements of an in-situ experimental study conducted in Southern France by Enrique et al. [46] and Calvet et al. [47]. Sensors were installed to monitor the surface water exchanges and the soil investigation was conducted to measure the soil properties at the site. The thermal and hydraulic behaviour of the ground at the site throughout the monitoring period (1 January 1995 to 7 March 1995) was simulated.

Initial and boundary conditions
A 2D simulation was carried out for validation. As shown in Fig. 3, the model domain has a depth of 4 m and a width of 0.1 m, and it is made up of five layers. The selected domain was discretized into 5718 triangular elements connected by 3091 nodes. Finer mesh size was set near the ground surface than the rest of the domain.
The soil at the site has been described as a typical hydromorphic 'boulbène' with a silt loam texture for the 1 m depth surface layer (14% sand and 28% clay). However, the clay fraction increased from 17% at the surface to 40% at 1 m depth. As shown in Fig. 3, five layers of soils were determined for the first 1.3 m by Enrique et al. [46] and Calvet et al. [47].  In the simulation, the initial ground temperature was approximated using an analytical expression proposed by Hillel [48], which provides the ground temperature profile based on soil parameters and climatic conditions, given as: where T(y, t) is the soil temperature at time t (d) and depth y beneath the ground surface (m), T a is the constant ground temperature ( • C), A 0 is the annual amplitude of the surface soil temperature ( • C), t 0 is the lag time from arbitrary start date to the occurrence of the minimum soil temperature in a year (d), and d is a damping depth, which is calculated as: where D h is thermal diffusivity (m 2 /d), and ω is 2π/365 (d − 1 ).
The parameters to determine the initial ground temperature are listed in Table 1. Although the ground temperature of the site in Southern France is not available in the literature [46,47], according to the on-site measurement elsewhere and literature data [26,49], the seasonal fluctuations in ground temperature decrease with the increasing depth. In addition, the lower the latitude, the higher the ground temperature at the same depth. Taking the UK as an example, mean annual temperatures at 1 m depth ranged from 12.7 • C in southern England to 8.8 • C in northern Scotland [49]. At the Wallingford weather station in the UK, a mean undisturbed ground temperature at 4 m beneath the ground surface was about 12 • C, with slight seasonal fluctuations [49]. Therefore, a fixed temperature with the reasonable value of 12 • C (285.15 K) was prescribed to the bottom of the domain for the current case.
Due to the lack of information on the initial pore-water pressure in the literature [46,47], several scenarios, as shown in Table 2, are made to calculate the initial pore-water pressure based on the value of h g and S l using the Van Genuchten [27] Model. Moreover, a water table was reported to be at a depth of 4 m below the ground surface during the inspection period [46], and therefore, the pore-water pressure was fixed as 0 Pa at the bottom of the domain.
The climatic variables of ambient air temperature, shortwave solar radiation, wind speed, air relative humidity, and rainfall were prescribed uniformly on the top of the domain (ground surface), and their variations can refer to Enrique et al. [46].

Material parameters
Each layer of soil was assumed to be isotropic in the simulation. The soil properties involving the hydraulic conductivity and SWCC have been summarized in Table 3 based on the measurements conducted by Enrique et al. [46] and Calvet et al. [47]. It should be mentioned that there were no material data available beyond the depth of 1.3 m, hence the material parameter for soils from the depth of 1.3 m-4.0 m was assumed to be the same as the fifth layer of soil. In addition, other soil parameters, such as the solid's density, specific heat capacity and thermal conductivity, as well as surface boundary-related parameters were gathered from the literature, as listed in Table 4.

Validation results
As shown in Table 2, three scenarios were adopted to determine the initial pore-water pressure for each layer of soil. The simulated daily evaporation and the measured data collected at the monitoring site are compared in Fig. 4.
From the figure, it can be seen that the maximum absolute error between the Base case and the measurement was 0.5 mm/day for the 65day simulation period. Although there are variations between the simulated Base case results and the measured data, which may be caused Table 1 Parameters to determine the initial ground temperature [26,48].
by the initial conditions assumed and the material parameters used, the modelling presented herein can capture the changing pattern of the daily evaporation. Furthermore, the daily evaporation simulated by the Base case, Case 1, and Case 2 shows no significant differences. Therefore, the daily evaporation calculation is not sensitive to the assumption of the initial pore-water pressure in the soils. Based on the aforementioned validation, it can be concluded that the surface ground boundary considering a range of climatic variables and mechanisms has been correctly implemented in the coupled TH model.

Application
The coupled TH model with the ground surface boundary has been applied to investigate the long-term (5 years) ground behaviour in the campus of the University of Warwick in response to seasonal climatic conditions. The findings of the test site will help to guide the utilization of the shallow ground source heat energy and the installation of the ground source heat system.   campus in Warwickshire County, UK. The site is located within the case study area of the EPSRC funded "Integrated Heating and Cooling Networks with Heat-sharing-enabled Smart Prosumers" project (Grant number: EP/T022795/1). The campus is a 24/7 town of 25,000 people, with 7000 student rooms, more than 150 buildings, 3 conference centres, 2 sport centres, retails, cafes, restaurants, offices and teaching buildings, industrial and research buildings. Peak electricity and heat demand in the campus are 10 MW e and 15MW th , respectively. The heat demand in the campus is primarily met by a district heating network supplied by natural gas CHP (5 CHP units with the total capacity of 8.2 MW e ) and auxiliary natural gas boilers (3 boilers with the total capacity of 14.7 MW th ). The heat network is comprised of 8.87 km supply pipes and 8.87 km return pipes and operates at 100 • C supply temperature and 50 • C return temperature. The University of Warwick's Estate Office is investigating options to decarbonise the heat system in the campus. One of the promising alternatives to the current natural gas CHP and boilers is large scale ground source heat pump.

Test site at Warwickshire, UK
The highlighted square with a dimension of 200 m × 200 m is the proposed test site investigated within the scope of this work. The soil profile layers are obtained from three borehole logs (as shown by red dots in Fig. 5) obtained from the British Geological Survey (BGS) [51] borehole database. The soil profile from the ground surface is categorized into three layers, that is, Layer 1: 0-0.3 m sandy clay loam, Layer 2: 0.3-2.4 m silty clay, and Layer 3: >2.4 m mudstone.

Initial and boundary conditions
In this study, a representative unit with three soil layers was employed. It was assumed that each soil layer was uniform. In addition, a uniform ground surface was assumed owing to the intricacy of the ground surface types in the test site. The 3D domain for the representative unit is shown in Fig. 6. Based on the soil profile and the computation efficiency, the domain depth was set to 4 m, and the length and width were taken as 2 m.
A pre-modelling exercise was conducted prior to the formal 5-year simulations. The pre-modelling lasted for a year, and the initial ground temperature was approximated using the analytical expression proposed by Hillel [48], and the initial pore-water pressure for each soil layer was computed using the Van Genuchten [28] Model, assuming the initial saturation was equal to 0.75.
In order to obtain the optimal mesh, various mesh sizes were used, which was reflected by the number of elements, and the related ground   temperature profiles at the end of the pre-modelling exercise were compared, as shown in Fig. 7. From the figure, it can be seen that there is a negligible difference in the ground temperature profiles as the element number of the mesh increases. Considering the computational efficiency and result accuracy, the discretized mesh consisted of 1125 hexahedron elements connected by 1656 nodes, as shown in Fig. 6, was adopted in the following simulations. Correspondingly, the ground temperature and pore-water pressure on the last day of the year became the initial ground temperature and pore-water pressure for the long-term simulations, as shown in Fig. 8.
As explained in the model validation, the mean undisturbed ground temperature at 4 m below the ground surface was about 12 • C [49], so a fixed temperature boundary (285.15 K) was prescribed at the bottom of the domain to consider the undisturbed ground temperature, and the pore-water pressure was assumed to be fixed at saturation of 0.75 at the domain bottom.
On the ground surface, climatic boundary conditions were evenly prescribed. The climatic data was obtained from a meteorological station near the University of Warwick (Church Lawford, Met Office [52]), and the most recent data of 2019 was used in the current simulations. Fig. 9 depicts the monthly average climatic variables.

Material parameters
Soil properties should ideally be measured using soil samples taken from the test site. In the current simulation, the soil properties for the three soil layers, as listed in Table 5, were derived from the literature [53][54][55] and the UNSODA Unsaturated Soil Hydraulic Database [56]. In addition, the values of model parameters in Table 4 were utilised and C y was given a typical value of 0.016 [26].

Simulated scenarios
A Base case was established using the aforementioned borehole logs, material parameters, and modelling procedure. Four more cases, as listed in Table 6, are designed to study the influences of the soil types and hydraulic drainage conditions on the ground temperature and saturation, and then the heat content of ground. Different soil profiles were employed in Case 1 and Case 2 compared to the Base case. Layers 1 and 2 in Case 1 are silty clay, whereas Layers 1 and 2 in Case 2 are sandy clay loam. In contrast with the Base case, the saturated hydraulic conductivities of Layer 2 in Case 3 and Layer 1 in Case 4 were increased by 1000 times.

Results of base case
Upon the collective inspection of the simulated results, both the Table 5 Material parameters for the three soil layers [53][54][55][56][57].    thermal and hydraulic ground behaviour reached a steady annual cycle state after approximately two years. In order to improve the clarity of plots, as shown in Fig. 6, the results of selected inspected points, including IPs 2-4, IP 11, IP 21, IP 25, IP 31, and IP41, during the 5-year simulation are illustrated. The observed cyclic behaviour was expected based upon the annually prescribed climatic data [26,58].
(1) Temperature and saturation As demonstrated in Fig. 10, the amplitude of the simulated temperature decreased with the increase of the depth, which is consistent with previous findings [11,12,20,26]. In comparison with other inspected points, significant changes in temperature of IPs 2-4 can be easily observed, indicating that the temperature of Layer 1 was most affected by the climatic conditions. The highest and lowest temperatures of IP 2 were 291.6 K and 275.1 K, respectively, then the amplitude of IP 2 was the greatest, reaching 16.5 K. Furthermore, curve peaks were off-set one another with respect to depth, with the deeper the inspected point, the later the time, as found in the literature [11,12,20,26].
As shown in Fig. 11, the simulated saturation of IP 2 and IP 3 almost coincided, ranging between 0.58 and 0.39 annually, with the peak value occurring in mid-February and the lowest value occurring in mid-August. Insignificant differences existed in the saturation of IP 11 and IP 21, both of which were in Layer 2. At IP 25, the interface between Layer 2 and Layer 3, the difference between the highest and lowest saturations in a year (0.37 and 0.13, respectively) was the biggest (0.24). The simulated saturation of IP 31 in Layer 3 varied slightly between 0.75 and 0.72 in a year.
(2) Heat content The development of the total heat content per unit area for the ground Λ total with time is shown in Fig. 12. By observing the figures, values of the total heat content per unit area for the whole domain and each soil layer repeated annually after approximately 2 years.
There are four seasons in the UK and in the Western Europe: Spring is from March to May, Summer is from June to August, Autumn is from September to November, and the remaining three months are Winter. Due to the periodic dynamics as shown in Fig. 12, the results in Year 3 (Day 730 to Day 1095) were studied in great detail and illustrated in Fig. 13.
As shown in Fig. 13(a), the value of Λ total for the whole domain varied from 9.05 × 10 7 J/m 2 to 1.15 × 10 8 J/m 2 , with the maximum value on Day 973, the final day of Summer, and the lowest value on Day 826, in the middle of Spring. Compared with the initial value on Day 730 (9.92 × 10 7 J/m 2 ), the Λ total generated from 1 January to 2 June in Year 3 (Day 730 to Day 883) was lower, indicating that the whole domain had been losing heat during this period. It was not until 3 June (Day 884) when a higher Λ total than the initial value on Day 730 was achieved, implying that heat was accumulated in the ground, and this trend continued until the last day of Year 3.
The total heat content per unit area for Layer 1 in Year 3 is shown in Fig. 13(b). The lowest value was 2.32 × 10 6 J/m 2 , which occurred on 31 January (Day 761), and the highest value of Λ total was 1.16 × 10 7 J/m 2 , which took place on 31 July (Day 942). In comparison with the initial total heat content per unit area on Day 730 which was 3.69 × 10 6 J/m 2 , a lower value of Λ total was only produced in the first 33 days in Year 3 (primarily in January), which indicated that Layer 1 was storing heat for about 91% of the year, mainly in Spring, Summer, and Autumn. Fig. 13(c) presents the total heat content per unit area for Layer 2 in Year 3. As can be seen from the figure, the highest and the lowest value of Λ total were 5.05 × 10 7 J/m 2 and 3.40 × 10 7 J/m 2 , respectively, on 2 September and 1 May (Day 975 and Day 851), respectively. In the first 170 days of Year 3, a smaller Λ total than the initial value on Day 730 (4.20 × 10 7 J/m 2 ) was produced, implying that Layer 2 had been losing heat for almost half of the year.
The total heat content per unit area for Layer 3 in Year 3 is given in Fig. 13(d). The value of Λ total in this layer changed from 5.07 × 10 7 J/m 2 to 5.50 × 10 7 J/m 2 . In Year 3, Layer 3 kept losing heat until 22 August (Day 964) when a higher Λ total than the initial value on Day 730 (5.34 × 10 7 J/m 2 ) was generated, which indicated that heat began to accumulate for the remaining time.  On the basis of the total heat content per unit area in Year 3 (Fig. 13) and Equation (38), the annual net heat content per unit area N HC for the whole domain and each soil layer can be calculated and summarized in Table A1 of the SI. The annual net heat content per unit area for the whole domain was positive and reached 9.52 × 10 8 J/m 2 , which indicated that complementary heat derived from solar radiation was stored in the ground. The 0.3 m-thick Layer 1 contributed the most, which was as high as 1.09 × 10 9 J/m 2 , while the 2.1 m-thick Layer 2 only provided 1.26 × 10 8 J/m 2 . Compared with N HC > 0 for Layers 1 and 2, the value of N HC for Layer 3 was negative with a value of − 2.61 × 10 8 J/m 2 , showing that heat was lost via Layer 3 during the year.
The seasonal and monthly net heat content per unit area N HC for the whole domain and each soil layer were computed and provided in Table A1 of the SI using a similar method. Fig. 14 illustrated changes in the seasonal and monthly net heat content per unit area for the whole domain and each soil layer.
According to the data and figure, for the whole domain, N HC > 0 can be obtained in Summer and Autumn. The value of N HC in Summer was 9.30 × 10 8 J/m 2 , which was very close to the annual net value (9.52 × 10 8 J/m 2 ), and most of it came from August with the value of 4.47 × 10 8 J/m 2 .
For Layer 1, N HC > 0 can be produced in Spring, Summer, and Autumn, and the value of N HC in Summer was the largest, reaching 6.45 × 10 8 J/m 2 , which accounted for 59.2% of the whole year (1.09 × 10 9 J/ m 2 ). Among the twelve months, the value of N HC in July was the highest and it was 2.37 × 10 8 J/m 2 . It should be noted that the current simulation's monthly N HC for Layer 1 was consistent with heat flux values for shallow depth soil in the literature [11,12,16,59].
For Layer 2, it can be observed from the table and figure that N HC > 0 took place in Summer and Autumn, and Autumn had a larger positive N HC equalled to 5.44 × 10 8 J/m 2 , much higher than the annul value of 1.26 × 10 8 J/m 2 . Throughout the year, N HC in September was the highest, reaching 2.40 × 10 8 J/m 2 .
Although the annual value of N HC was negative for Layer 3, N HC > 0 can still be found in Autumn with a value of 1.08 × 10 8 J/m 2 . The value of N HC in November was the greatest in the year and it was 4.32 × 10 7 J/ m 2 .
Based on the net heat content per unit area listed in Table A1 of the SI and considering the 200 m × 200 m inspected site in the University of Warwick, the annual/seasonal/monthly net heat storage N HS for the whole domain and each soil layer can be estimated following Equation (39). A large amount of heat from the solar radiation can be stored in the ground. For instance, the annul net heat storage for the whole domain reached 3.81 × 10 13 J, which can be potentially used as shallow ground source heat for building heating as environment-friendly energy. Among which, the annual net heat storage for Layer 1, Layer 2, and Layer 3 was 4.35 × 10 13 J, 5.06 × 10 12 J, and − 1.04 × 10 13 J, respectively.

Influence of soil typesresult comparison of base case, case 1, and case 2
The ability of the soil layer to store and conduct heat is related to the type of soil [11,16,18]. In addition to the soil profile in the Base case, Case 1 and Case 2 with different soil profiles were carried out (see Table 6), and their results were compared.
(1) Temperature and saturation Variations of the temperature of inspected points in Year 3 of the Base case, Case 1, and Case 2 are illustrated in Fig. 15. Although various soil profiles were adopted in the Base case, Case 1, and Case 2, the temperature changes for IPs in Layer 1 and Layer 2 were not significantly different, as shown in Fig. 15(a) and (b), respectively. Fig. 16(a) and Fig. 16(b) depict the saturation variations for inspected points in Layer 1 and Layer 2, respectively. The variations in soil types are reflected in the disparities in Ips' saturation. Saturation was lower for Ips in Layers 1 and 2 in Case 1 than that for Ips at the same depth in the Base case. In contrast with the Base case, there was a generally higher saturation for Ips in Layer 1 in Case 2, especially during Winter. Meanwhile, Case 2 exhibited greater saturation than the Base case for Ips at the same depth in Layer 2.
(2) Heat content Fig. 17 shows the total heat content per unit area for the whole domain and each soil layer obtained by the Base case, Case 1, and Case 2 in Year 3. According to Fig. 17(a), which presents the value of Λ total for the whole domain, the development trends of Λ total in the three cases were comparable, while from the perspective of the value of Λ total for the whole domain, Case 2 > Base case > Case 1, as a result of the differences in Λ total for each soil layer. Inspecting the variations of Λ total for Layer 1, Layer 2, and Layer 3, as shown in Fig. 17(b), (c), and Fig. 17(d), respectively, it can be found that: (1) the values of Λ total for Layer 1 in the Base case and Case 2 were close, and both were higher than that in Case 1; (2) the values of Λ total for Layer 2 in the Base case and Case 1 were similar, and both were much lower than that in Case 2; and (3) there were small differences in the values of Λ total for Layer 3 in the three cases. The initial value, maximum value, and minimum value of the total heat content per unit area in Year 3 of the Base case, Case 1, and Case 2 are compared in Fig. 18.
The values of the monthly net heat content per unit area N HC in Year 3 of the Base case, Case 1, and Case 2 are compared in Fig. 19. The values of the net heat content per unit area in Year 3 of Case 1 and Case 2 are provided in Table A2 and Table A3 of   Superimposing the monthly net heat content per unit area shown in Fig. 19, the annual/seasonal net heat content per unit area for the whole domain and each soil layer in Case 1 and Case 2 can be obtained. Compared with the Base case, the annual N HC equals to 9.52 × 10 8 J/m 2 , the values of annual N HC in Case 1 and Case 2 were 5.95 × 10 8 J/m 2 and 9.43 × 10 8 J/m 2 , respectively.
Based on the aforementioned analyses, it can be concluded that the spatial and temporal distributions of heat in soil layers of the ground can be significantly affected by the combination of the soil types (soil profile). It highlights the possibility of maximizing shallow ground source heat resources by improving soil types.

Influence of hydraulic drainage conditionsresult comparison of base case, case 3, and case 4
The saturated hydraulic conductivity of the soil would affect the drainage conditions in the ground. Higher hydraulic conductivity results into faster drainage and vice versa. Using the same model parameters as in the Base case but specifying a higher saturated hydraulic conductivity for Layer 2 in Case 3 and Layer 1 in Case 4 (see Table 6), the results of the Base case, Case 3, and Case 4 were compared.
(1) Temperature and saturation In Year 3, the saturation variations of inspected points in the Base case, Case 3, and Case 4 are illustrated in Fig. 21. It can be found that compared with the Base case, the saturation distributions in Case 3 were slightly changed, while Case 4 was basically unchanged.
(2) Heat content The total heat content per unit area for the whole domain and Layers 1-3 in Year 3 of the Base case, Case 3, and Case 4 is presented by Fig. 22 (a), Fig. 22 Table A4 and Table A5 of the SI, respectively. According to Fig. 23, except the monthly N HC for the whole domain (0-4 m) and Layer 2 (0.3-2.4 m) in Case 3 was slightly lower than that in the Base case, insignificant differences can be found in the three cases.
Based on the above analysis, while keeping other model parameters unchanged, when the saturated hydraulic conductivity of the first and second layers increased by three orders of magnitude in the simulations, changes in the temporal and spatial distributions of the shallow ground heat were not significant.

Conclusions
In this paper, a framework to calculate and to predict theoretical heat capacity, temperature, and soil-moisture behaviour of a ground for potential exploitation of shallow ground source heat (SGSH) is presented. Atmospheric data e.g., solar radiation, rainfall, humidity, air temperature, wind velocity is considered in the study, to estimate temperature and ground heat content of individual subsurface soil layers. A good agreement between the model predicted results and the in-situ experimental data is observed from the model validation exercise.
The model is applied to investigate ground temperature, heat capacity, and moisture distribution of a site in Warwickshire County, UK for a period of 5 years. The SGSH reserves of the test site, such as the total heat content per unit area Λ total and the annual/seasonal/monthly net heat content per unit area N HC are assessed. From the base case  scenario, it was found that the Λ total of the 4-m-deep domain varied from 9.05 × 10 7 J/m 2 to 1.15 × 10 8 J/m 2 and the annual N HC of the 4-m-deep domain reached 9.52 × 10 8 J/m 2 . Therefore, for the test site with an area of 200 m × 200 m, the annual net heat storage N HS can be as high as 3.81 × 10 13 J. This large amount of heat can be potentially used as ground source heat for building heating as environment-friendly energy.
Several scenarios where the influences of various soil layer profiles and hydraulic drainage conditions on the ground temperature, saturation, and heat content were considered, it was found that the soil types and conditions had greater impact on ground thermal potential than hydraulic drainage conditions. Improvement of soil types, in this context, could perhaps lead to enhanced ground heat capacity. For example, horizontal ground heat exchangers are not recommended in sandy soil as heat imbalance is prone to occur due to low specific heat capacity and high thermal conductivity of sand. Mixing or replacing sand with other soils with a higher specific heat capacity (e.g., clay or mudstone) can increase its heat storage potential. In addition, introducing phase change materials into the ground can also enhance its inter-seasonal heat storage capacity.
An in-depth understanding of the thermo-hydraulic behaviour of the ground is critical for realistically quantifying the annual renewable heat supply and the potential for inter-seasonal heat storage. Therefore, accurate estimation of the ground behaviour can significantly affect the optimal design of potential heating and cooling networks and their technical and economic feasibility. In addition to the thermo-hydraulic behaviour of the ground, different types of the ground loop i.e., shallow horizontal heat collector vs. deep vertical borehole, can present significantly different characteristics. The shallow ground source heat resources of a specific site analysed using the assessment method proposed in this study can be used as a basis for the feasibility of a horizontal ground source heat pump system. The horizontal ground source heat pump system is not suitable if the shallow ground displays a negative net heat content.
Furthermore, the ground condition could have a significant impact on the CoP of a ground source heat pump supplying a district heating network, or on the supply temperature of an ambient loop heating and cooling network (i.e., 5th generation heating and cooling networks). This, consequently, affects the optimal design (e.g., size of pipes and circulation pumps) and operation of the heating and cooling networks and their life-cycle cost.
Therefore, incorporating the thermo-hydraulic behaviour of the ground in the design and operation modelling of heating and cooling networks is essential to ensure best technical and economic performance will be achieved.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.