## Abstract

**Abstract** Previous research has derived a general 2D analytical solution, using stream functions, to the case of an abstraction–injection well doublet in a regional groundwater flow field. An easily applied analytical solution is presented for the special case where the injection well lies immediately down-gradient of the abstraction well and where the doublet provides ground source heating or cooling. The resulting system of equations can rapidly be used to find a first approximation for (1) the risk of internal thermal feedback in the doublet and (2) the time at which it will occur, (3) the future equilibrium abstraction temperature, (4) the width and (5) the length of the thermal plume migrating down-gradient from the reinjection well.

The use of the ground as a resource for heating and cooling, employing both open-loop (groundwater-based) and closed-loop heat exchangers, is expanding rapidly in many European countries (Banks 2009*a*). A common configuration for an ‘open-loop’ ground source heating or cooling scheme (Banks 2008) is an abstraction–reinjection well doublet. According to Banks (2009*b*), such a scheme typically comprises three elements: (1) an abstraction well, from which water is abstracted at a rate *Q* and a temperature θ_{gwabs}; (2) a heat-transfer system (a heat exchanger or a heat pump), which either extracts heat from, or rejects heat to, the groundwater flux; (3) one (or more) reinjection well(s), at some distance from the abstraction well, where water is reinjected at a rate *Q* and temperature θ_{gwinj}; for space-cooling schemes, θ_{gwinj} > θ_{gwabs} and for heating schemes θ_{gwinj} < θ_{gwabs}.

To minimize the risk of internal thermal feedback of reinjected water, the reinjection well is normally located directly down the hydraulic gradient from the abstraction well. It is in the interests of the scheme's operators to minimize the risk of thermal feedback. This implies, however, that they are maximizing the amount of heat being ‘lost’ from the doublet in a thermal plume that migrates down-gradient with regional groundwater flow.

In many nations, including the UK, the environmental regulatory authorities are interested in being able to predict the size and speed of this thermal plume (Todd & Banks 2009; CZC 2010), to ensure that down-gradient users of the aquifer are not adversely affected by changing groundwater temperature. This paper proposes an ‘introductory’ tier of risk assessment that can be applied rapidly by both regulators and scheme operators to predict the extent of thermal plumes (subject to a number of key assumptions) and the extent and magnitude of thermal feedback.

## Internal risk of thermal feedback

If we assume that the doublet spacing is 2*d*, that the wells fully penetrate the aquifer and that the aquifer is homogeneous with a uniform thickness *m*, it can be demonstrated (by balancing two times the Darcy velocity owing to pumping (*Q*/2π*dm*) against the regional Darcy velocity *U*) that there is minimal risk of internal feedback (Banks 2009*b*; Fig. 1) if (1) where *U* is the regional Darcy velocity (specific groundwater flux). The thermal ‘plume’ emanating from the reinjection well is separate from the capture zone of the abstraction well and two stagnation points are located on the axis between the two wells (i.e. on the *x*-axis in Fig. 1). The maximum down-gradient width of the thermal plume (and also the width of the up-gradient capture zone; see, e.g. Javandel & Tsang, 1986) *W*_{pl} is simply given by (2)

If the well separation is less than the critical distance given by equation (1), then there is potential for thermal feedback (Fig. 2) of water from the reinjection well to the abstraction well, potentially compromising the long-term viability of the system. Banks (2009*b*) has considered the internal risk of thermal feedback in this scenario and concluded that the hydraulic (*t*_{hyd}) and thermal (*t*_{the}) travel times from the recharge well to the abstraction well along the shortest flow path are given by the formula of Lippmann & Tsang (1980) as (3) (4) where instantaneous thermal equilibrium between water and aquifer matrix are assumed. *S*_{VCaq} and *S*_{VCwat} are the volumetric heat capacities of the bulk saturated aquifer and of water, respectively, *n*_{e} is effective porosity and (5) The dependence of thermal breakthrough time on β is shown in Figure 3: as β increases, the doublet spacing (2*d*) decreases and the breakthrough time (*t*_{the}) sharply decreases.

## External risk: dimensions of a thermal plume

Banks (2009*b*) also noted the potential ‘external’ risk to other aquifer users located down-gradient of the doublet, in the path of the thermal plume emanating from the recharge well, but did not consider the shape of the thermal plume in detail. It is the shape and velocity of this thermal plume that this paper will consider, together with the proportion of the reinjected water ultimately recirculated back to the abstraction well.

Both Shan (1999) and Luo & Kitanidis (2004) have already provided a general 2D analytical solution to the problem of a well doublet within a regional flow field. Luo & Kitanidis (2004) considered cases where an abstraction well and injection well exist within a regional groundwater flow field, where the regional Darcy velocity vector is at any angle to the *x*-axis. Solutions were derived by considering discharge potentials (Φ) and stream functions (Ψ); these concepts have been fully discussed in several hydrogeological textbooks (e.g. Bear 1972; Raudkivi & Callander 1976). In the case where the well doublet is aligned along the *x*-axis (Figs 1 and 2), with the abstraction well at (*x* = -*d*, *y* = 0) and the injection well at (*x* = *d*, *y* = 0), Ψ and Φ can be obtained additively, as follows: (6) (7) where and are both in the range −π to + π radians and *U _{x}* and

*U*are the components of regional Darcy flux in the

_{y}*x*- and

*y*-directions.

Both Javandal & Tsang (1986) and Luo & Kitanidis (2004) have provided explicit solutions to these equations for various arrangements of wells, but not for the specific case of the thermal well doublet where the reinjection well is located directly down-gradient from the abstraction well. It is this special case that this paper addresses, where there is no component of regionalgroundwater flow in the *y*-direction (i.e. *U _{y}* = 0 and

*U*=

_{x}*U*). Thus, equation (7) reduces to (8) According to Luo & Kitanidis (2004), the location of the stagnation points (

*z*

_{s}in complex number space) of the flow system can be found by (9) where is the complex conjugate of the regional Darcy flow vector. However, in our special case, the regional flow is in a direction parallel to the

*x*-axis and equation (9) simply becomes (10) Thus, referring back to equation (1), if β < 1 , then the root has a real solution, yielding locations of the stagnation points (

*x*

_{s}) on the

*x*-axis (Fig. 1): (11) If β > 1 , then the root has an imaginary solution, yielding locations of the stagnation points (

*y*

_{s}) exactly on the

*y*-axis (Fig. 2): (12) where i = . Thus, the flow system has two stagnations points (0,

*y*

_{s}) and (0,

*−y*

_{s}) and these stagnation points must lie on the streamline that defines the boundary of the thermal plume. Thus, the bounding stream function Ψ

_{s}is given by (from equation (8)) (13) Luo & Kitanidis (2004) further demonstrated that the flow that escapes down-gradient in a thermal plume (

*Q*

_{pl}) is given by (14) and the recirculated (feedback) flow

*Q*

_{r}is given by (15) The dimensionless nomogram in Figure 3 shows the dependence of the ratios

*y*

_{s}/

*d*and

*Q*

_{pl}/

*Q*on the parameter β. Equations (14) and (15) allow us to calculate not only the asymptotic down-gradient width of the thermal plume (ignoring dispersion),

*W*

_{pl}, (16) but also the asymptotic temperature of the abstracted water (θ

_{abseq}) from the doublet (17) where θ

_{gwinj}is the constant reinjection temperature (in reality, of course, the reinjection temperature will probably increase in parallel with the increasing abstraction temperature, but this effect is ignored here). θ

_{o}is the ambient undisturbed groundwater temperature.

Finally, the Darcy velocity (*v*_{D}) profile along the *x*-axis (*y* = 0) is given by (18) This equation could be solved for small discrete intervals along the streamline to yield approximate travel times and travel distances from the well along the *x*-axis. The Darcy velocity typically returns relatively rapidly towards the value of *U* (Fig. 4) and thus the length of the thermal plume (*L*_{pl}) for long elapsed distances and times (*t*) can be approximated by Darcy's Law with an appropriate thermal retardation factor (Banks 2009*b*): (19) However, as Prof J. Barker points out (J. Barker, pers. comm., 2010), it is possible to derive an exact (although not uncomplicated) solution to the travel time from the well along the *x*-axis from equation (18).

## A worked example

As an example, a homogeneous sandstone aquifer will be considered, of 50 m thickness and transmissivity 250 m^{2} day^{−1}, implying a hydraulic conductivity of 5 m day^{−1}. Let us further assume a regional hydraulic gradient of 0.01, yielding a Darcy flux of 0.05 m day^{−1}. Our thermal well doublet is separated at a distance of 60 m and pumps water from the aquifer at an initial temperature of 10 °C and at an average rate of 10 l s^{−1}. Water is reinjected at a constant temperature of 20 °C. The volumetric heat capacity of the saturated aquifer (*S*_{VCaq}) is 2.2 MJ m^{−3} K^{−1} and that of water (*S*_{VCwat}) is 4.19 MJ m^{−3} K^{−1}. Thus, *d* = 30 m, *Q* = 864 m^{3} day^{−1}, *U* = 0.05 m day^{−1}, θ_{o} = 10 °C, θ_{gwinj} = 20 °C, *m* = 50 m, *Q*/*m* = 17.28 m^{2} day^{−1} and β = 3.667.

The stagnation points are given by equation (12): as (0, 49) and (0, −49), in metres. This allows the calculation of the bounding stream function from equation (13), and it is found to have values of 5.47 and −5.47 m^{2} day^{−1}. Thus, the recirculated quantity of water is given by equation (14) as and the volume of water (*Q*_{pl}) migrating down-gradient in the thermal plume is 547 m^{3} day^{−1}. This implies that: (1) the thermal plume would tend towards a width of *W*_{pl} = 219 m (equation (16)); (2) the thermal plume would have a length of *c*. 870 m after 25 years (9131 days); (3) the abstracted water would tend towards an equilibrium temperature of 13.7 °C. The velocity profile along the *x*-axis is shown in Figure 4.

Table 1 illustrates the effect of altering the doublet spacing. The critical distance 2*d* between the wells of the doublet is found (from equation (1)) to be 220 m. Thus, at a spacing of 250 m, no solution is found to the set of equations described above and the thermal plume separates from the capture zone of the abstraction well.

## Applicability to real hydrogeological scenarios

The reader's attention is drawn to the main assumptions that underlie this approach, namely:

constant pumping rate and injection temperature;

a homogeneous, isotropic aquifer of constant thickness, with Darcian flow;

confined aquifer (although the approach will apply to unconfined aquifers if head variations are small relative to total thickness);

fully penetrating wells;

instantaneous thermal equilibrium between groundwater and matrix (for equations (4) and (19)); this assumption may be poorly valid for fractured or dual-porosity aquifers;

injection well immediately down-gradient of abstraction well;

there is no vertical conduction (or other transfer) of heat to overlying or underlying strata (the aquifer is thermally confined); in other words, groundwater and heat fluxes are treated as wholly 2D;

the model considers an ‘equilibrium’ plume development;

dispersion is neglected.

In real well-doublet ground source heating and cooling scenarios, the pumping rate and injection temperature will usually vary diurnally and seasonally inresponse to heating and cooling demand. This does not hinder the utility of the system of equations presented, however, provided average pumping rates and temperatures over an appropriately long period can be calculated. Of the remaining assumptions, three will often be significantly violated in the context of real well-doublet scenarios, as follows.

### Three-dimensionality

In real hydrogeological scenarios, there will always be the possibility for heat to migrate via conductive processes into overlying and underlying strata. Moreover, if the aquifer is relatively close to the ground surface, the atmosphere can be approximated to a constant temperature (i.e. annual average soil or air temperature) source or sink for heat. 3D finite-element numerical modelling work suggests that the vertical conduction of heat away from the aquifer will always be highly significant for the heat balance around a well doublet, effectively further retarding the lateral migration of heat and restricting lateral thermal plume development by typical (CZC 2010) factors of at least 60–70%. It might thus reasonably be argued that this observation negates the simple 2D analytical approach described above, necessitating the use of complex 3D finite-element models for the realistic simulation of subsurface heat transport. The author argues, however, that the 2D analytical model is many times simpler and faster than such complex numerical models and is also highly conservative in predicting thermal plume evolution and breakthrough. Thus, if the use of the 2D analytical model at an early stage of risk assessment can demonstrate the lack of significant downstream thermal impacts and doublet breakthrough, one may be confident that a fully 3D numerical model will reach the same conclusion. If possible adverse impacts are predicted by the 2D analytical model, then a ‘tiered’ risk assessment may progress to more costly and time-consuming approaches to ascertain whether such impacts are indeed likely to occur.

### Instantaneous thermal equilibration

In porous, unlithified, clastic aquifers (e.g. sands, silts and gravels), the assumption of instantaneous thermal equilibrium is believed (de Marsily 1986) to be valid, as the groundwater is in intimate thermal contact with small-sized aquifer clasts. In aquifers where preferential flow pathways are dominant (e.g. fractured hard-rock aquifers, fissured karstified aquifers, or dual-porosity aquifers such as the British Chalk), the assumption of instantaneous thermal equilibrium becomes less tenable. This is because groundwater conduits of limited dimension and volume (e.g. Chalk fissures) are separated from each other by relatively large blocks of aquifer matrix. The time taken for heat to be conducted into the large matrix blocks from the groundwater in the fissures will be long and one cannot assume instantaneous thermal equilibration (i.e. the aquifer matrix and the groundwater will not be characterized by identical temperatures). This implies that there is a risk that the analytical model described in this paper will underestimate thermal travel velocities in dual-porosity aquifers, such as the Chalk, at least over short distances and time periods. Barker (2010) has proposed a somewhat more complex analytical approach for evaluating heat transport in such fractured or dual-porosity aquifers.

It can be argued that, over longer time periods (of years or decades) and distances (of hundreds of metres), the time for thermal equilibrium between fractures and matrix becomes short, relative to overall travel time, and thus that standard analytical models (such as described in this paper) and numerical models (which also typically assume instantaneous thermal equilibrium) become increasingly valid. Over short periods or during short-term testing, especially in dual-porosity or fractured aquifers, the applicability of the assumption of instantaneous thermal equilibrium needs to be carefully considered.

However, Todd & Banks (2009) have described a well doublet in the Sherwood Sandstone in Selby, Yorkshire, where simple analytical approaches give conservative prognoses, despite the fact that the Sherwood Sandstone is known to display a dual-porosity nature. Clarkson *et al*. (2009) have described a well doublet in the London Chalk, where very rapid (4.5 h) breakthrough of a conservative chemical tracer was proven during a 29 day constant rate test, but where no thermal breakthrough was detected during the entire test. Gropius (2010) has described another example from the London Chalk where an equivalent porous medium numerical model was successfully utilized to simulate thermal breakthrough in a well doublet, without explicit consideration of dual-porosity or fracture pathways.

### Equilibrium plume form

Finally, it should be noted that the systems of equations described in this paper consider the equilibrium situation that a plume will evolve towards. In cases where the doublet pumping rate is high compared with the regional Darcy velocity, the equilibrium plume may be very wide, with heat migrating rather slowly down-gradient with natural groundwater flow. It may therefore take many decades or centuries before sufficient heat has entered the aquifer to ‘fill’ the plume. In such cases, or in cases where the time-period being considered is very short, the equilibrium plume shape may not be a good representation of the actual plume evolution at finite times.

## Conclusion

A system of equations has been described in this paper and in an earlier paper (Banks 2009*b*) that allow a hydrogeologist or thermogeologist to construct a mathematical spreadsheet that provides a 2D analytical solution for an open-loop ground source heating or cooling well doublet in a regional groundwater flow field. It is thus possible to predict:

whether there is a risk of internal thermal feedback in the doublet (equation (1));

the time at which the thermal feedback will occur (equation (4));

the eventual equilibrium temperature of the abstracted water (equation (17));

the width of the thermal plume migrating down-gradient from the reinjection well (equation (16));

the length of the thermal plume migrating down-gradient from the reinjection well (equations (18) and (19)).

It is emphasized that this approach does not represent any new mathematical breakthrough. It is a special case of the general solution described by Luo & Kitanidis (2004), but one that is especially relevant to the rapidly growing ground source heating and cooling industry. It illustrates the great power of the stream function in complex number space to provide solutions to multiple well configurations. It should be noted that the system of equations described in this paper contains only three independent variables describing the hydraulic behaviour of the doublet system: the well discharge per unit thickness of aquifer *Q*/*m*, the regional Darcy velocity (specific groundwater flux) *U* and the well doublet spacing 2*d*. Porosity, transmissivity and hydraulic conductivity do not appear as explicit variables. The rate of migration of the plume along the *x*-axis is controlled by one additional factor (assuming instantaneous thermal equilibration between groundwater and aquifer matrix), namely, the ratio of volumetric heat capacity of the saturated aquifer to that of water.

The 2D approach described in this paper is affected by the questionable validity of at least two fundamental assumptions: (1) an equivalent porous-medium aquifer is assumed, with instantaneous thermal equilibrium between groundwater and matrix; (2) no 3D heat transfer is considered (i.e. no heat conduction into overlying and underlying strata). Nevertheless, it is argued that the simple 2D analytical approach can be regarded as conservative (at least when relatively long time-frames of years to decades are considered), and can thus be applied as a rapid, early stage of risk assessment of well doublet systems.

## Acknowledgements

The author wishes particularly to thank Professor J. Barker for his valuable and instructive criticism of this article and for providing additional solutions that he hopes he will turn into a discussion paper in due course.

- © 2011 Geological Society of London