Deep geothermal resource assessment of early carboniferous limestones for Central and Southern Great Britain

.


Introduction
Geothermal energy has the potential to provide a stable supply of clean heat energy for Great Britain, helping to reduce carbon emissions, diversify from fossil fuels and improve domestic energy security.Exploration for geothermal resources has traditionally focused on development of high-temperature geothermal systems for power generation, usually associated with large igneous intrusions or areas with high volcanic or tectonic activity (e.g., Iceland and New Zealand).Identification of geothermal plays in conductive intracratonic basins, (defined as play type CD-1 by Moeck (2014)) as in the case of the UK and western Europe, has demonstrated growing interest in the use of low--medium temperature geothermal systems (50-150 • C) as an alternative and reliable source of heat for direct-use applications.Such uses can range from district heating schemes to sustainable agriculture and food production (IRENA, 2022), as well as bathing and swimming.Technological advancements may also enable the utilization of low to medium temperature resources for electricity generation through the use of binary power plant technologies (Stober, 2021).
Low to medium temperature geothermal resources in the UK were previously assessed as part of a national 'Geothermal Energy Programme' conducted between 1977 and1994 which targeted deep onshore Mesozoic sedimentary basins (>3 km depth) containing Permo-Triassic and Palaeozoic sandstones (Downing and Gray, 1986).An initial geothermal resource assessment estimated that the deep sandstone aquifers within these basins could contain up to 328 × 10 18 J or 328 EJ of heat in place (Busby, 2014).
The UK bedrock also contains large volumes of (often deeply buried) Early Carboniferous limestones (Fig. 1), which, relative to Permo-Triassic sandstone aquifers, have been somewhat overlooked as a source of low-medium temperature geothermal energy (Downing and Gray,1986).In mainland Europe, equivalent limestones have been developed with working hydrothermal systems in both Belgium and the Netherlands (Bos and Laenen, 2017;Mijnlieff, 2020;Reijmer et al., 2017).The Jurassic Malm limestone in the Molasse Basin of SW Germany although not contemporaneous in geological age is a proven analogue of how fractured limestones can be used to generate heat using hydrothermal well doublet systems (Agemar and Tribbensee, 2018;Agemar et al., 2014;Stober, 2014).The most known example of an exploited deep geothermal carbonate reservoir is from the Dogger limestone located in the central part of the Paris Basin, in Ile-de-France.This reservoir has been well-developed since the 1960s and consists of oolitic limestone and dolomites of Mid-Jurassic age, at depths and temperatures ranging from 1450 to 200 m and 56 to 80 • C respectively (Lenoir, 1992 andLopez et al., 2010).The UK Early Carboniferous limestones (ECL) form a productive aquifer (Abesser et al., 2005;Abesser and Smedley, 2008), and host several deep hydrothermal circulation systems which are associated with thermal springs at Bath, Bristol, the Taff Valley (South Wales) and the Peak District (Downing and Gray, 1986).Elsewhere, it presents a potentially important geothermal target, as the geographical extent of the ECL is coincident with populous regions of the UK, where domestic heat energy demand is high.Unlike sandstone aquifers, permeability in these carbonates is highly heterogenous.Fluid flow is thought to largely concentrate along faults, fractures, and bedding planes, and is significantly influenced by karstification and alteration.Successful development of geothermal resources will therefore require the identification of areas with suitable lithology at sufficient burial depth and temperature, that are likely to be affected by extensive faults, fracture networks and karst development.
Regional heat-in-place (HIP) resources assessments of the ECL have been calculated for the Netherlands/Belgium and Ireland as part of the GeoERA HotLime project, co-financed from the European Union's Horizon 2020 (https://geoera.eu/projects/hotlime6/.A preliminary assessment of the UK ECL geothermal potential was made by Narayan et al. (2021).This study extends this work by characterising the depth, thickness, and extent of the ECL to develop new volumetric geological models of the onshore ECL for central and southern England.These models are used as input in to established volumetric and statistical workflows (Garg and Combs, 2015;Muffler and Cataldi, 1978;Piris et al., 2021) to estimate the heat-in-place and recoverable heat resources for central and southern England.

Regional geology
In Great Britain, Mississippian (early Carboniferous mostly Tournaisian-Visean) platform and reef-limestones are thick enough for geothermal exploitation (> 50 m thick) are present in two provinces (Fig. 1 and 2), south and north of the Wales-Anglo-Brabant Massif (WABM).A detailed geological description of the two provinces can be found in Pharaoh et al. (2021).In this paper we refer to these simply as the northern and southern province.
In the northern province, there is a polygonal mosaic of deep-water basins, ramps and platforms developed in response to extensional stress in early Mississippian times (Fraser and Gawthorpe, 1990).The interaction of local tectonic controls, e.g., reactivation of Caledonide basement structures, with glaciation enhanced eustatic sea level change, led to the development of complex carbonate system tracts on these ramps and platforms (Davies et al., 2008).The northern basins lay at a greater distance from the developing front of the Variscan Orogen and therefore were less strongly deformed than the southern shelf.The dominant tectonic response to orogenic compression was inversion, particularly along favourably orientated structures (Corfield et al., 1996).
In the northern province, depth of the ECL ranges from surface outcrop to >5500 m sub-ordnance datum (Pharaoh et al., 2021).Areas where ECL platforms are likely to be more prospective for geothermal heat resources, and in many cases, underlie major conurbations with high heat demand (e.g.Greater Manchester).Burial beneath thick Permian-Mesozoic sequences, as in the Cheshire Basin, and the flanks of the East Irish Sea and Southern North Sea basins; or beneath thick later Carboniferous strata, as in the Stoke-on-Trent area, are regarded as potential enhancing factors.Several basins containing generally thin developments of ECL occupy the northern edge of the WABM.These include the Stafford, Needwood and Knowle basins in the west; and the Hathern Shelf, marginal to the East Midlands basin complex, in the east.Throughout the province, platforms are juxtaposed with deep basins of the same age, in which carbonate facies may be partially or completely absent and clastic facies are present (Fraser and Gawthorpe, 1990).
In the southern province, the ECL were deposited in Mississippian times upon a shelf occupying the northern margin of the Rhenohercynian Basin, and extending from the Rhineland via Belgium, across southern England, to Ireland and likely, beyond.This shelf sloped into deep water in the south and was subsequently incorporated into the Rhenohercynian (externide) Zone of the Variscan Orogen in late Carboniferous time (Pharaoh et al., 2021).Consequently, it experienced significant thrusting, folding and post-Carboniferous erosion.
ECL outcrops around margins of the coalfield basins in South Wales, the Forest of Dean, Bristol-Somerset, and in the Mendip Hills (Fig. 1).This region is notable for the development of thermal springs, of which the thermal aquifer at Bath which have temperature of 45-46 • C is perhaps the best-known example.The presence of palaeo-karstification which formed in the Triassic in these limestones likely facilitates the passage of water to the Bath occurrence (Gallois, 2007).Unfortunately, few of the existing geophysical 2D seismic reflection data, optimised to image the hydrocarbon-prospective Mesozoic section of the Wessex Basin provide a clear image of the Carboniferous structure and there are limited borehole penetrations (Pharaoh et al., 2021).Here the ECL is buried at depths up to 2.5 km and form a relatively continuous E-W trending belt comprising fold-thrust nappes, comparable to those of the Dinant Syncline of Belgium (Busby and Smith, 2001;Pullan and Butler, 2018), close to the Variscan Orogenic Front (Fig. 1).Underlying Devonian strata, perhaps as much as 2000 m thick are thought to exist but have low porosity and permeability and therefore not considered prospective for geothermal exploration (Downing and Gray, 1986).

Datasets
Key subsurface datasets used in this study included legacy 2D seismic reflection data, borehole stratigraphic and lithological information, and borehole bottom hole temperature profiles.

Seismic data
The Carboniferous Basins of England have been extensively explored for hydrocarbons and are covered by a grid of 2D seismic surveys, mostly acquired 1970-1990 (Fig. 3).Over 24,000 line-kilometres of 2D seismic data are available from the UK Onshore Geophysical Library (UKOGL) across the northern province study area.Across the southern province, 20,700-line kilometres of seismic data are available from the same source.
A large library of existing BGS (British Geological Survey) legacy interpretations based on this data were also available to this study, originating from a diverse range of projects that include production of the BGS Subsurface Memoirs (Chadwick et al., 1995;Pharaoh et al., 2011;Plant et al., 1999;Smith et al., 2005), regional geological mapping activities (Fraser and Gawthorpe, 2003), and assessments of shale gas potential across England (Andrews, 2013).This study repurposes and builds on this database for the purpose of characterising the depth, thickness, and extent of the potential Early Carboniferous deep geothermal reservoir.Interpretation data was compiled in Schlumberger's Petrel™ software package, and the interpretation extended as needed to fill any significant gaps in coverage.
Interpretation data for the Base Permian, Top Tournaisian-Visean (top ECL), and Base Carboniferous horizons were converted from twoway travel time into depth to supply inputs for each 3D geological model.Velocity data in the form of checkshots and vertical seismic profiles (VSP) were used to derive velocity functions for each study area.For the northern province, a two-layer velocity model was used, based on separate V 0 +k velocity functions for the Mesozoic and Palaeozoic intervals.Values for the initial velocity (V 0 ) were interpolated between wells, and a constant velocity gradient (k) was applied for each interval.This approach was chosen to account for the highly variable thickness of siliciclastic sediments preserved within the high-relief Carboniferous block and basin system.In the southern province, Upper Carboniferous siliciclastic sediments are absent, and the Base Permian Unconformity therefore coincides with the top of the Tournaisian-Visean ECL limestones.This allowed for a simpler approach, whereby the average velocity to the Base Permian Unconformity surface was calculated at each well, and then interpolated to generate an average velocity map for the Mesozoic interval.This map was then convolved with the interpretation data in two-way time to generate a depth map for the Top Tournaisian-Visean limestones.

Borehole and well data
The UK Carboniferous is penetrated by numerous boreholes related to a long history of exploration for coal, groundwater, and hydrocarbons (Fig. 3).Formation tops were compiled from released composite well logs, BGS borehole records, and the UK Geothermal Catalogue.Lithofacies-scale information was primarily derived from geophysical log parameters and checked against cuttings returns and core where available (Fig. 4).
Within the northern ECL province study area, 241 boreholes have been identified as encountering strata of Tournaisian-Visean age (Fig. 3).Of these, 72 prove a continuous succession of limestones more than 50 m thick before terminating within the formation.A further 39 boreholes encounter >50 m of mixed siliciclastic rocks of Tournaisian-Visean age, which represent the equivalent basinal (non-prospective) facies.Only 24 boreholes record the full thickness of Tournaisian-Visean age rocks, terminating in Pre-Carboniferous or igneous rocks.
Across the southern ECL province, the seismic interpretation is constrained by 129 boreholes that reach Carboniferous rocks.The majority of these are hydrocarbon exploration wells that targeted Mesozoic formations, and most do not prove the full thickness of the Carboniferous Limestone Supergroup.The Cannington Park borehole in the Somerset Basin proves the thickest Carboniferous limestones at around 600 m.Elsewhere the ECL has been deeply eroded beneath the Variscan unconformity (Pharaoh et al., 2021).

Fig. 2. A generalised stratigraphic chart
showing the lithologies, facies and the marker horizons for the 3D geological model across both northern and southern provinces (adapted from (Waters et al., 2007a;Waters et al., 2007b).Epigraphs shown are labelled as important geological ages for these provinces; Q -Quaternary, N -Neogene, Pl -Pliocene, K -Cretaceous, J -Jurassic, Tr -Triassic, P -Permian, C -Carboniferous, D -Devonian.The colours of the of the stratigraphic horizons reflect the colours used in the 3D model.

Temperature database
The BGS UK Geothermal Catalogue provides a database of subsurface temperatures compiled from over 1300 onshore boreholes, together with measurements of thermal conductivity and heat flow.The database was originally compiled in the early 1980 ′ s from a series of BGS reports, and subsequently updated through a series of several revisions, most recently by Rollin (1987).The current version contains additional well data from the 1990 ′ s and 2000 ′ s, including from wells at Sellafield Nuclear Power Plant and the Newcastle Science Centre geothermal well.
Temperature data held in the database have been compiled from a variety of measurement types, which include bottom hole temperatures (BHT), wireline temperature logs, and coal field temperature measurements.BHT measurements are the most abundant measurement type in the catalogue and are typically recorded several hours to several days after drilling is completed.Due to the cooling effect of the circulating mud, BHT values tend to underestimate the host rock temperature, and require a correction value to be applied according to the amount of elapsed time between the cessation of mud circulation and the temperature measurement (the 'elapsed time').Corrected temperature values are available in the catalogue based on a set of empirically derived correction curves which are valid within defined temperature, depth, and elapsed time ranges.Where measurements fall outside of these ranges, corrected temperature values are not available.The correction curves are most reliable where the elapsed time is long, with uncertainty ranging from ±5.5 • C after 6 h, to ±1 • C after 20 h.A full description of the correction methodology is provided in Burley et al. (1984).For this study, a subset of wells located within each study area were extracted from the catalogue.Corrected temperature values are available for 209 wells in the Northern England area, and for 120 wells in the Southern England area.As direct temperature measurements of Carboniferous rocks are derived almost exclusively from BHT records, these values have some uncertainty associated with both the measurement equipment and the correction method applied.Estimates of the geothermal gradients are based on linear regression analysis of these wells for each study area.

Creating the 3D geological grid
Separate 3D geological models were created for the northern and southern ECL province using borehole horizon markers, depth converted seismic picks and polygons defining the erosional limits of formations at outcrop from BGS surface mapping (Fig. 5).They were created in SKUA-GOCAD using the Structure & Stratigraphy workflow.A digital terrain model (OSTerrain50 -Ordnance Survey, 2022) formed the upper surface of the models.Due to the complexity of the geology and the variable spatial distribution of the control points, an implicit geological modelling method (e.g.Cowan et al., 2003) was used where all geological units were modelled simultaneously using all the available data held and interpolated within a 3D framework.Rules were applied to ensure that stratigraphic relationships such as onlap and truncation at unconformities were honoured.A simplified fault network was used by selecting the largest structures where interpreted as necessary.The result of this process was a 3D stratigraphical (irregular) grid discretised into 'regions' which correspond to broad Stratigraphical units (formations or groups) and/or approximate chrono-stratigraphical time slices (Fig. 5).This grid then attributed with additional properties such as lithology and density that were needed for the analysis.

Lithological attribution of the geological models
Where geological regions within the 3D gridded model were broadly homogeneous, they were assigned to a single lithological category (e.g., limestone, Table 1).Where geological regions included significant lithological variation, they were further refined to include this heterogeneity using a stochastic lithological modelling process.Control data for this process was derived from a new lithological classification with integration of lithological interpretation of composite and geophysical logs, of which gamma ray and sonic transit time were the most widely available.To minimise any inconsistencies that arise from differences in tool calibration or downhole conditions (e.g., composition of the drilling Fig. 4. Selected wells show lithological variability in the Tournaisian-Visean strata.Geophysical logs were converted into a lithological classification using a cut-off analysis method.fluid), gamma ray and sonic logs were first normalised to a common reference range using standard techniques (Shier, 2004).The result was then compared to the chippings logs to access their accuracy.
Geophysical logs were converted into a lithological classification using a cut-off analysis method (Fig. 4).Eight broad lithological types were create using gamma-ray and sonic log parameters.Whenever possible, results were verified against borehole core and cuttings returns.
Optimal results were achieved by applying different combinations of cut-off criteria according to the age of the strata: Westphalian and Namurian post-rift strata are dominated by prograding clastic deltaic environments that contain few limestone beds, and low gamma ray values tend to reflect clean sandstone units (the 'grits' of the Millstone Grit Group).In the Visean-Tournaisian succession thick sandstones are rare, and a low gamma ray response is more typical of carbonate rocks.Several boreholes are known to contain thick volcanic tuff deposits within the Namurian succession.Gamma ray and sonic logs are not sufficient to distinguish these tuff deposits from low-gamma sandstone/ grit units, and instead the tuff deposits were manually corrected based on existing BGS interpretations and company composite logs.
Identifying the location of the limestone at depth is critical as only they have known high secondary permeabilities.While there are many different paleogeographic basin and block reconstructions of ECL (e.g.Smith et al., 2005, Andrews 2013, Pharaoh et al., 2021) none of these show how the limestone vary vertically.To understand this we use a stochastic method in this study.
The lithology logs computed for each borehole were transferred to the 3D geological grid, using upscaling criteria of most frequently occurring lithology in a voxel where required.The lithology data were then simulated throughout the geological grid using a stochastic geostatistical technique (Sequential Indicator Simulation) available in GOCAD-SKUA.This was carefully controlled by modelling within specified grid regions (e.g., a particular formation), providing constraints on the expected number and relative proportions of different lithologies, and replicating the known vertical and lateral spatial distribution within the formation.Experimental variograms to control the lithologies were created.These were calculated from the well data; one for the Limestones; Sandstones and Conglomerates and one for the Mudstones.Separate sets of variograms were calculated for the Tournaisian-Visean and Namurian intervals to reflect the changes in facies and provenance between these unconformities creating a total of 6 variograms.The simulations were run one hundred times and the most frequently occurring lithology was used.This was used because the heat in place calculator needs each voxel to be attributed with a single lithology.

Density attribution of the geological models
The gridded geological model was attributed with appropriate bulk density values (Table 1) based on average values that are known from the two areas (Rollin, 1987).These were mapped to the 3D geological grid using a combination of lithostratigraphic unit and any internal lithological variation that had been modelled.The attributed geological grids were then converted to regular grids with a 2500 × 2500 metre horizontal resolution and 50 metre vertical resolution to be imported into the volumetric 'Heat in Place' (HIP) calculator.

Computing an average geothermal gradient
An average geothermal gradient for the northern and southern provinces was calculated by plotting depth against corrected bottom hole temperature for 199 wells in the northern province and 117 wells in southern province, using temperature data from the UK Geothermal Catalogue (Rollin, 1987).Mean annual surface temperatures vary slightly between study areas due to differences in latitude and elevation.Modal values of 10.1 • C and 10.9 • C were therefore used for the northern and southern provinces respectively.Seventeen wells with anomalously high or low gradients were excluded from the analysis.This method of analysis reveals a clear temperature-depth relationship, with geothermal gradients of 28.7 • C km − 1 and 31.3 • C km − 1 derived for the northern and southern provinces respectively (Fig. 6).

Resource calculation
A geothermal resource assessment for the UK ECL was calculated using the '3DHIP Calculator' tool developed by Piris et al. (2021) and freely accessible for download from the website of the Cartographic and Geological Institute of Catalonia (ICGC, 2021).This tool is based on the USGS volumetric 'Heat in Place' method, developed by Muffler and Cataldi (1978) and later revised by various authors (Allen et al., 1997;Garg and Combs, 2015;Williams, 2007;Williams et al., 2008).The tool uses a Monte Carlo simulation to calculate a probability distribution function for the heat in place (HIP) and recoverable heat fraction (Hrec) for a supplied geothermal reservoir volume.The volumetric HIP method calculates the heat energy (in joules) stored in the both the rock mass and the formation fluid and remains the most common method for estimating resource potential in deep geothermal reservoirs.HIP is given by Eq. ( 1): where V is the reservoir volume (m 3 ), ϕ is the porosity, ρ F is the fluid density and ρ R is rock density (kg/m 3 ), C F is the fluid specific heat capacity and C R is the rock specific heat capacity (kJ/kg.• C), T r is the reservoir temperature ( • C), and T i is the re-injection temperature ( • C).The Recoverable Heat (Eq.( 2)) provides an estimate of the producible thermal power (in kilowatts), based on assumptions regarding the conversion efficiency of the heat exchanger (C e ), a recovery factor (R), the expected lifetime of a geothermal project (Tlive), and the proportion of time a plant is likely to be operating (plant factor, Pf).
Reservoir parameters are input as probability distribution functions and are listed in Table 2. Geothermal gradient and mean surface temperature were defined as constants for each model (see Section 4.5).Reservoir volume and temperature are derived from the 3D geological models for each study area, which represent the P50 distribution of ECL lithofacies.P50 is referred to as "probable", whereas P10 is "proved" and P90 as "possible".P10, P50 and P90 correspond to the 10th, 50th and 90th percentiles of the calculated cumulative distribution function.A total of 10,000 simulations were run for each model using an XYZ spatial resolution of 2500 × 2500 × 50 m.An upper depth cut-off of 1000 m and 1200 m was applied to northern and southern models respectively, which corresponds to a reservoir temperature of 50 • C and is judged to be the minimum temperature required for direct-use applications of geothermal energy.The resulting heat-in-place values are an estimate of the total geothermal energy contained within the volume of the ECL, below 1000 m depth, and represent the Accessible Resource Base (sensu Muffler and Cataldi, 1978).

Results
The probability distribution outcomes for total HIP and recoverable heat are summarised in Table 3.A breakdown of resources by combined local authority area, along with predicted ECL depth, thickness, and temperature ranges, are presented in Table 4 Maps of the depth, temperature and distribution of HIP are shown in Fig. 7 and Fig. 8.
The limestone lithology of the ECL is present across much of the northern ECL province, generally forming in shallow marine basins and across intra-basinal highs (Fig. 1).Limestones are well developed beneath parts of the East Midlands at a depth of around 1200 m, reaching up to 1600 m in thickness (Fig. 7).Limestones are also present across large parts of Cheshire and Merseyside, Greater Manchester, West Yorkshire and Harrogate, at depths ranging from less than 400 m to greater than 3000 m.Here, vertical thickness of the ECL reaches a maximum of around 2750 m in the Huddersfield Basin, adjacent to the Lymm High (see Fig. 1 for names of Early Carboniferous highs and basins).Limestones are absent/thinly developed (<50 m) within deep Carboniferous basins such as the Bowland Basin, Gainsborough and Edale Half-grabens and the Widmerpool Gulf, which underlie much of Lancashire and South Cumbria, South Yorkshire and Bassetlaw and the western edge of the East Midlands.Across parts of the West Midlands, to the south of the Hathern Shelf, well data show a high proportion of siliciclastic sediments within the Tournaisian-Visean succession, and limestones are thought to be comparatively thin.Thick limestone deposits have not been proven within the Cheshire Basin and are assumed to be absent for this study.
In the southern ECL province, the ECL is present at subcrop beneath the Variscan Unconformity, and is preserved mainly within a series of discrete, east-west trending sub-basins (Fig. 1).In the north-west of the province, the ECL is found at the surface forming the Mendip Hills and Wye Valley in Somerset, Gloucestershire and South Wales.Burial depth increases to greater than 1000 m across parts of Wessex, Surrey and Sussex, and reaches a maximum of 2500 m beneath Winchester and Crawley, before shallowing again to less than 500 m along the Kent and Medway coast (Fig. 7).Across much of the south-east of England, preserved thickness of the ECL is relatively thin, ranging from 150 to 350 m.In the Mendip Hills, where true stratigraphic thickness is measured to be approximately 800 m, the strata are folded into deep fold-thrust nappes, and the vertical thickness of the ECL exceeds 2500 m; this is reflected in the HIP calculation.
The calculated geothermal resource is strongly dependant on, and proportional to the depth and thickness of the predicted ECL reservoir, and this is reflected in both the predicted temperature at top ECL, and in the spatial distribution of total HIP (Fig. 8).In the northern ECL province, temperature trends reflect the depth to top ECL, ranging from approximately 50 • C across the East Midlands, and West Yorkshire and Harrogate, to over 100 • C beneath Greater Manchester and west Cheshire.Total HIP more closely reflects the thickness trend, with values of over 200-350 PJ/km 2 in areas with the greatest ECL thickness, near to Lincoln and north-west of Manchester, contrasting with a background value of 40-100 PJ/km 2 .In the southern ECL province, temperature at top ECL is greatest in the centre of the basin, reaching 85-95 • C near to Winchester and Crawley, and up to 80 • C south of Bath.As a result of the reduced thickness of ECL across much of the southern ECL province, high values of HIP are limited to the fold-thrust nappes of Somerset and Wiltshire, where values exceed 250 PJ/km 2 .

Discussion
The mean annual demand for domestic heating is shown in Fig. 9. Comparison with the calculated total HIP highlights several areas where high HIP density corresponds with high heat demand.In the northern ECL province this includes parts of Merseyside, Greater Manchester, West Yorkshire and Harrogate.In the southern ECL province, the areas around Bath and Bristol have high potential.Table 4 shows the total HIP and estimated recoverable heat (Hrec) by combined local authority area.
A tentative estimate of recoverable heat (Hrec), based on an assumed project duration of 30 years, yields a value of 133 GW for the northern province, and 19.1 GW for the southern province.This estimate is based on various scalar factors to account for recovery and production losses, but still considers the full volume of the ECL reservoir.In practice, successful exploitation of heat from the ECL is likely to be limited to areas with remarkably high flow rates.Such flow rates may be associated with fault damage zones and/or areas affected by palaeo-karsts.Williams (2007) showed that fracture spacings of less than 30 m may be required to effectively extract heat from the rock mass; with larger fracture spacings, the circulating fluids increasingly bypass the rock mass and less thermal energy is extracted.An improvement on the current method would be to consider only the volume of rock predicted to be associated with sufficient fracturing and fracture connectivity.Further work is needed to calculate the volume of rock associated with the damage zones of major faults, and to consider their flow properties in relation to the present-day stress regime.This assessment would also need to consider the formation and distribution of palaeo-karsts and whether these are likely to provide enhanced flow at depth.Such an assessment would yield far smaller values of recoverable heat.Based on known large scale regional understanding areas of intense faulting and ideally where Variscan and post-Variscan faulting intersect are likely to produce enhanced areas of transmissivity and flow rates sufficient for geothermal production.

Uncertainty and accuracy of predictions
The sensitivity of the HIP calculation to reservoir depth and thickness is evident when comparing the maps presented in Fig. 7 and Fig. 8: high values of HIP north-east of Manchester and on the East Midlands Platform near to Lincoln are related to the particularly thick vertical sections of ECL in these areas.Conversely, areas in the southern province such as Winchester and Crawley, where elevated temperatures of between 50 and 100 • C are predicted at top ECL but reservoir thickness is low, show correspondingly low values of total HIP.The difficulties involved in distinguishing the extent of limestone deposits at depth and in identifying the base of the ECL succession directly from seismic data mean that the estimation of total reservoir thickness remains the main source of uncertainty for this study.Thickness of the ECL is influenced by both the depth to the top and base of the Tournaisian-Visean interval (base Carboniferous), and the stochastic simulation of limestone lithofacies within this interval.Depth to the top Tournaisian-Visean interval is well-constrained by seismic interpretation data and well picks.The depth error in the final model was assessed by calculating the mismatch between the final models and the available well picks.This indicated a median error of ±47 m for the northern ECL province and ± 126 m for the southern ECL province on the top Tournaisian-Visean (and higher on the base).This error is highest in areas with fewer wells.Depth to the base of the interval is constrained by interpretation data across approximately 80% of the study area but is encountered in only 14 wells, so the error will be larger for that horizon.No data is available to constrain the base of the interval across the northeast of the study area (indicated in Fig. 7), which includes much of the Market Weighton High (Humber Coast and Vale area); here the ECL is encountered proving up to 90 m of limestone before terminating within the limestones.To limit over-estimation of ECL in this area, the maximum ECL thickness has been constrained to 490 m, reflecting the average thickness of ECL proven by wells across similar platform areas.
As a result, resource estimates in this area are poorly constrained and should be treated with caution.Uncertainty associated with the location of the limestone lithofacies also affects the total reservoir thickness and is greatest where wells with available geophysical logs are widely spaced, or where limestones form a minority component of the succession.In the West Midlands area, six wells have encountered the Tournasian-Visean interval and prove a  succession containing mixed carbonates and siliciclastic sediments.Limestones make up an appreciable part of the succession but are characterised by relatively thin beds (<15 m) separated by thick intervals of siliciclastic sediments.Due to the limited vertical resolution of the geological model, these areas are classified as siliciclastic dominated.The predicted ECL across this area (indicated in Fig. 7 and Fig. 8) is represented as discontinuous and thin and may under-represent the volume of limestone (especially the occurrence of reef knolls), and consequently the total geothermal resource in this area.

Comparison with other resource assessments of the carboniferous limestones
Various approaches for calculating heat resources have been done for assessing the geothermal potential of the ECL in the Netherlands and Ireland (Veldkamp et al., 2021).The Netherlands assessment, which was a country-wide assessment, estimating the recoverable heat of 3183 EJ.
Other recent estimates for the UK include an assessment by Narayan et al. (2021), which yielded a resource estimate for HIP of 129 EJ (P50 case) for the UK ECL.This assessment used a fundamentally different approach based on constant (averaged) values for reservoir depth and thickness for each regional area, and considered only the fluid fraction of the reservoir, not including the heat stored in the rock mass.Additional differences in the method Narayan et al. (2021) used is the surface temperatures as the reinjection temperatures whereas this study used 21 • C as stated in Table 2 taken from Limberger et al. (2018).

Conclusions
The Early Carboniferous limestones may offer significant potential as a resource for deep geothermal energy across large parts of central and southern Britain, where porosity/permeability is enhanced by faulting/ fracturing and karstification.In this assessment we integrate seismic and well subsurface data to model the 3D geometry of limestones and calculate a total HIP geothermal resource of 1415-1528 EJ for the limestone succession.For an assumed 30-year extraction period, a tentative maximum value of 106-222 GW of heat potential for direct use applications is calculated, though it should be noted that only a small percentage is likely to be recoverable in practice, in zones of enhanced permeability.Based on geothermal gradients calculated here, direct use of geothermal energy for heating could be possible where the ECL reaches depths greater than 2 km and power generation may be feasible at depths greater than 4 km.The largest HIP resources are found in Central Britain with greatest potential across parts of the East Midlands and Greater Manchester.In southern Britain, the largest resources are in the Somerset, Wiltshire, Avon and Gloucester regions.HIP calculations are sensitive to the estimated thickness of the ECL reservoir, and this is a significant uncertainty in the geological modelling, which relies heavily on interpretation of legacy seismic data with relatively few well data penetrating the base ECL and limited associated time-depth data points for depth conversion.Reducing this uncertainty will require investment to acquire new geophysical datasets, particularly over the Humber, Coast and Vale, and West Midlands regions where current data are sparse or of inadequate quality.Further work is required to assess the distribution of permeability in fractured carbonate reservoirs, and to identify areas where fracture networks and palaeo-karst are likely to provide enhanced flow characteristics.These conditions may be associated with the damage zones of major faults.

Fig. 3 .
Fig. 3. Seismic and well database for the northern and southern areas of interest.For wells that reach the Visean-Tournaisian succession and that have sufficient data, an assessment of the dominant lithofacies is indicated (by symbol colour).

Fig. 6 .
Fig. 6.Corrected bottom hole temperatures for 199 wells within the northern province and 117 wells within the southern province.Different symbols indicate the age of the strata that the temperature measurements are associated with geothermal gradients of 28.7 • C km-1 and 31.3 • C km-1 are derived for the northern and southern provinces, respectively, with R 2 values of 0.87 and 0.83.Based on mean annual surface temperatures of 10.1 • C for the northern province and 10.9 • C for the southern province.

Fig. 7 .
Fig. 7. Maps of the northern and southern provinces showing (a) depth to Top Early Carboniferous limestones and (b) thickness of Early Carboniferous limestones, overlain with ECL outcrop, major Variscan faults and combined local authority areas.

Fig. 8 .
Fig. 8. (a) Map showing the temperature distribution at top of the Early Carboniferous limestones (b) Map showing the distribution of HIP (PJ/km 2 ) for the ECL, based on the P50 resource outcome.Areas of high uncertainty areas are indicated.

Fig. 9 .
Fig. 9. Map showing the annual domestic heat demand for 2009, after Taylor et al. (2014).Contains data made available by the School of Civil and Building Engineering, Loughborough University.

Table 1
Lithostratigraphy, bulk lithology and bulk density components of the 3D models created for the Northern and Southern Provinces.

Table 2
Petrophysical and operational parameters used for heat-in-place and recoverable heat calculations.

Table 3
Total heat-in-place (PJ) and recoverable heat (kW) for the UK Early Carboniferous limestones (ECL).

Table 4
Depth, thickness, temperature, HIP, tentative recoverable heat, and domestic annual heat demand for each combined local authority area.