A simplified method to evaluate geothermal storage in an aquifer with consideration of heat transfer between aquifer and caprock/baserock

. Storing and extracting heat during different seasons of the year is possible through the utilization of a ground aquifer with an open loop Ground Source Heat Pump (GSHP) system. Being able to predict the hydrothermal performance of geothermal storage is required for an efficient operation of the system for cooling and heating of buildings. Complex 2D and 3D hydrothermal numerical models can simulate the thermal performance of geothermal storage accurately but often lack the desired computational speed for conducting large number of simulations for performance optimization. Instead, a 1D radial model can be used to conduct fast evaluation. However, it is important that the model computes the amount of heat loss from an aquifer into the overburden and underlying layers accurately to evaluate the amount of geothermal storage in the aquifer at different times. In this study, a source term is introduced into a 1D model to simulate the heat transfer between the aquifer and caprock/baserock in the vertical direction. The following two heat loss models are introduced in the heat advection-conduction equation: (i) Newton’s heating/cooling law, which leads to a closed form solution, and (ii) a conduction-based semi-analytical model, which requires a 1D finite element solution. When compared to a full 2D axisymmetric simulation result, it was found that the Newton’s heating/cooling law model with a constant heat transfer coefficient works well in cases of fast heat flow rate in thick aquifers of around 100 meters. But large errors in estimating heat dissipation are observed in cases with low heat flow rate in thin aquifers, especially for simulations exceeding two to five years. On the other hand, the model with the conduction-based semi-analytical solution gives a better match for these conditions.


Introduction
A centralized ground source heat pump (GSHP)-based heating and cooling system that is directly interacting with the groundwater through a geothermal borehole is an open-loop system, while a closed-loop system consists of heat pumps with ground heat exchangers [1]. In general, during summer, the excess heat generated inside buildings by cooling is dumped into the ground and in winter that stored excess heat is utilized for heating through GSHPs. However, in reality, GSHPs operate in a cyclic manner depending on the heating and cooling demands that vary daily, weekly and monthly. Because of the variation in the actual heating and cooling demands, one of the challenges in the GSHP industry is to accurately predict the thermal performance of the overall system.
Complex hydrothermal 2D and 3D numerical models are able to simulate the thermal performance of geothermal storage accurately but often lack the desired computational speed for conducting large number of simulations for performance optimization. Instead, a 1D radial model can be used to conduct fast evaluation. For a closed loop system, a 1D conduction model that incorporates the radial heat transfer is available [2] (e.g. g-function by Claesson and Javed, 2012). For an open loop system, the performance of geothermal storage is governed mainly by the loss of heat from the aquifer layer into the overburden and underlying layers. In this study, a sink term is introduced into a 1D heat advectionconduction model to simulate the heat transfer between the aquifer and caprock/baserock. The following two heat loss models are introduced in the governing heat transfer equation: (i) Newton's heating/cooling law, which leads to a closed form solution, and (ii) the conduction-based Vinsome and Westerveld [3] model, which leads to a 1D finite element solution.

Solution methods
The governing equation for the heat advection and conduction in radial coordinates is given as: where θ is the temperature within the aquifer, D is the thermal diffusivity, H is the thickness of the aquifer, ρ is the density of the aquifer, c is the specific heat capacity of the aquifer, Q is the fluid flow rate, r is radius, w0 is the heat flux at the interface between the aquifer and overburden layer. In this study, two models are implemented for w0.

Newton's heating/cooling law based analytical solution
A Newton's heating/cooling law-based analytical solution was developed by Lu [4]. The boundary heat flux term w0 is equal to h×(θ-Tc) in which h is the convective heat transfer coefficient and Tc is the ambient temperature of the caprock/bedrock. The source term for the heat loss calculations contains an empirical expression for the heat transfer coefficient, which was derived through compiling various deep geothermal system simulations with realistic flow rates and heat conduction properties (Lu, 2019). The aquifer was assumed to be deep and thick (more than 100 meters) and the flow velocity was relatively large (between 0.001 to 0.01 m 3 /sec per meter depth). The interface convective heat transfer considers both the heat transfer between the aquifer and the caprock. The dimensionless analytical solution T(r,t) (r=radius, t=time) for Eq. (1) based on Newton's heating/cooling law is shown by Eq. (2). where T0 is the initial aquifer temperature, Tw is the injected water temperature, n is the porosity, R is the delay factor, λ is the thermal conductivity of the aquifer, ρ is the density of the aquifer, c is the specific heat capacity of the aquifer, s is the solid, f is the fluid, H is the thickness of the aquifer, and h * is the dimensionless heat transfer coefficient. Lu (2019) performed a series of 2D axisymmetric finite element simulations of full-scale overburdenaquifer models to relate the heat transfer coefficient h to a given set of heat injection rate and ground thermal properties. Results show that the dimensionless heat transfer coefficient h* depends on the terms of dimensionless heat injection rate and dimensionless thermal conductivity and the empirical relationships obtained through curve fitting of the finite element simulations are given by the following equations.
To incorporate the interaction of injection and extraction wells (i.e., a doublet well system), Lu (2019) modified the dimensionless injection rate α* as follows. * = + 0.0914 . . * , 0.6 ≤ ≤ 600 (5) where L* is the dimensionless distance, which is defined as L/H and L is the distance between the injection and extraction wells. Figure 1 shows the performance of the 1D analytical model when compared to a 2D axisymmetric FE model with the dimensions of a rectangle of 200 meters in length and 50 meters in aquifer thickness (actual aquifer thickness of 100 m is halved considering axisymmetric conditions). In these simulations, the initial temperature of the aquifer is assumed as 14.8°C with the hot-water injection with a temperature of 25.8°C is performed. The purpose of those simulations is to understand how the temperature front moves with time under different injection rates. The thermal properties used in the COMSOL model are provided in Table 1.
The upper and lower boundaries were set as convective heat transfer boundaries with a heat transfer coefficient being equal to h=5 W/(m 2 K). Simulations of two relatively fast injection rates (0.01m 3 /sec and 0.001 m 3 /sec per meter depth) are performed up to 1.6 years and the two computed temperature variations inside the aquifer match well.  When the aquifer thickness is increased to 100 m (Fig. 2b), the match is good after 1 year of operation but not for the 10 and 40-year operations. Close examination of the results show that the error originates from large heat dissipation by slow conduction from the aquifer to the overburden at larger times. Hence the Newton's Heating/Cooling Law Based Analytical Solution, developed for a deep geothermal aquifer system, may not be applicable for a shallow aquifer system because the heat dissipation is governed by slow thermal conduction from the aquifer into the overlying/underlying layers.

Vinsome & Westerveld (1980) conductionbased numerical solution
The conduction-based semi-analytical model by Vinsome and Westerveld (1980) is often used for heat loss in thermal reservoir simulations. The semi-analytical model is implemented into a finite element method-based coupled model. In this case, the source term w0 of Eq. (1) utilizes the following expression derived by Vinsome and Westerveld (1980) where θ is the interface temperature, which is assumed to be temperature T in Eq. (1), d is the diffusion length, and λ2 is the thermal conductivity of the caprock. p is the fitting parameter for the following expression which is assumed to be distribution of the temperature in caprock. The fitting parameter p is formulated by Vinsome and Westerveld (1980) after they have chosen a twoparameter flexible function that would fit the diffusivity equation and satisfy the initial and boundary conditions of the problem.

(7)
The diffusion length d is calculated using the following equation: where κ is the thermal diffusivity that is related to the thermal conductivity and capacity of the caprock, while t is the time calculated from the beginning of the simulation. The thermal diffusivity can be expressed in terms of caprock parameters as follows: where λcap and ρcap are the thermal conductivity and heat capacity of the caprock, respectively.
Considering the heat flow at the interface between the aquifer and caprock with the application of the boundary conditions results in the following expressions for the curve fitting parameters I N , p N and q N , where N is the superscript denoting the value from the old time step. The explicit forms of the curve-fitting parameters are given in the following equations. The temperature calculation for the aquifer and caprock-aquifer interface can be explained in a sequential way. At the beginning of each time step, the diffusion length d is calculated using the total amount of time passed since the beginning of the simulation. After having the diffusion length, the curve-fitting parameters p and q are obtained using Eq. (11) and Eq. (12), respectively. The value for I N term for the next time step is calculated using Eq. (10) after having the curve-fitting parameters in hand.
In order to calculate the amount of heat loss, Eq. (6) is utilized in order to have the forcing vector for the finite element solution. By completing the whole calculation scheme, the temperature calculation problem for the current time step is done through the application of the finite element method. The computations give the temperature value within the aquifer for the current time step, and the same procedure explained above is used for the next time step.
In order to investigate the performance of this model for shallow geothermal storage problems, a benchmark problem consisting of a 100-meter thick caprock overlying a 20-meter thick aquifer is considered. These site conditions are adopted from a case study in Portland, Oregon where a feasibility study for a direct use geothermal application was conducted [5]. For the simulation, hot water is pumped at a rate of 0.015 m 3 /sec with a temperature of 25.8 degrees Celsius. The initial temperature in the aquifer is 14.8 degrees Celsius. A 5year long operation is simulated. A 2D axisymmetric numerical model that includes both aquifer and overburden layers developed using COMSOL was used for the comparison.
The temperature profiles in the middle of the aquifer computed from the two models are shown in Fig. 3. In the early stage of hot water injection (year 1), the frontal location computed by the 1D FEM code is more than that of the 2D model. However, at the end of the tenth year, the shapes of the temperature fronts are close to each other. This is because the heat loss to the overburden layer is more conduction dominated as the frontal location moves away from the injection well. Hence the Vinsome & Westerveld (1980) based 1D numerical model can be a good alternative for initial estimation of the temperature profile calculations for long operational conditions.  Fig. 4 shows the aquifer temperature profiles with and without heat loss to the caprock from the aquifer at 5, 10 and 20-year operations. As expected, the temperature profiles are sharp in space when there is no heat loss to the overburden layer, but they are more spatially distributed when there is heat loss to the overburden layer. The front locations which show an excessive temperature change in a small distance are similar after 5 and 10 years of operation, but they become different after 20 years of operation. Results show large heat loss when the heat is injected for long time in a thin aquifer typical utilized by an open loop GSHP system.

Conclusions
When heat transfer in an aquifer has large advection component (like near the injection well), the heat transfer from the aquifer to the overburden/underlying layers is dominated by the advection characteristics and results show that the Newton's Heating/Cooling Law based analytical solution is applicable for such conditions (thick aquifer, fast injection rate and short operation time). The advantage of this solution is that it is a closed form solution and the computational speed to find the solution is very fast.
When heat transfer in an aquifer occurs with slow moving groundwater flow, the heat transfer mechanism from the aquifer to the overburden/underlying layers becomes conduction dominated. This condition occurs in a thin aquifer for a long operational time when the heated/cooled water moves far away from the well. In such conditions, the Vinsome & Westerveld (1980) based numerical solution appears to be applicable to evaluate the temperature profile inside the aquifer and the amount of heat stored in the aquifer at different times. Since it is a numerical-based model, the computational demand can be large when large number of simulations are required for system optimization purpose. However, the computationally efficient 1D FE model developed in this study may be used for such purpose.
Further work includes the study on the range of thermal properties and injection rate where either Newton's heating/cooling law or the conduction-based Vinsome and Westerveld (1980) approach work in different geothermal storage cases.