## Abstract

**Abstract** The Cooling the Tube Programme was set up to implement London Underground's intention to mitigate future warming of the underground system and to control tunnel temperatures. As part of this programme the use of groundwater from the Chalk in open-loop systems was considered for cooling a number of stations and a groundwater model was determined as a requirement by the Environment Agency to assess impacts on the Chalk aquifer. It seemed plausible that the plumes of injected, warmed water would interact hydraulically and possibly also thermally with each other at some locations. Furthermore, London Underground wished to investigate any longer-term warming effects and potential loss in cooling benefit. This paper provides an overview of the hydrogeological modelling approach that was adopted to assess the potential hydraulic and thermal effect of the proposed schemes and summarizes some of the findings that may be of broader interest. A staged approach was used to guide and refine the modelling. Initially, analytical solutions were used to investigate the processes of heat transport in the Chalk and to derive parameters from a field-scale test. The results indicated that, even for low fissure-to-matrix contact areas, heat is likely to be conducted far into the Chalk matrix. Over longer time scales, significant heat loss occurs into the under- and overlying formations and the available analytical solutions are of less use in scoping heat calculations for open-loop schemes at this scale. Distributed, finite-element, numerical models were constructed to simulate the interactions within groups of ground source cooling systems and to investigate specific operating conditions at single schemes. Extensive sensitivity analysis was carried out to assess the level of uncertainty in the model predictions, and the relative sensitivity to the various parameters. The greatest uncertainty was associated with detailed aspects of the conceptual model of heat flow. Model results were sensitive to the assumed vertical distribution of permeability in the Chalk. Additional uncertainty was associated with the selection of appropriate thermal boundary conditions and the extent of interaction with the under- and overlying formations. Uncertainty in the model hydraulic and thermal material parameters had less effect on the simulated performance of the systems. Operational conditions, such as the amount of heat injected and the separation of the abstraction and injection boreholes, were assessed and this paper provides some quantification of the significance of these factors. This work also demonstrates that equivalent MT3D and FEFLOW models give reasonable agreement and both approaches are feasible for modelling heat transport, as density-driven flow is not significant for the temperature differences associated with these cooling schemes.

Summer platform temperatures in the London Underground (LU) are typically between 30 and 32 °C during peak hour operations, and are still a relatively warm 24 °C during the winter months. LU services are planned to increase to meet the predicted growth in passenger demand; new trains are on order and more trains, moving more customers, require more energy, which creates more heat. Transport for London (TfL) has identified the need to invest £150 million to address the issue of heat on the LU Victoria Line and restrict platform temperatures to a maximum of 29 °C. This is a huge engineering challenge and TfL established the Cooling the Tube Programme (CTP) to provide solutions to control air temperatures in the underground network and prevent temperatures rising to undesirable levels.

Keeping the Tube's customers cool involves developing new technologies and making the best use of more traditional approaches. Every effort is being made to use environmentally friendly methods wherever practical and several cooling options are being considered. Cooling benefits may be provided by the reduction of energy loss from trains and train traction systems, and from using mid-tunnel and station ventilation shafts. Wet cooling systems using water from the River Terrace Gravels and groundwater from the Chalk aquifer are also proposed. Chilled water can also be obtained from central cooling plants or from air-cooled chillers. Ventilation solutions (that offer both temperature and airflow improvements) are generally preferred but, where these systems are not viable (usually because of high costs of new shafts), there is recourse to water-based cooling systems. At sites where the groundwater quality and flow from the River Terrace Gravels are unsuitable, the use of borehole cooling solutions abstracting water from the Chalk are considered.

Groundwater cooling solutions were initially considered as one of the options at up to 20 LU stations, mostly in central London (Fig. 1). This would potentially be one of the largest schemes to use a major aquifer for cooling purposes in the world. As the relevant regulatory body (the Environment Agency of England and Wales) considers the water resources of the Chalk aquifer in London to be broadly in balance at present (Fry 2009), the proposals were for non-consumptive use of the water. The cool abstracted groundwater would pass through heat-exchanging pipe networks within the stations and tunnels, and sensible heat would be removed into the circulating water to chill the surrounding air within the LU network. The warmed water would be returned to the aquifer via injection boreholes.

At LU stations where borehole cooling systems are proposed, several other practical considerations determine the viability of any given scheme. The station must be located where it is suitable for the installation and associated infrastructure and where it is practical to route pipework between the abstraction and re-injection wells (with a desired separation of around 200 m). The selection of a suitable location is also determined by the likelihood of obtaining the desired borehole yield, where the desired cooling load can be met and where there is little likelihood of any thermal interference between the abstraction and re-injection borehole or with other nearby schemes.

The use of groundwater, abstracted via open-loop systems, is also increasingly being considered as a means to provide heating (and cooling) of large developments, particularly in London. At present, aquifers can be readily exploited for net extraction or disposal of heat, but as the number of developments using this technology increases so too does the number of heat sources and sinks in an aquifer. Where a dense concentration of such systems is used, injected water may ‘thermally pollute’ the abstraction borehole, and the groundwater abstraction wells may start to abstract water that has already been warmed (or cooled), potentially reducing their efficiency or affecting other users. LU therefore commissioned a programme of groundwater modelling to assess the potential hydraulic and thermal effect of the proposed schemes as well as the feasibility of the schemes for the expected cooling loads.

This paper provides an overview of the approach to modelling and summarizes some of the findings that may be of broader interest. In particular, this paper identifies critical factors that need to be taken into account when modelling heat transport in the groundwater system under London and key parameters that are likely to control the long-term efficiency of open-loop groundwater cooling schemes.

This paper is one of a thematic set entitled Hydrogeology and Heat Engineering. To avoid unnecessary repetition between papers in the set, the background to the science of thermogeology and the exploitation of ground source heat has been described by Banks (2009). The context of ground source heat in the UK regulatory environment and a comprehensive description of the regional hydrogeological conceptual model of the London Basin have been provided by Fry (2009).

## Approach to modelling

The objectives of the hydrogeological modelling were to: (1) make an initial assessment of the feasibility and performance of cooling schemes to support design decisions; (2) identify key parameters and help to design trials and experiments to reduce uncertainty; (3) provide a tool to support regulatory approval for sites for which ground source cooling was identified as the optimum solution.

The modelling presented here was carried out in advance of site investigation and most of the input parameters are based on literature data and limited site-specific data. The models developed at this stage would need to be further refined once the designs of any cooling schemes are finalized to meet the overall modelling objectives (a pilot scheme at Green Park is currently under development to provide more detailed site-specific characterization data). At this initial stage, the focus was on assessing the sensitivity of the model results to aspects of the conceptual model and parameterization.

A variety of modelling tools were used in a phased manner. Initially, simple analytical solutions were employed to gain an appreciation of the processes of heat transport within the Chalk and key uncertainties that control the feasibility of the groundwater cooling schemes. Preliminary design options were investigated and the sensitivity of results to parameters was assessed. These results were taken forward to distributed numerical models.

A total of four sub-regional groundwater models were constructed (see Fig. 1 for locations). FEFLOW (Diersch 2005) and MODFLOW/MT3D (Zheng & Wang 1999) were tested in parallel. The first models constructed (the ‘Stockwell’ and ‘Central London’ models) had rectangular outlines to allow equivalent FEFLOW and MODFLOW models to be built. Subsequent model domains were not constrained by rectilinear grids once a decision was taken to rely on FEFLOW, and the accuracy of the translation between the conceptual and numerical models was improved by using groundwater contours and flow lines as model boundaries.

## Conceptual model

A conceptual model of the relevant aspects of the hydrogeology of the London Basin has been presented by Fry (2009). The key elements of relevance to this paper are as follows (see Fig. 2).

(1) The strata beneath London form a shallow synform known as the London Basin. The beds generally dip towards the lowest point in the basin, which is broadly aligned along the River Thames.

(2) Formations of note include the Cretaceous Chalk, the Palaeogene Lower London Tertiary units and the Eocene London Clay. In central London, the Lower London Tertiary units comprise the Woolwich and Reading Beds with a variable thickness up to 30 m and a thinner layer of underlying Thanet Sands. There is a substantial thickness (up to 100 m) of low-permeability London Clay throughout the central areas of the London Basin, which forms a confining layer for the aquifer system. Above the London Clay, Quaternary Superficial Deposits mainly comprise River Terrace Gravels and other Alluvial Deposits.

(3) The Chalk forms the main aquifer of interest for the proposed open-loop ground source cooling systems, which is largely a result of the fractured nature in the uppermost tens of metres and the development of significant secondary porosity and permeability (Ellison *et al*. 2004). The Chalk below this fractured zone is of low productivity and low permeability and is considered to act as an aquitard. In central London the Chalk is typically between 170 and 200 m thick.

(4) The Chalk aquifer of central London was originally fully confined. However, high rates of abstraction starting in the 19th century resulted in a large fall in groundwater levels such that some parts of the aquifer have become unconfined. The decline in abstraction rates in central London in the second half of the 20th century resulted in some recovery of groundwater levels, but new management strategies have now achieved water table elevation stabilization (Fry 2009). Where the upper, more permeable, horizons of the Chalk aquifer are still unsaturated, this is likely to significantly decrease the transmissivity of the aquifer (in comparison with saturated conditions).

## Parameterization

The Chalk is a dual-porosity aquifer in which groundwater in the matrix blocks is essentially immobile because of the small size of the pore throats. Significant permeability is generally developed within the upper 50 m of the aquifer and is associated with fractures or joints (primary fissures) through which the majority of groundwater flow occurs (Allen *et al*. 1997; Ellison *et al*. 2004). Groundwater in the Chalk is normally considered to flow in corridors of high-permeability Chalk (fissures) separated by low-permeability blocks (matrix) (Downing *et al*. 1993). The permeability of the Chalk is highest in areas of greatest fracturing, normally where the overburden is thinnest and groundwater flux greatest. At greater depths within the Chalk, the frequency and aperture of fissures declines as a result of the increasing overburden and a general reduction in circulating groundwater.

To provide an estimate for the Chalk transmissivity in the model, data from 491 pumping tests in the Upper Chalk were obtained from Allen *et al*. (1997). Data obtained from pumping test scale or above were used, because smaller-scale tests represent a volume of rock smaller than the minimum representative elemental volume of the system under consideration here. Pumping test information was also collected from borehole logs and transmissivities were derived from the interpretation of yield tests using Logan's Rule (Logan 1964), which provides a crude estimate of the Chalk transmissivity from abstraction rate and drawdown. The data imply that there is significant variation in transmissivity over the study area (Fry 2009). In addition, there is an absence of yield data in some parts of the model area, suggesting that there is low potential yield and that transmissivity may be significantly lower here (resulting in fewer boreholes and fewer pumping tests being performed).

From work undertaken in this study and previous work (Allen *et al*. 1997; Ellison *et al*. 2004), it is considered that a value of transmissivity of around 300 m^{2} day^{−1} is a representative average for the confined Chalk in central London. Although pumping tests gave variable values of transmissivity across the model area, it is not clear how spatially representative these results are. It was therefore considered more appropriate to use a constant transmissivity value across the model area. The estimated likely range of the transmissivity is 50–500 m^{2} day^{−1} for boreholes penetrating 50 m into the Chalk and this was tested during the sensitivity analysis.

To estimate values of hydraulic conductivity of the Chalk it is necessary to consider the likely distribution of the transmissivity throughout the formation. Two plausible end-members for the distribution of transmissivity in the Chalk aquifer in central London are considered, as follows.

(1) The majority of the transmissivity is concentrated in the top few metres of the Chalk: if 80% of the transmissivity is in the top 5 m, this would give a hydraulic conductivity of 48 m day^{−1} for the top 5 m and 1.33 m day^{−1} for the lower 45 m of the Chalk, for a borehole penetrating 50 m into the Chalk.

(2) The majority of transmissivity in the Chalk is distributed evenly throughout the upper 50 m: this gives a typical hydraulic conductivity of 6 m day^{−1} if it is assumed that the transmissivity (300 m^{2} day^{−1}) is uniformly distributed over this interval.

Gropius (2010) used a similar approach in testing different conceptual models of a doublet system installed in the Chalk of London. He found that concentrating most of the transmissivity in a 10 m thick layer gave model results that matched the observed behaviour of the doublet system better than focusing it in a 1 m thick layer. The relevant hydraulic properties of the formations simulated by the numerical groundwater models are summarized in Table 1. Below the first 50 m of Chalk there is little information on the hydraulic conductivity as few boreholes penetrate to this depth. Hydrogeologists understood that the top 50 m or so of the Chalk was the zone of the majority of groundwater flow and drilling deeper would see decreasing benefits in borehole yield.

No site-specific data on porosity, thermal conductivity or heat capacity were available for the various geological formations and literature values for equivalent lithologies or regional data were used (see Table 2). As a first estimate, the effective thermal conductivities and heat capacities used (those of saturated rocks) were based on dry rock type or mineral specific heat capacity and conductivities, calculated as bulk parameters for an aquifer matrix assuming a total water-filled porosity of 30%. These calculated values are typical, or within the range, of the literature data for similar lithologies (see, e.g. Busby *et al*. 2009; Headon *et al*. 2009). A wide range of effective heat capacity and thermal conductivity values (as well as hydraulic properties) were assessed during the modelling work (see Table 3).

The Superficial Deposits, London Clay, and the Woolwich and Reading Beds were provisionally given the same thermal properties as each other, although the Woolwich and Reading Beds were considered separately during sensitivity analysis. The River Terrace Gravels are likely to have different thermal parameters from the underlying London Clay, but the separation of abstraction and injection points in the Upper Chalk and the thickness of the London Clay meant that this is not expected to have significant effects on the simulated heat distribution.

Primary (or matrix) porosity of the Chalk is fairly large (around 25–40%), although very small pore throats result in the primary matrix permeability being low, *c*. 10^{−3} m day^{−1} (Price *et al*. 1993). Typically, the matrix porosity (including non-interconnected and dead-end pores) does not contribute significantly to groundwater flow. The porosity through which flow can occur is the effective porosity and, in the Chalk, a good first approximation is generally to use a value similar to the specific yield, or one calibrated from site-specific data. In the unconfined Chalk, specific yield commonly varies between 0.01 and 0.02. For the confined Chalk a value of 0.01 or less is generally used. In the upper, high-permeability Chalk and the underlying lower-permeability Chalk, effective porosity has been calculated based on an average Chalk effective porosity of 0.01.

The value of longitudinal dispersivity was chosen as 1/10 of the distance between the injection and the abstraction borehole (this had been chosen as the relevant transport distance; see Gelhar 1993). Transverse and vertical dispersivities of 1/100 and 1/1000 of the transport distance were used respectively (FEFLOW does not distinguish between transverse and vertical dispersivity and only the transverse value is used).

## Analytical model

### Scoping calculations

A key issue is the degree of interaction of heat with the Chalk matrix. Within the flow conducting fissure system to be addressed in the modelling, two extremes may be considered: one in which heat acts as a conservative tracer and no conductive transfer of heat occurs, and one in which heat fully equilibrates with the Chalk matrix. The real system is likely to be somewhere between these extremes, with equilibrium between fissure flow temperature and intact Chalk matrix temperatures taking some time to develop. If the equilibrium time scale is relatively short, then an approach that includes the effects of heat conduction to the Chalk matrix is likely to be the most accurate. If the flow conducting fissures are widely spaced, or the contact between flow in the fissures and the rock matrix is reduced (by channelling, for example), then not all of the rock matrix will be heated and a conservative heat tracer representation may be more appropriate.

A number of analytical solutions exist for heat transport based on solutions to the diffusion equation and are useful for scoping and quick sensitivity calculations (see, e.g. Tang *et al*. 1981; Sudicky & Frind 1982; Lever *et al*. 1983; Goebes & Younger 2004). Initial scoping calculations were undertaking using the analytical solution of Tang *et al*. (1981). This analytical model applies to 1D heat transport along a discrete fissure from which an infinite porous matrix extends. Advection only is considered in the fissure, and diffusion occurs within immobile water within the matrix.

Table 4 summarizes the parameters used to populate the model and the results are shown in Figure 3. Figure 3a shows the temperature breakthrough at a point in the fissure 100 m downstream from the heat input. The value of the wetted fraction is varied (between 0.01 and unity) to investigate the effects of the fissure-to-rock matrix contact area. This represents the accessible fraction of the matrix porosity and relates to the proportion of (infinite) fissure plane that actually conducts flow. Figure 3b shows a temperature profile in the rock matrix at this location with a wetted fraction of fissure face of 0.01. An effective bulk diffusion coefficient (diffusivity) is used to represent matrix diffusion, determined from the quotient of the effective thermal conductivity and heat capacity. The diffusivity is the rate at which heat can be conducted through the saturated porous medium: in this case 5.82 × 10^{−7} m^{2} s^{−1} (thermal conductivity and heat capacity data presented for the Chalk by Headon *et al*. (2009) give a slightly higher diffusivity value at 7.27 × 10^{−7} m^{2} s^{−1}).

The results show that, after a period of 10 years, the breakthrough temperature (*T*) is predicted to be about half the temperature of the injected water (*T*_{0}) at 100 m with a wetted fraction of fissure face of 0.01. As the wetted fraction increases, temperature breakthrough rapidly declines (at distances in excess of 500 m, heat breakthrough does not occur, even for a wetted fraction of fissure face of 0.01). The results suggest that, for this simple system, and given the high thermal diffusivity of the Chalk matrix (relative to standard tracers), thermal diffusion allows heat to be conducted into the Chalk matrix away from the fissure flow contact areas (channels) even with a low wetted fraction of fissure face.

Given the expected fissure frequency, time scales and borehole separation, it seems likely that diffusion into the rock matrix continues over large distances in excess of the real fissure spacing (the fissure separation is used in the analytical solution to convert the hydraulic conductivity into an average fissure transmissivity, not to provide a boundary condition limiting diffusion into the matrix away from the fissure face). This is clearly unrealistic: the analytical model does not demonstrate that there will be no breakthrough of heat at an abstraction well, rather that all of the rock matrix between flow conducting fissures is likely to be heated, supporting models that incorporate full interaction with the matrix.

### Tracer testing

Results from the scoping calculations suggested that heat diffusion within the Chalk matrix was important. The tests at the Royal Festival Hall (RFH) in June–July 2007 (Clarkson *et al*. 2009) provided a first opportunity to calibrate the simulation of heat transport and to provide some quantification of the effects of the Chalk matrix and confirm or deny the key findings of the scoping calculations. Fluorescein, salt (NaCl) and heat tracer tests were conducted between a test abstraction and injection borehole with a separation of 144 m at an average flow rate of 8.3 l s^{−1} over a period of 29 days. Injected water was heated to between 8 and 10 °C above the abstracted water temperature (over the duration of the test assuming an average 8 °C temperature rise, this equates to a total heat input of 194.1 MWh).

Data from the RFH tests were obtained and modelling of the fluorescein test results was conducted to obtain information on the hydrogeological and solute transport properties of the Chalk. This was undertaken by using an analytical model for heat and solute breakthrough developed by Barker (2010). This analytical model applies to a regularly fractured medium in a borehole dipole flow system and represents a more realistic simulation of the test results than the Tang *et al*. (1981) solution. The analytical model assumes a steady-state pumping rate and that water pumped from an abstraction borehole is injected at an equal rate into an injection borehole at a given separation. The injected water has a fixed increase in temperature (or concentration) above background and the formation is considered to have constant hydrogeological properties. Groundwater can flow only along the fractures (i.e. no matrix flow); heat (or solute) can diffuse into the formation around the fracture.

#### Fluorescein model

The results of the laboratory testing of samples from the tracer test were used to calibrate the analytical solution, as these were more reliable than those from the in-line fluorimeter (Clarkson *et al*. 2009). The input parameters for this simulation are given in Table 5.

The breakthrough time is given by the following equation (Barker 2010):

Because the distance between the boreholes (*D*) and the pumping rate (*Q*) are well constrained, this means that the product of the number of fractures and their aperture (*N*_{f}*a*) is constrained by the observed breakthrough time; that is, the first time at which elevated fluorescein concentrations are observed, 0.26 days.

It was found that the results were not sensitive to the formation thickness (*B*) for the expected range of thickness for the Chalk (for *B* greater than 0.2 m, the breakthrough curve is unchanged). This is because the fluorescein does not diffuse far enough into the matrix for it to penetrate the full half-thickness of a matrix block and start being affected by solute diffusing from the adjacent fractures. The results were also insensitive to the tracer input period over the probable range (a range of 15 s to 20 min made no discernible difference to the results).

To calibrate the curve to the observed fluorescein concentrations, the parameters matrix porosity and number of fractures were varied (fracture aperture was calculated from the observed breakthrough time and the number of fractures using the above formula). It was assumed that the diffusivity of fluorescein and the pumping rate were well constrained. The results are shown in Figure 4 for the parameter set shown in Table 5. There is a degree of non-uniqueness of the calibration as a higher number of fractures and a lower matrix porosity result in similar breakthrough curves. For the possible range of Chalk matrix porosity, say 0.2–0.5, the number of fractures would be between five and 10.

The recovery rate in the field test was 40% compared with a simulated recovery of just over 50%. It is not clear whether the fluorescein has been sorbed onto the aquifer or degraded, if there was analytical error, or whether the unrecovered tracer was simply not captured by the pumping borehole. If sorption or degradation has occurred, then the values assigned to the matrix porosity, number of fractures and fracture aperture required to calibrate the breakthrough curve would change. For comparison, in a 3D distributed model, to obtain a good calibration with the fluorescein tracer experiment, 100% of the Chalk transmissivity was attributed to the top 5 m of the Chalk and a low effective porosity of 0.0025 was applied (Clarkson *et al*. 2009).

#### Heat model

The heat tracer test was also modelled using the analytical solution, and the input parameters for this simulation are given in Table 5 and the results shown in Figure 5. The injection temperature should be measured in the bypass loop, which monitors the temperature of the hot injection water prior to any dilution effects before it enters the borehole. However, because of data gaps, it was necessary to use a combination of data including monitoring data from the injection borehole. An average injection temperature of 23.6 °C was estimated for the test, which is 9.6 °C above the ambient groundwater temperature at this locality. The peak temperature simulated at the abstraction borehole using the analytical solution was 0.37 °C above background, which occurs after 500 days.

Although the results of the analytical solutions were consistent with the test results (i.e. a rapid breakthrough of the fluorescein tracer was observed, yet there was no heat breakthrough on the time scale of the test), it has not been possible to quantify the significance of diffusive heat transport within the Chalk matrix from the data available. The analytical solution assumes a steady-state pumping rate (i.e. for an infinite time). In this case, the fluorescein breakthrough is fast and therefore this assumption is valid, as the pumping period is long compared with the breakthrough time. For a heat breakthrough to be seen, a longer test would be required. For example, at these flow rates and borehole separations (for the described parameterization) the analytical model predicts a peak temperature rise of 3.5 °C after 700 days for a 1 year test, and a peak temperature rise of 4.5 °C after 1000 days for a 2 year test (see Fig. 5). These time scales are longer than would be practically feasible for a pre-commissioning test, which suggests that constraining the relevant heat parameters of the Chalk will require monitoring of operational systems.

#### Comparison with distributed model results

The analytical solution was used to attempt to simulate the breakthrough of heat from LU cooling schemes, which were analysed over a 60 year time scale, to compare the results with the fully distributed models. Model results for a cooling scheme option at Camden Town are presented here for this comparison (Camden Town station is isolated from the hydraulic and thermal effects of other cooling schemes; see below).

In Camden Town, the top of the Chalk is unsaturated and advective flow cannot occur here, or within the overlying Thanet Sands (although it is still possible for heat to diffuse upwards through the unsaturated zone into the overlying formations). As a first approximation (to allow a comparison with the distributed model) the geological sequence in the analytical model was simplified to an equivalent thickness of Chalk (59 m) and London Clay (52 m); weighted average properties for the Chalk and Clay were used.

The number of fractures and the aperture size determine the half block size (this is the distance from the centre of a block of aquifer to the nearest fracture). Heat is transported into the London Clay only by diffusion from the top of the Chalk, so a single fracture would represent this heat transport (the equivalent half block size for the Clay would therefore be around 60 m). There are likely to be more fractures in the Upper Chalk, assumed to be evenly distributed throughout the aquifer thickness in the analytical solution (as such a half block size for the Chalk would typically be lower than for the London Clay). Breakthrough curves for varying half block sizes in comparison with the model results are shown in Figure 6.

The results show that the initial breakthrough of the distributed model was coincident with a half block size of around 18.5 m. For a single fracture, breakthrough occured at earlier times and increased with the number of fractures represented; as the formation thickness reduced, the temperatures increased. Although the results are not presented here, the solution was also insensitive to either the heat capacity or diffusivity within the ranges expected for Chalk and Clay. It was not possible to match the analytical solution to the model results by varying the input parameters.

The results suggest that heat transport on these longer length or time scales was not dominated by advective fracture flow in the top part of the Chalk and that the influence of the underlying and overlying strata was important. The geological sequence in central London cannot be simplified into a single regularly fractured medium and the analytical model was not suitable for modelling the tunnel cooling schemes. It was necessary to represent heat diffusion into the formations above and below more accurately, and a distributed model that can represent the detail of the geological setting was required.

## Numerical model set-up

### Model structure

Detailed geological information was taken from over 500 geological logs of boreholes within the model area. Five geological units have been incorporated in the models: the top Superficial Deposits (Layer 1), underlain by the London Clay (Layers 2–4), underlain by the Woolwich and Reading Beds (Layer 5), underlain by the Thanet Sands (Layers 6 and 7, each half of the total thickness), which in turn are underlain by the top 50 m of the Upper and Middle Chalk (Layers 8–18, 10 layers of 5 m thickness plus one layer below that of 50 m thickness), the lowest geological stratum considered in the model.

Given the hydraulic isolation of the River Terrace Gravels from the Chalk, the River Terrace Gravels were not modelled in detail (i.e. recharge and surface flows were not represented). The River Terrace Gravels and superficial deposits were included as the top layer in the model. Where River Terrace Gravels or high-permeability superficial deposits were inferred to be greater than 5 m thick, a high hydraulic conductivity was specified, and elsewhere the superficial deposits were specified as a low hydraulic conductivity clayey layer.

Very fine grids were required to demonstrate grid convergence of the results because there were very steep gradients close to the boreholes when fully resolved, leading to very small gradients in the region between the boreholes: coarse grids failed to fully resolve the gradients adjacent to the boreholes and similarly overestimated gradients at the central region between boreholes.

### Boundary conditions: flow

In each sub-regional model, either constant head or no-flow boundaries were defined along lateral boundaries in the Chalk, Thanet Sands, and Woolwich and Reading Beds, based on information taken from the Environment Agency groundwater level map (Environment Agency 2009). In the study zone hydraulic gradients were assumed to be between 0.01 and 0.001 with groundwater flow directions towards central London (locally there were variations in hydraulic gradient from model to model). Lateral boundaries for the London Clay and superficial deposits were set as no-flow.

The River Thames was included as a fixed head boundary at sea level in the top model layer and no information on the river bed elevation or properties was included. A significant thickness of low-permeability London Clay is present beneath the River Thames, which limits the effect of the River Thames on groundwater levels in the Chalk.

Up to six pairs of cooling scheme abstraction–re-injection boreholes are represented in each model. A pumping rate of 2160 m^{3} day^{−1} (25 l s^{−1}) is assigned to the abstraction and injection wells based upon provisional CTP scheme designs. Both abstraction and injection boreholes are represented as open over the upper 50 m (10 layers) of the Upper Chalk. FEFLOW simulates the boreholes as a highly conductive linear element connecting all the neighbouring boundary condition nodes for the well. During computation, the hydraulic head was uniform throughout the linear element (because the well conductivity is very large) and the abstraction rate from each layer was calculated according to the surrounding material properties.

There are about 20 licensed groundwater abstractions within the central London model area, although actual abstraction rates and borehole construction details were not available for this project. In view of this, only the largest abstractions (with annual licensed rates above 1000 m^{3} day^{−1}) were included in the models, pumping at their maximum permitted annual rate from the Chalk aquifer. It was clear that the abstractions close to the arbitrary domain boundaries would influence the groundwater contours at the boundaries. The influence of abstractions on the regional flow pattern was accounted for, to some extent, by using fixed head boundaries based on measured groundwater head contours.

### Boundary conditions: heat

Monitoring data provided by the Environment Agency indicated that annual groundwater temperatures at depth in the Upper Chalk in central London were typically between 12 and 14 °C. Recent work undertaken by Headon *et al*. (2009) has demonstrated that groundwater temperatures in the London area exhibit a reasonable amount of lateral variation. A median temperature value of 13.1 °C was derived at depths of 100 m, with an overall variation in temperature at the top of the Upper Chalk observed from east to west from <11 °C to >15 °C respectively.

In the UK the upper 15 m of the ground tends to exhibit a seasonal temperature variation in response to air temperatures (Busby *et al*. 2009). At greater depths it is generally assumed that these surface effects diminish and that temperatures increase with depth at a rate determined by the geothermal gradient. In contrast to temperature increasing along an ideal geothermal gradient, many temperature logs also exhibit negative temperature gradients, where temperatures initially decrease with distance from the surface but the gradients subsequently become reversed at depth (Banks 2009; Headon *et al*. 2009).

It is likely that observed temperatures are affected by the dense development in the city and local anthropogenic influences. Buried services, cable tunnels, sewers and downward conductive heat leakage from buildings, as well as groundwater abstractions, artificial recharge and historical groundwater levels, will cause local temperature variations. Natural groundwater gradients, the effect of variable recharge and groundwater surface water interactions are also likely to play a role. For these reasons there is considerable uncertainty associated with the assignment of thermal boundary conditions and initial temperature distributions.

Typically, over long time scales, it is reasonable to assume that air temperatures are translated to ground temperatures and that a constant temperature boundary condition on the upper surface of the model equal to the mean annual air temperature is a valid assumption. In London, the air temperature varies from about 6 to 19 °C with an annual average of around 11.7 °C (Headon *et al*. 2009). However, the role of built structures in the London Clay and other anthropogenic influences is unclear and these may give rise to higher (or lower) surface temperatures than would otherwise be assumed. It is likely that surface losses from the cooling schemes in largely urban areas would be lower than from comparable systems in areas where fewer buildings are present, where buildings are generally maintained at temperatures in excess of air temperatures, and where downward conduction of heat has been demonstrated to occur (Banks 2009). In central London it may be appropriate to assume that there are no surface heat losses, or that surface temperatures are generally in excess of air temperatures.

Boreholes in which the groundwater temperatures have been monitored (described above) tend to be completed at a depth of about 50 m into the Chalk. Few data are available beneath this depth in central London (at depths in excess of 100 m below ground level). Preliminary modelling indicated that the effects of downward heat flux beneath the level of the boreholes were important, and that underlying formations should be included. As such, a 50 m layer of the Middle Chalk is included below the base of the abstraction and injection wells, and this assumption was tested in sensitivity analysis. Estimates of temperatures at the base of the Middle Chalk can be made from an assumed geothermal gradient. Barker *et al*. (2000) estimated the temperature at 7 km depth in London to be around 130 °C; this equates to a geothermal gradient of around 1.66 °C per 100 m depth (0.0166 °C m^{−1}). This would suggest an increase above the ambient temperature at the model base of *c*. 0.83 °C. Other estimates of the geothermal gradient of 0.026 °C m^{−1} and 0.035 °C m^{−1} (Busby *et al*. 2009) provide temperature rises of 1.3 °C and 1.75 °C respectively. The uncertainty in the geothermal gradient in the Chalk amounts to *c*. 2 °C per 100 m without considering the effects of anthropogenic influences that may affect this gradient locally. For example, downward heat conduction would mean that the geothermal gradient is comparatively lower than in areas not subject to such influences (negative gradients near the surface indicate that temperatures at depth have not been increasing at the same rate as those at the surface).

It was not immediately apparent whether uncertainty associated with initial temperature distribution or the assignment of the upper and lower thermal boundary conditions would have a larger impact on the simulated abstraction temperatures of the cooling schemes. To address some of this uncertainty the effects of alternative boundary conditions were investigated during the development of the conceptual model.

The initial temperature distribution in the model was assumed to be a uniform 14 °C. This has largely been a pragmatic decision based on the available data: given the overall uncertainty, this appears to be a generally representative value of groundwater temperatures in central London. In assuming a constant and uniform temperature distribution the imposition of any fixed temperature boundary conditions is therefore constrained. Assigning an increase of temperature at depth or decrease of temperature at the surface would require that the initial temperatures are set at a constant gradient and that a steady-state (or transient) model would need to be implemented to calculate these initial conditions. Given the uncertainty, this is an additional model complexity that is hard to justify and, at the sub-regional scale, there are too few data against which the model can be calibrated. Adopting a constant initial temperature distribution in the model also has the advantage of being a stable baseline temperature upon which the impact of the cooling schemes can be imposed and from which the relative impacts of the various schemes can be determined.

The River Thames was specified with a constant temperature in the river of 14 °C although this assumption was tested in sensitivity analysis. The presence of high-permeability deposits above the London Clay, which may influence the heat loss at the surface (as a result of flushing with surface water or rainfall), was accounted for by using a fixed temperature boundary condition at 14 °C where high-permeability deposits were present. No-diffusive-flux boundaries were set elsewhere on the top surface of the model (this is a boundary condition that allowed heat to be advected out of the solution domain). As there was no advection out of the upper surface, the assignment of a no-diffusive-flux boundary essentially imposed a no-heat-flux boundary where heat may accumulate. This may be a conservative assumption for the impact of the cooling schemes. To address this uncertainty, a fixed temperature boundary condition was also applied to the model surface in sensitivity analysis. Results demonstrated that the effects of this boundary condition may be significant only at longer time scales (>60 years), when heat is conducted this far from the Chalk injection borehole through the significant thickness of overlying London Clay.

The lower boundary was specified as a fixed temperature boundary at 14 °C to prevent heat from accumulating at the base of the model, and this assumption was further tested in sensitivity analysis. Heat can therefore pass out of the lowest geological unit represented in the model, the Middle Chalk, into the underlying formations. However, this is located 50 m below the base of the boreholes, at a sufficient distance that, within the time scale of the simulations, there was not a significant amount of heat leaving through the model base. No-diffusive-flux boundaries were set at the lateral edges of the model domain. Heat can be advected out of the solution domain should the plume extend that far, but will generally not have an influence on the transport of heat near the centre of the model domain in the region of interest.

### Heat injection

Heat generation within the tube tunnels occurs from a variety of sources. In the carriages the largest proportion of heat is generated by the passengers; loadings are variable, but it is estimated that heat generation may be up to 15 kW for a peak loading of up to 100 passengers per car (Gilbey & Thompson 2009). The major source of heating within the stations and tunnels is from the trains themselves from frictional losses, losses through drive inefficiencies and heat generated during braking. For a typical section of track at an interchange station, losses have been estimated at 393 kW. This varies considerably according to the characteristics of each station, train type, speed, loading and service frequency; a proportion of these heat losses is also removed by ventilation and via dissipation of the heat into surrounding strata.

Detailed transient loading requirements and seasonal temperature and cooling demands for each station are determined using detailed tunnel ventilation models. These models are calibrated against historical data and predict an excess of up to 6 °C above the desired platform and tunnel temperatures during the summer months allowing for the future growth of LU services (Gilbey & Thompson 2009). The amount of heat removal required ranges between 200–300 kW peak load for a single-line station and 600–1000 kW for an interchange station. Taking into account diversity and load fluctuations, this might be of the order of 750–1000 MWh per year for a single-line situation. Operating at maximum theoretical efficiency for 1 year, assuming a reduction of platform temperatures of 6 °C and a concomitant increase in the temperature of injected water, the open-loop schemes would provide a cooling capacity of 5500 MWh.

It is expected that the demand for cooling may change over the longer term as a result of changes in the cooling requirements. Preliminary models (using an annual sinusoidal signal around the average) demonstrated that seasonal fluctuations in the borehole re-injection temperatures were unlikely to result in substantial variations in abstraction temperatures at this separation and scale because of the thermal capacity of the Chalk and overlying formations. An annual (seasonal) variation in thermal load was therefore not represented and a single annual average re-injection temperature was applied. It has also been assumed that the injection and abstraction boreholes operate continuously throughout the model simulations.

An approximate annual thermal demand upon the aquifer is based on an estimated fraction of the calculated peak thermal demand established for each borehole cooling system derived from the tunnel ventilation models described above; this is represented by the modelled water re-injection temperature, given by where *T*_{inj} and *T*_{abs} are the injection and abstraction temperature (°C), *T*_{max} is the maximum injection temperature (°C), and *f* is a pumping regime factor that defines the increase in temperature of the re-injected water. At the time of this study, the maximum permissible injection temperature imposed by the Environment Agency was a 10 °C rise above the ambient water temperature (i.e. *c*. 24 °C; see Fry 2009).

A pumping regime factor of 5/6 is the CTP best estimate of the annual average cooling demand (5/6 of the summertime peak cooling demand). Reductions in the pumping regime factor (e.g. 4/6, 3/6 in Table 6) represent a reduced load or smaller borehole cooling systems. Although seasonal variations are not explicitly represented, the input boundary condition (using the above equation) allows for the injected water temperature to increase with time in response to an increasing abstraction water temperature; this would imply a gradual loss of cooling capacity. The effect of increasing injection temperature with time can be modelled in FEFLOW using a solution-dependent boundary condition. This was implemented in station-specific models to investigate the sensitivity of results to various pumping–heat scenarios.

The predicted average injection temperature rise for the 5/6 scenario is fairly constant and can be modelled using a steady injection of water at a constant temperature with little loss of accuracy. This is a convenient representation of the conditions at the injection well and was used to assess the combined impact of the cooling schemes at the sub-regional scale (as well as to allow comparison between MODFLOW/MT3D and FEFLOW). For an initial abstraction temperature of 14 °C, the increase in injection temperature will be 8.33 °C. Over the lifetime of the schemes (60 years), scoping calculations suggest that the increase in abstraction temperature (at a separation of around 200 m) will be around 3 °C (i.e. up to 17 °C). Assuming that the abstraction temperature has increased by this amount, then the final injection temperature will be 8.83 °C above the ambient groundwater temperature. This gives an average temperature increase of 8.6 °C over the 60 year simulation (i.e. a constant injection temperature for these simulations of 22.6 °C).

## Numerical model results and sensitivity analysis

### Best estimate model results

This section presents the results of the simulations undertaken with the FEFLOW model of the central part of the study area (Fig. 1). The results are for the Best Estimate Model (BEM), which simulates the Chalk aquifer with an upper layer of increased hydraulic conductivity and effective porosity. The sensitivity of the model to various assumptions and parameter values is discussed below.

Figure 7 shows a cross-sectional view of the heat plumes around the injection boreholes after 60 years of continuous pumping. The plumes are fairly symmetrical about the injection wells in the vertical plane, which suggests that heat is transported mainly as a result of heat conduction within the Chalk layer. In the uppermost Chalk layer, the increased hydraulic conductivity enhances the rate of heat flow from the borehole to the aquifer and the influence of heat transport by advection becomes more significant in the vicinity of abstraction boreholes, where groundwater flow velocities are high. Nevertheless, the dominant heat transport mechanism remains heat diffusion through the rock matrix.

Figure 8 shows a plan view of the simulated extension of the heat plumes in the uppermost Chalk layer for the BEM after 60 years. This illustrates the potential impact on the surrounding aquifer and hence, potentially, on other users.

The simulated temperature at the abstraction boreholes for each dipole pair is presented in Figure 9. The figure shows breakthroughs in the uppermost Chalk layer (Layer 8) at the abstraction boreholes. The average temperature breakthrough is generally slightly lower as a result of the increased advective transport of heat along this higher permeability layer.

The distance to the nearest injection borehole for the most part determines the time to temperature breakthrough, and to a large extent the eventual steady-state temperature rise. The feasibility of schemes at varying separations was examined in more detail as part of the sensitivity analysis, which focused on four main aspects of the modelling, as discussed below: conceptual model, parameterization, operational aspects and model code.

### Sensitivity analysis: conceptual model

In the absence of site-specific data, two conceptual models of the distribution of hydraulic conductivity in the Upper Chalk were considered: a ‘heterogeneous’ model (BEM), in which 80% of the transmissivity is in the top 5 m of the Chalk; and a ‘homogeneous’ model, where the transmissivity is distributed throughout the Chalk formation thickness. The effect of a fixed temperature on the upper and lower model surface of the BEM was examined. The effect of the position of the lower boundary condition was also investigated by subtracting a 50 m layer of Chalk below the base of the boreholes.

The results of these sensitivity runs are shown in Figure 10 as ‘breakthrough’ curves of the average temperature at one of the abstraction boreholes, located *c*. 300 m from the injection well. In the BEM (Model 1), the final temperature simulated after 60 years is 18.5 °C.

A sensitivity run was carried out without a fixed temperature of 14 °C at the lower boundary (Model 2). At early times, the difference in the basal boundary makes little difference, as the heat is going into storage within the basal layer of the model. In the BEM, at later times, heat is transported through the base of the model, which leads to a decrease of about 1 °C in the final temperatures in comparison with Model 2. The removal of the lower Chalk layer (Model 3) increases the final simulated temperature by around 2 °C in comparison with the BEM. There is no heat flux across the boundary at the base of the boreholes into underlying formations and this temperature difference can be attributed to the accumulation of heat at the model base.

The addition of a fixed temperature boundary at the top of the model (Model 4) makes little difference to the results because it is separated from the boreholes by around 50 m of low-permeability and lower thermal conductivity deposits above the Chalk. Only a small amount of heat is lost from the surface at later times.

The ‘no unsaturated zone’ run (Model 5) removes a low-permeability zone in the top layer of the Chalk representing an area of unsaturated Chalk (Fig. 5 and Table 1). In the ‘homogeneous’ Chalk model (Model 6), the breakthroughs appear later, and the final temperatures are much lower as the heat will have accessed the heat capacity of the full Chalk formation.

### Sensitivity analysis: model parameters

Modelling of the central London area included detailed sensitivity and uncertainty analysis. This was used to determine the models' sensitivity to all the input parameters. The highest sensitivity was to the volumetric heat capacity, thermal conductivity and hydraulic conductivity. These were then combined into a set of ‘best case’ and ‘worst case’ parameters to provide the range of uncertainty for the subsequent modelling work. A range of effective porosity values was used to determine the range of *C*_{v} and κ_{v}, but was not varied in the model itself. The parameter values used in the ‘best’ and ‘worst’ case scenarios are given in Table 3.

The resultant final temperatures are shown in Table 5 (the operational scenarios shown in this table have been described previously). The model results suggest that, in general, the uncertainty in results owing to parameter uncertainty is of the order of ±0.5 °C. This is less than the range of uncertainty associated with the conceptual model (total range about 2 °C).

### Sensitivity analysis: operational aspects

At one site the effect of increasing injection temperature with time for different pumping regimes and variable re-injection distances was examined to assess different scheme feasibility options. Four re-injection temperature scenarios, which use feedback of the abstraction temperature to control the re-injection temperature, were tested (i.e. as the abstraction temperature increases, so does the injection temperature). Each pumping regime was assessed at each of four proposed injection locations, at distances of 150, 255, 350 and 500 m from the abstraction borehole.

The final temperatures for the 16 scenarios are summarized in Table 7. These results are typical of the reduction of the abstraction temperature seen with increasing borehole separation. For boreholes with 500 m separation, the increase in temperature above background at the abstraction borehole is 2.1 °C after 60 years, compared with 4.6 °C for a borehole separation of 255 m.

### Comparison between FEFLOW and MODFLOW/MT3D results

For the Stockwell and Central London models, a parallel model was constructed in MODFLOW/MT3D. These were in most aspects the same as the FEFLOW model and the differences are listed below. The relationships used to simulate heat transport in MODFLOW/MT3D are given in Table 8.

MODFLOW/MT3D does not simulate heat transport above the water table, which was artificially raised by 100 m, to ensure saturated conditions across the model. This does not affect the groundwater flow field, as all boundary heads and initial heads have been increased consistently and the transmissivity is unchanged as the layers are specified as confined. MODFLOW/MT3D also assumes that the heat difference between the ambient groundwater temperature and the expected temperature of the injected water does not involve a significant change in density or viscosity of the groundwater, and that the groundwater flow remains unaffected by buoyancy-driven flow. For most cases, it is assumed that a temperature difference of up to 10 °C does not involve a significant change in density or viscosity of the groundwater (between 12.5 and 22.5 °C the density of water changes only by about 0.01%; the viscosity decreases by about 17% with a reciprocal rise in hydraulic conductivity).

Figure 11 shows a comparison of the MODFLOW/MT3D and FEFLOW results for constant density and viscosity simulations and the results of a coupled FEFLOW model representing density- and viscosity-dependent flow. This confirms the results of previous modelling, which had shown compatible results between the two codes and that density and viscosity changes over the simulated temperature range have only a very small effect on the transport of heat.

## Summary and conclusions

This paper presents a structured approach to modelling the potential effects of large-scale open-loop ground source cooling schemes in central London. A tiered approach of increasing levels of complexity has been used to guide and refine the approach to modelling.

Initially a simple 1D analytical solution (Tang *et al*. 1981) was used to investigate the processes of heat transport in the Chalk. The results suggested that heat would be conducted into the Chalk matrix even for low fissure-to-matrix contact areas as a result of high Chalk thermal diffusivity, supporting the use of heat transport models that incorporate full interaction with the matrix. A more complex analytical solution for heat and solute transport in a regularly fractured medium in a borehole dipole flow system (Barker 2010) supported these preliminary findings. The analytical solution was also used to derive parameters from a field-scale test (Clarkson *et al*. 2009) and compared with results obtained from distributed numerical models. This approach is useful for characterizing aquifer properties relevant to tracer tests at this scale, but because of the longer time required to observe heat responses, which brings into play significant interaction with vertically adjacent strata, it is of less relevance to scoping heat flow in open-loop systems for the confined Chalk in central London.

Distributed numerical groundwater models were constructed to simulate the interaction within groups of ground source cooling systems using the well-established conceptual model of the hydrogeology of the London Basin (Fry 2009). Because of the lack of site-specific data, extensive sensitivity analysis was carried out to assess the level of uncertainty in the model predictions. Sensitivity was measured relative to a best estimate model, considered to represent the best estimates of the parameters and conceptual model. The greatest uncertainty was associated with detailed aspects of the conceptual model of heat flow in this system.

The best estimate scenario took the rather conservative assumption that 80% of the transmissivity for a 50 m deep borehole is located in the top 5 m of Chalk. This represents a somewhat extreme concentration of the transmissivity at the top of the Chalk on the regional scale. This assumption would imply the presence of a continuous layer of high hydraulic conductivity over the entire model domain with a hydraulic conductivity of about 50 m day^{−1}. If flow is concentrated within a small volume of the aquifer, the potential for fast paths between the borehole pairs is increased, thus increasing the risk of thermal breakthrough and reduced efficiency of the cooling scheme. However, on the scale of the proposed cooling schemes, the dominant heat transport mechanism was diffusion through the rock matrix.

The role of over- and underlying strata is important, and it is necessary to characterize these layers and incorporate them in the models with sufficiently detailed grids to resolve the heat transport through them. There is significant uncertainty associated with the assignment of appropriate thermal boundary conditions. At the base of the model there is little detail regarding the geological system below the Upper Chalk, or temperatures at depth. Similarly, surface conditions will be complicated by the built environment at the surface, and utilities in the near surface. There were no data available to the project on the temperature distribution at shallow depths or within the London Clay.

Operational factors such as the amount of heat injected and the spacing between abstraction and injection boreholes are clearly important in terms of predictions of future efficiency of these schemes and the potential for impacts on other users. The modelling presented here provides some quantification of the significance of these factors.

This work has demonstrated that equivalent MODFLOW/MT3D and FEFLOW models give good agreement and are both feasible for modelling heat transport, as density-driven flow is not significant for temperature differences of the scale used by these cooling schemes. The unsaturated zone is important for heat transport. Although this can be a limitation if MODFLOW/MT3D is used, this can be avoided by using a water table approximation to ensure saturated conditions.

The work presented here summarizes the results obtained from an initial modelling phase as part of extensive investigations being undertaken by the CTP. As part of these continuing investigations, a pilot scheme is currently in development at Green Park to provide more detailed site-specific characterization data. The CTP constantly assesses the suitability of cooling schemes as new information becomes available. Future modelling work will become increasingly refined as additional data become available to meet the overall objectives of the project scope.

## Acknowledgements

The work reported here was commissioned and funded by London Underground Ltd., whose permission to publish this paper is gratefully acknowledged. The project was managed by A. Groves of Parsons Brinkerhoff and valuable advice was provided by M. Gilbey of the Cooling the Tube Programme; the authors are grateful for their input.

- © 2010 Geological Society of London