Seismicity at Newdigate, Surrey, during 2018–2019: A Candidate Mechanism Indicating Causation by Nearby Oil Production

Abstract

During 2018–2019, oil was intermittently produced from the Late Jurassic Upper Portland Sandstone in the Weald Basin, southeast England, via the Horse Hill-1 and Brockham-X2Y wells. Concurrently, a sequence of earthquakes of magnitude ≤3.25 occurred near Newdigate, ∼3 km and ∼8 km from these wells. The pattern, with earthquakes concentrated during production from this Portland reservoir, suggests a cause-and-effect connection. It is proposed that this seismicity occurred on a patch of fault transecting permeable Dinantian limestone, beneath the Jurassic succession of the Weald Basin, hydraulically connected to this reservoir via this permeable fault and the permeable calcite ‘beef’ fabric within the Portland sandstone; oil production depressurizes this reservoir and draws groundwater from the limestone, compacting it and ‘unclamping’ the fault, reaching the Mohr-Coulomb failure criterion and causing seismicity. In principle this model is fully testable, but required data, notably the history of pressure variations in the wells, are not currently in the public domain. Quantitative estimates are, nonetheless, made of the magnitudes of the variations, arising from production from each well, in the state of stress on the seismogenic Newdigate fault. The general principles of this model, including the incorporation of poroelastic effects and effects of fault asperities into Mohr-Coulomb failure calculations, may inform understanding of anthropogenic seismicity in other settings.

1. Introduction

Highlights.

Earthquakes at Newdigate in 2018–2019 correlate with oil production from Portland sst.

A cause-and-effect relation via a high-permeability hydraulic connection is proposed.

The model seismogenic fault incorporates internal structure with asperities.

Testing by Mohr-Coulomb failure analysis incorporates effects of poroelasticity.

A ‘swarm’ of earthquakes with magnitudes up to ∼3, starting on 1 April 2018, has affected the Newdigate area of Surrey, in the Weald Basin of southeast England (Figures 1–3). As is detailed in the online supplement, a potential connection with local oilfield activities, in the nearby Brockham-X2Y (BRX2Y) and Horse Hill-1 (HH1) wells, was immediately suspected, but dismissed by petroleum developers (e.g., [2, 3]). Concerns about the possibility that activities in these oilfields were indeed causing these earthquakes were raised through correspondence in The Times newspaper in August 2018 [4]. A workshop, convened by the Oil & Gas Authority (OGA), followed on 3 October 2018, the OGA being a UK government body with responsibilities that include the licencing of exploration and development of onshore oil and gas resources in England, including managing the risk of seismicity from such operations. A summary of the proceedings of this workshop was reported [5], including the statement that ‘the workshop participants concluded that, based on the evidence presented, there was no causal link between the seismic events and oil and gas activity although one participant was less certain and felt that this could only be concluded on “the balance of probabilities” and would have liked to see more detailed data on recent oil and gas surface and subsurface activity’ ([5], p. 1). It has subsequently been argued that there is indeed no such cause and effect connection (e.g., [1, 6]), developers having repeatedly issued strong public statements to this effect (e.g., [3, 7, 8]). However, a major issue, not noted in any of the above-mentioned works, is the clear temporal pattern of earthquake occurrence (Figure 4), with earthquakes strongly concentrated at times when oil is being produced from the Upper Portland Sandstone via the HH1 and/or BRX2Y wells. Production will reduce the fluid pressure in the reservoir being pumped. Fluid pressure changes within faults are well known as a cause of anthropogenic seismicity (e.g., [9, 10]); however, rather than a decrease, the causative change is usually an increase in fluid pressure as, for example, for the Preese Hall earthquake sequence in 2011, caused by injection of water under pressure during ‘fracking’ for shale gas (e.g., [11]).

Figure 1.

Map of the study area, modified from Figure 2(a) of Hicks et al. [1]. The original geographical (latitude-longitude) co-ordinate system has been retained, but ‘greyed out’, with a new co-ordinate system added, indexed to the British National Grid (BNG). Inset shows location. As is discussed in the main text, the original scale bar by Hicks et al. [1] is much too small and has also been ‘greyed out’. Faults are identified thus: BF, Brockham fault; BHF, Box Hill fault; CF, Crawley fault; COF, ‘Collendean fault’; FGF, ‘Faygate fault’; HF, Holmbush fault; HHF, Horse Hill fault; HKF, Hookwood fault; HWF, Holmwood fault; KFF, Kingsfold fault; LHF, Leigh fault; NGF, Newdigate fault; OKF, Ockley fault; WB1F, and Whiteberry-1 fault. Most of these structures are depicted as shown by Hicks et al. [1], although some are misplaced or misidentified, as discussed in the text. The Crawley and Holmwood faults, not recognised by Hicks et al. [1], are shown schematically (from [15]) where they cross seismic line TWLD-90-15, the southward continuation of which (beyond the excerpt shown in Figure 2) is also shown schematically. The ‘Faygate fault’ is a mistaken concept by Hicks et al. [1], so is shown ‘greyed out’ (see text and online supplement). Horse Hill 1 well track is from https://ukogl.org.uk/map/php/well_deviation_survey.php?wellId=3041. Positions of seismic lines, including line TWLD-90-15, are from the schematic location map provided by the UK onshore geophysical library (https://ukogl.org.uk/), which is indexed to the BNG, and was transformed to geographical co-ordinates by Hicks et al. [1]. Seismograph station GATW ceased operation on 17 May 2019 due to equipment theft. It was replaced by station GAT2, ∼230 m northwest, from 6 June.

Figure 2.

2-D seismic section along seismic line TWLD-90-15, modified from Figure 6 of Hicks et al. [1]. The original was indexed by Hicks et al. [1] to the British National Grid (BNG), rather than the geographical co-ordinates used for Figure 1; it has been geolocated for this study using documentation provided by the UKOGL. The labelled horizons were not explained by Hicks et al. [1], but appear to be the top Portland group, top Kimmeridge clay formation, and top coralline Oolite formation (cf. Table 4). Faults designated by Hicks et al. [1] are identified thus: COF, Collendean fault; LHF, Leigh fault; and NGF, Newdigate fault. CF denotes the Crawley fault. The depth scale from Hicks et al. [1], which appears to be based on their velocity model for earthquake location (Table 2), is ‘greyed out’, being now considered inaccurate. The new depth scale, using the seismic interval velocities from well HH1 (Table 1) and now considered more accurate, is emphasised. Additional interpretation has also been added, including the interpreted top Penarth group / base Lias group reflector and its offset by the main strand of the Newdigate fault, and some of the additional lesser fault strands forming the multi-stranded Newdigate fault zones, other strands being evident in Figure 3 and in the uninterpreted version of this seismic section provided by Hicks et al. [1] in their online supplement.

Figure 3.

Excerpt from the record section for seismic line TWLD-90-15 between shot points 1400 and 2070, as provided by UKOGL, illustrating the Newdigate fault. As geo-located using documentation provided by UKOGL, this excerpt extends between BNG references TQ 21618 38952 and TQ 20381 46852, a distance of ∼8 km. The location of this seismic line adjoining the Newdigate fault is illustrated in Figures 1 and 7.

Figure 4.

Time series of Newdigate earthquakes and other activities affecting the Horse Hill 1 and Brockham X2 wells, based on Figure 3 of Hicks et al. [1]. Note that their original Fig. 6(d) has been omitted as it depicts an incorrect timeline for production at Brockham up to 2016: The correct timeline is shown in Fig. 11 of Angus [27]. (a) Installation dates of stations forming the local temporary seismic monitoring network. From Fig. 3(a) of Hicks et al. [1]. (b) Detected earthquakes, cumulative number of events, and inferred variations in the completeness threshold magnitude M_C. From Fig. 3(b) of Hicks et al. [1]. Note that some of the magnitudes M_L depicted here differ slightly from those listed in supplementary Fig. S2 of Hicks et al. [1], which feature in discussion in the text. (c) Summary timeline for activities at the Horse Hill 1 well, indicating (based on the information sources discussed in the text and details in part (d)) the phases of production from each reservoir. Notes refer to details discussed in the text, thus, regarding the HH1 well: 1, first known intervention affecting the well, 5 April 2018; 2, removal of bridge plug that had isolated the Portland reservoir from the surface, 4 July 2018; 3, production from KL3; 4, production from KL4; and 5, ‘co-mingled’ production from both KL3 and KL4. Regarding the BRX2Y well, based on the timeline reported by Hicks et al. [1]: 6 denotes the restart of production on 22 march 2018; 7 denotes a later resumption, with injection of water starting on 25 June and (net) production restarting on 28 June (but with both injection and production occurring on 27 June); and 8 denotes the end of production on 15 October 2018. (d) More detailed operations timeline for activities at the Horse Hill 1 well, with flow-period averaged production and cumulative production over time. From Fig. 3(c) of Hicks et al. [1], with further details, including dates, provided in their supplementary Table S4. The information provided by Hicks et al. [1] is much more detailed than that which has been otherwise released into the public domain, and must have been obtained from the developer. However, their reporting of the information does not identify the hydrocarbon reservoirs to which the activities relate (see part (c)), which is essential to reveal the pattern of correlation between seismicity and activities affecting the Portland reservoir (see also the online supplement).

Given the geology of the Weald Basin [12, 13], a conceptual model can be envisaged whereby pressure reduction the Portland reservoir might bring nearby faults to the Mohr-Coulomb condition for slip, as illustrated in Figure 5. As summarised in the figure caption, this model also provides a natural explanation for why production from the deeper Kimmeridge reservoir does not have an equivalent effect. Nonetheless, testing this model is difficult, for several reasons. The map and cross-section reported by Hicks et al. [1] provide the most detailed documentation of the Newdigate seismicity available (Figures 1 and 2), and thus serve as a basis for further discussion. However, a first reason why model testing is difficult is that use of these outputs is problematic because of mistakes in their preparation; before they can be used their geolocation has to be improved (as discussed in the present online supplement, also below). A second reason is uncertainty in the hydraulic properties of elements of the proposed model; this includes the distribution of the permeable ‘calcite “beef”’ fabric within clay-dominated lithologies that are otherwise impermeable. Each of these aspects will be investigated in this study. A third reason why testing the proposed model is difficult is that key operational data, such as pressure variations in oil wells and logs of wellsite activities that might affect reservoir conditions, have been found to be unavailable. Indeed, preparation of this manuscript was delayed pending an attempt to obtain such data from the OGA under UK law using a Freedom of Information (FOI) request; this request was unsuccessful on the basis that the OGA did not hold such data, notwithstanding the scope of their statutory duties. In the absence of pressure data, testing the proposed model will be limited to investigating the magnitudes of pressure variations (and associated changes to the state of stress) that can be anticipated in the model fault and in each well, as a result of operations in the other well, and the time delays for their propagation, subject to the assumed hydraulic properties for each lithology, which are informed by the limited available data.

Figure 5.

Cartoons summarising the proposed conceptual model linking the Newdigate seismicity to reductions in fluid pressure in oil wells. (a) Cross-section showing large-scale processes. Production of oil (1) will reduce the pressure within the Portland reservoir near the production well. This will cause flow of oil towards the well from more distal parts of the reservoir (2). This will be accompanied by flow of groundwater into the volume vacated by the oil, which is inferred to be hydraulically connected (presumably by fractures) to the ‘hyper-permeable’ calcite ‘beef’ fabric in the underlying lower Portland sandstone a short distance below (3), which will cause flow within this fabric (4). This ‘beef’ fabric is inferred to be continuous as far as the Newdigate fault zone. The flow within it will therefore draw groundwater from greater depths (5), up one or more strands of this fault zone, which is assumed permeable, reducing the pressure in the section where this fault transects the Dinantian limestone. This pressure reduction will act to draw groundwater from the permeable Dinantian limestone into the fault (6). The associated compaction of the Dinantian limestone will cause separation of its two surfaces across the seismogenic fault (7). Surface interventions affecting the wells, such as bleeding pressure following shut-in, will reduce the pressure inside the well and have a similar overall effect. In contrast, the Kimmeridge reservoir is below the ‘beef’ layer, so no corresponding downward hydraulic connection from this reservoir exists, therefore production from this reservoir causes no equivalent effect on the hydrology and state of stress within the Dinantian limestone. (b) Plan view of a patch of the seismogenic fault showing processes on a micro-structural scale on the seismogenic fault, where separation of the fault surfaces by a small distance Δx, from configuration (i) to configuration (ii), affects three model asperities (1, 2 and 3). After this change, at asperity 1 the fault surfaces are no longer in contact, at asperity 2 what was an interference fit between the fault surfaces has become a clearance fit, and the rocks forming at asperity 3 have decompressed elastically, so they remain in contact but with a reduced normal stress and thus a reduced limiting shear stress that can maintain fault stability.

2. Geological structure and stratigraphy

The study area is in southeast England, near the boundary between the counties of Surrey and West Sussex, ∼40 km WSW of central London, on the northern flank of the Weald Basin (Figures 1 and 6). The outcrop geology and shallow subsurface structure of this area are documented by Dines and Edmunds [14] and Gallois and Worssam [15]; Trueman [16], DECC [17], and others have discussed the history of petroleum exploration. Many authors have discussed the origin and structure of the Weald Basin, or Weald sub-basin of the wider Wessex Basin (e.g., [12, 13, 18, 19, 20, 21, 22, 23, 24, 25]). As these and many other works demonstrate, this basin has developed near the northern margin of the Variscan orogenic belt, Variscan reverse faults having been reactivated as normal faults during the Mesozoic. Chadwick [19] resolved two phases of Mesozoic extension, during the Early Jurassic (Hettangian to Toarcian; extension factor, β, 1.12) and Late Jurassic and earliest Cretaceous (late Oxfordian to Valangian; β = 1.10). The Jurassic and Cretaceous sedimentary formations that accompanied and followed this extension are documented in many works and summarised in the British Geological Survey (BGS) stratigraphic lexicon (https://www.bgs.ac.uk/lexicon/); Table 1 summarises the local stratigraphy, based on the HH1 well record. This basin experienced Cenozoic inversion, when some of the Mesozoic normal faults were reactivated as reverse faults (e.g., [23]; Figure 6). As a result of this history, some faults have normal offsets within the syn-rift succession but show reverse slip in younger sediments, as illustrated in Figure 2.

Figure 6.

Map of the structure of British National Grid 100 km × 100 km quadrangle TQ, showing the depth of base Jurassic (in feet below O.D., with contours at 200 ft. or ∼ 60 m intervals) and locations where the base Jurassic is offset by faults. Faults in the vicinity of the present study area are named: C and M correspond to the Crawley and Maplehurst faults [15]; B and H appear to denote the Box Hill and Holmwood faults (cf. Figure 1), although the latter is misplaced. Modified from part of Fig. 4(a) of Butler and Pullan [12].

Most oil reservoirs in the study area are in the Upper Portland Sandstone, a shallow marine sandstone of Upper Jurassic age, sealed above by the overlying impermeable Purbeck Anhydrite, deposition of which reflects isolation of the Weald Basin interior from the sea (e.g., [21]; Table 1). The oil–water contact for the Horse Hill reservoir has been inferred as ∼580 m TVDSS [26], thus roughly at the mid-point of the Upper Portland Sandstone. The modelled extent of this ∼3 km square reservoir is illustrated in Figure 7; to the north and south it is sealed by normal faults (the Horse Hill Fault and Hookwood Fault) across which the Purbeck Anhydrite is downthrown, whereas in other directions it is sealed by the dip of the impermeable cap rock. This reservoir, hydraulically connected to the HH1 and HH2Z wells, from which production began in mid 2018 and late 2019, respectively (see supplement), is thus much larger than the hydraulic ‘radius of influence’ depicted by Hicks et al. [1] in Figure 1. The neighbouring Collendean Farm reservoir is to the north, separated by the Collendean Farm Fault, was reached by the Collendean Farm-1 (CF1) borehole (drilled in 1964; BGS ID TQ24SW1; at TQ 2480 4429; Figure 7). Oil has also been produced by the HH1 well from a deeper reservoir in the Kimmeridge Limestone (e.g., [26]) (the ‘Kimmeridge micrites’ in Table 1). However, the production from this reservoir shows no correlation with seismicity (Figure 4), indicating that this deeper reservoir is not hydraulically connected to the seismogenic Newdigate Fault.

Figure 7.

Structure contour map of the top Portland sandstone in the vicinity of the Horse Hill oilfield. Redrawn after part of Fig. 6.6 of Xodus [26]. The contours within the range of reservoir depths are at 530.4, 536.4, 542.5, 548.6, 554.7, 560.8, 566.9 m, 573.0, and 579.1 m TVDSS. Faults are labelled thus: CFF, Collendean farm fault; CFF2, northern strand of Collendean farm fault; HKF, Hookwood fault; NGF, Newdigate fault; and NWF, Nalderswood fault. The key earthquakes depicted are those that initiated the first and third ‘bursts’ of seismicity, for which analysis of subsurface pressure changes is undertaken: For the first at 11:10 on 1 April 2018 (B; TQ 21992 41976; depth 3.08 km; M_L 2.66; M_W 2.76); and for the third at 07:43 on 14 February 2019 (H; TQ 21959 41543; depth 2.05 km; M_L 2.47; M_W 2.52). The event that initiated the second ‘burst’ of seismicity occurred at 12:28 on 27 June 2018 (TQ 23230 42421; depth 2.39 km; M_L 2.52). Earthquake details are after Hicks et al. [1]. Locations of seismic lines are from UKOGL.

The smaller (∼2 km × ∼1 km) Brockham reservoir (Figure 8) has a permeability of up to ∼200 mD (e.g., [27]), its base at 582 m TVDSS [27]. As Figure 8 shows, it is sealed to the southeast by the Brockham Fault and in other directions by the dip of the Purbeck Anhydrite. However, although the Brockham Fault has ∼40 m of throw and acts as an impermeable reservoir seal, ∼1 km farther SW, where the top Portland Sandstone is deeper than the reservoir base, this fault has low throw and the permeable sandstone is juxtaposed on both sides. This aspect of the fault geometry (discussed in more detail in the supplement) is key to the proposed conceptual model (Figure 5), as it makes possible a permeable hydraulic connection between the Brockham reservoir and the seismogenic Newdigate Fault to the southeast. In contrast, Hicks et al. [1] depicted the Brockham Fault as sealing the reservoir to the SE for ∼3 km distance (Figure 1). Together with the different geometry of other faults in the area, this portrayal would preclude the possibility of any permeable hydraulic connection between this reservoir and the Newdigate Fault. The Angus [27] fault reconstruction is favoured in the present study, being assessed as more soundly based than either the Hicks et al. [1] version or the earlier interpretation by Europa [28] (see supplement). As detailed in the supplement, during prolonged production from Brockham well BRX2Y from 1998 to 2016, amounting to 36,900 m³ [11], the reservoir pressure had decreased from ∼900 to ∼500 psi. Angus [29] reported that at the end of production in 2016 the reservoir pressure was ∼490 psi or ∼ 3.4 MPa, indicating a decline by ∼2.7 MPa below the expected hydrostatic pressure of ∼6.1 MPa at its ∼622 m depth.

Figure 8.

Structure contour map of the top Portland sandstone in the vicinity of the Brockham oilfield. Redrawn after part of Fig. 4 of Angus [27]. Faults are labelled thus: BF, Brockham fault; BHF, Box Hill fault; and HWF, Holmwood fault.

The geometry of the Newdigate Fault also differs radically between the Xodus [26] and Hicks et al. [1] interpretations (Figures 1 and 7). Hicks et al. [1] depict it as a continuous structure with total E-W length ∼ 13 km, extending from the Newdigate area to the vicinity of Horley, its eastern part thus forming the southern seal to the Horse Hill reservoir. In contrast, Xodus [26] depicted it as dying out eastwards ∼3 km SW of the HH1 well, at the eastern end of the Newdigate seismicity (Figure 7). In this interpretation the fault that forms the southern seal of the Horse Hill reservoir (here designated the Hookwood Fault, HKF; Figure 7) is separated from the Newdigate Fault (NGF) by a ∼ 1 km ‘gap’ around Charlwood. To investigate this difference, the records from seismic lines C79–36 and TWLD90–21, which run N-S through this area, were obtained from the UK Onshore Geophysical Library (OGL; https://ukogl.org.uk/) archive. Neither seismic line shows any stratigraphic offset transecting this area, favouring the Xodus [26] interpretation, thus also casting doubt on other aspects of the Hicks et al. [1] analysis where they differ from other interpretations (such as for the aforementioned Brockham Fault), and initiating thorough checking of their article, which identified other problematic aspects (see Figure 1 caption, also below).

The base of the Jurassic sequence lies at ∼2100–2200 m depth in the study area (e.g., [12, 13]; Figure 6). This sequence is locally underlain by thin Triassic deposits overlying pre-Variscan (Palaeozoic) ‘basement’ at depths of > ∼ 2200 m [30]. The uppermost ‘basement’ in much of this area is known from borehole records to be Dinantian (Early Carboniferous) limestone [13, 30]. Thus, the HH1 well log (Table 1) indicates that the Jurassic Lias Group is underlain by ∼60 m of latest Triassic Penarth Group (‘Rhaetic’) rocks, then ∼50 m of the Triassic Mercia Mudstone, then ∼10 m of dolomitic conglomerate of uncertain age, then ∼70 m of Dinantian limestone, above a mudstone-dominated Upper Devonian succession. Busby and Smith [30] estimated using gravity modelling that these Devonian rocks are typically ∼1–2 km thick, their base at a typical depth of ∼3.5–4 km, being underlain in the central Weald Basin by many kilometres of Lower Palaeozoic metamorphic basement. Around the northern margin of the basin and the southern margin of the adjoining London Platform the Dinantian limestone and underlying Devonian rocks are well imaged seismically at <1 s two way time (TWT), indicating depth < ∼1500 m (e.g., [18]). Both these subdivisions are locally several hundred metres thick, the limestone being relatively unreflective and the Devonian succession highly reflective. Moving southward, as the overlying Mesozoic succession thickens, the Dinantian limestone gradually becomes thinner and its boundaries become more difficult to interpret seismically (e.g., [18]). Busby and Smith [30] noted reports of this limestone in many boreholes beneath the Weald Basin; in their view it persists southward beneath most of the basin, almost to the English Channel coastline.

More recently, Pullan and Butler [13] have presented a new map (their Fig. 21) showing the pre-Variscan subcrop beneath much of the Weald Basin (including most of the area of Figure 1) as Devonian, the Dinantian limestone being inferred to be absent. As interpreted by these workers, this limestone dies out ∼2.5 km SW of the HH1 well, indicating that it is absent in the vicinity of the Newdigate fault and the associated seismicity. However, seismic lines in this area (e.g., that in Figure 2, which is part of line TWLD90–15), as well as lines C79–36 and TWLD90–21 that have also been scrutinised as part of this study, do not clearly resolve whether this limestone is present or not; as Pullan and Butler [13] showed, there is no well control for ∼20 km distance SW of the HH1 well, so no direct evidence either way. Pullan and Butler [13] noted dipmeter evidence that in the HH1 well this limestone dips northward at ∼20–30°; their inference that it dies out not far away seems based on structural projection given its thickness (Table 1) and assuming continued northward dip. However, it is clear from other seismic sections (e.g., [18]) that in other parts of the basin the Dinantian limestone is folded. At this stage it is unclear whether this lithology is present across the study area or not. The proposed conceptual model (Figure 5) requires a highly permeable lithology, such as this, beneath the Mesozoic sediment in this area.

The two key issues already noted will now be addressed. The distribution and properties of calcite ‘beef’ will first be discussed. Second, the geolocation of features, mislocated by Hicks et al. [1], will be considered.

2.1 Calcite ‘beef’ and its significance

Calcite ‘beef’, first reported in 1826 by Webster [31], consists of bedding-parallel veins of diagenetic calcite (e.g., [32, 33]). In 1835, Buckland and De la Beche [34] adopted this nomenclature for veins of fibrous calcite within claystone beds in what is now known as the Purbeck Group in Dorset, the term ‘beef’ having originally been used by quarry workers on account of similarity to the fibrous structure of meat. This fabric (illustrated by many authors, including [35, 36, 37]), is now recognised in mudstone formations worldwide (e.g., [35]). Following the above-mentioned early reporting its mode of origin was widely debated; the view has become accepted relatively recently that ‘beef’ develops by natural hydraulic fracturing associated with overpressure during hydrocarbon maturation and migration (e.g., [32, 33, 35, 36, 38, 39, 40, 41]; cf. [36, 42]). This fabric is indeed sometimes designated as ‘hydrocarbon-expulsion fractures’ (e.g., [43]). The conditions for calcite ‘beef’ development include palaeo-temperature ∼ 70–120 °C [35]. In the central Weald Basin, such conditions are expected throughout the Jurassic succession, given the estimated ≥2 km of burial during the Cretaceous, before the Cenozoic denudation (e.g., [18]). The idea that the properties of calcite ‘beef’ enable the Newdigate Fault and neighbouring oil reservoirs to be hydraulically connected was suggested by Geosierra [44], but the present study proposes a different physical mechanism as the cause of the hydraulic connection (Figure 5).

In southern England, calcite ‘beef’ is best known in the Early Jurassic Shales-With-Beef Member (https://www.bgs.ac.uk/lexicon/lexicon.cfm?pub=SHWB) of the Charmouth Mudstone Formation, part of the Lias Group, which crops out around Lyme Regis in Dorset (e.g., [36, 45, 46, 47]). Calcite ‘beef’ is also well known in the Late Jurassic of the Weald Basin from both outcrop and borehole sections (e.g., [48]). In the Howett [48] stratigraphy, this fabric occurs within the ‘Shales with Beef and Clay-ironstone’ unit, which occurs at the top of the Middle Purbeck succession and is typically ∼20 m thick.

This fabric (reported as ‘calcite veining’) is also known from older Late Jurassic deposits, for example in core recovered between 701 and 710 m depth (below ground level 80.3 m O.D., so at 621–630 m TVDSS) in the Collendean Farm borehole near the Horse Hill site (Figure 1), in glauconitic sandstone forming the lower part of the Portland Group. Gallois and Worssam [15] placed this stratigraphic level in what they regarded as the sandy upper part of the underlying Kimmeridge Clay Formation. Nonetheless, in recent petroleum exploration reports (e.g., [26]), as in Table 1, this glauconitic sandstone with calcite ‘beef’ is reinstated within the Lower Portland Sandstone. Its inclusion within the Portland Group explains why this group is portrayed as much thicker in the recent petroleum-oriented literature (e.g., ∼130 m thick in Table 1) than by Gallois and Worssam [15], who stated its thickness as only 54 m at Collendean Farm. As these latter authors noted, the Portland Group in the Weald Basin is not well correlated with the ‘type’ Portlandian of the Portland area of Dorset, which is in the Portland – South Wight Basin (e.g., [21]). The ‘type’ Portlandian includes the Portland Limestone (now known as the Portland Stone Formation), an important building stone, the sediments of this age not being sandstone-dominated as in the Weald Basin.

The significance of all the above for the present study is as follows. It has previously been noted that the processes responsible for ‘beef’ formation will create permeability anisotropy, permeability being far greater parallel to the fabric and bedding than in the perpendicular direction (e.g., [39, 49]). Various workers have estimated the permeability of such bedding-parallel fractures, the highest estimate identified during the present work, ∼900 mD (∼9 × 10⁻¹³ m²), being by Carey et al. [50] for the Ordovician Utica Shale of eastern North America. This is many orders-of-magnitude higher than the expected nanodarcy permeability of shale perpendicular to bedding, and is quite a high value for rocks in general.

2.2 Geolocation

The study area has been illustrated using the map (Figure 1), and seismic cross-section (Figure 2) from Hicks et al. [1]. However, the original versions of both these figures have required significant amendment regarding accuracy issues. This map was originally geolocated using geographical co-ordinates; to make it easier to use British National Grid (BNG) co-ordinates have been added. This map also shows seismic lines and faults. The information source for seismic lines, including line TWLD90–15 that is illustrated in Figures 2 and 3, was not reported by Hicks et al. [1]; it is evident that they are from the UKOGL location map. Hicks et al. [1] also explained that (in lieu of using the existing literature) they identified faults in the study area through their own interpretation of seismic lines. As already noted, regarding key aspects of the structure the existing interpretations are favoured over these revisions by Hicks et al. [1].

The seismic section in Figure 2 clearly has higher resolution than older ones, including those which informed earlier fault maps such as that by Butler and Pullan [12] (Figure 6). Some of the faults depicted in Figure 2 are, thus, recognised for the first time. However, additional faults are also evident in the uninterpreted version provided by UKOGL (Figure 3). It is thus evident that in the lower part of the Jurassic sediment and upper part of the underlying Palaeozoic basement, the Newdigate Fault consists of multiple fault strands distributed across a zone with N-S width approaching ∼2 km. Careful inspection of supplementary Figure S13 of Hicks et al. [1] and Figure 3 indicates broken and offset seismic reflectors which delineate these subsidiary strands of the Newdigate fault zone, some evidently near the limit of seismic resolution (cf. [51]), which merge upwards by ∼0.5 s two way time (TWT).

A significant issue to have emerged from checking the Hicks et al. [1] analysis concerns their velocity model used for earthquake location (Table 2). As Figure 9(a) shows, this is significantly slower than is expected from the sonic logs for wells BR1, CF1 and HH1, and from the recent depth-conversion analysis by Pullan and Butler [13]. The Hicks et al. [1] velocity model is also significantly slower that that obtained from moveout analysis during the processing of seismic line TWLD90–21 (Figure 9(b)). Hicks et al. [1] explained that they based their velocity model on existing regional models, but ‘then improved [it] using sonic logs from the Brockham and Horse Hill wells’. However, as Figure 9(a) shows, despite this workflow, their velocity model is inconsistent with the records from these wells; they evidently somehow made a mistake over this aspect. Being too slow, their velocity model causes hypocentres to be mislocated too shallow. Figure 2 includes a depth scale using the preferred velocity model in Table 1, based on the HH1 log. At the position of the earthquakes, this depth scale is ∼400 m deeper than that provided by Hicks et al. [1]. Furthermore, as is discussed in the supplement, Hicks et al. [1] appear to have depth-converted this seismic section using the same velocity model (Table 2) that they used for earthquake location. As is also discussed in the supplement, location of the earthquakes using the velocity model in Table 1 would adjust their hypocentres downward by an estimated ∼400 m. The hypocentres are thus positioned more-or-less correctly relative to the detail in the seismic section, but they and this detail are now placed ∼400 m deeper than Hicks et al. [1] envisaged. This adjustment moves the earthquake population from within the Jurassic sedimentary section to within the Palaeozoic ‘basement’. These earthquakes presumably occurred on one or more of the steeply north-dipping subsidiary strands of the Newdigate fault zone, given the steeply north-dipping nodal planes, identified as the fault planes, of the focal mechanisms (Figure 1 and Table 3), rather than on the main Newdigate Fault that dips south.

Subdivision	MD (m)	TVDSS (m)	TWT (s)	VI (m s^−1)	Notes
Younger subdivisions (Early Cretaceous; Berriasian to Barremian)
Weald Clay	7.6	-66.9	NR	ND
Hastings Beds	158.5	84.0	NR	ND
Grinstead Clay	211.8	137.3	NR	ND
Lower Tunbridge Wells Sands	234.7	160.2	NR	ND
Wadhurst Clay	245.4	170.8	NR	ND
Ashdown Beds	298.1	223.6	NR	ND
Purbeck Group (latest Jurassic and earliest Cretaceous; Tithonian and Berriasian)
Purbeck Durlston Beds	396.8	322.3	NR	ND	1
Purbeck Carbonates	464.8	390.3	NR	ND
Purbeck Main Anhydrite	604.7	530.2	0.370	2500
Portland Group (Late Jurassic; Tithonian)
Upper Portland Sandstone	622.4	547.7	0.384	2531
Lower Portland Sandstone	708.4	632.5	0.451	5011
Ancholme Group (latest Middle Jurassic and Late Jurassic; Callovian to Tithonian)
Kimmeridge Clay	755.9	677.6	0.469	2787
Kimmeridgian Micrite 1	851.3	765.4	0.532	2861
Kimmeridgian Micrite 2	939.7	835.5	0.581	2961
Top Corallian	1359.1	1139.0	0.786	3289
Corallian Limestone	1523.7	1272.2	0.867	3743
Oxford Clay	1539.2	1285.3	0.874	3540
Kellaways Beds	1666.0	1403.9	0.941	3725
Great Oolite Group (Middle Jurassic; Bathonian and Callovian)
Cornbrash	1681.6	1418.8	0.949	2600	2
Main Great Oolite	1682.8	1420.1	0.950	5095
Fuller’s Earth	1732.8	1468.5	0.969	4886
Inferior Oolite Group (Middle Jurassic; Aalenian and Bajocian)
Inferior Oolite	1767.5	1502.7	0.983	5584
Lias Group (Early Jurassic; Hettangian to Toarcian)
Upper Lias	1941.6	1675.8	1.045	4244
Middle Lias	2048.0	1781.9	1.095	4796
Lower Lias	2158.3	1892.2	1.141	4301
Older subdivisions (Triassic and older)
Rhaetic	2470.1	2204.0	1.286	5318
Mercia Mudstone	2528.6	2262.5	1.308	4434	3
Dolomitic Conglomerate	2581.7	2315.6	NR	ND
Carboniferous Limestone	2593.2	2326.8	1.337	ND
Upper Devonian	2659.4	2393.3	NR	ND	4
TD	2686.8	2420.7	NR	ND	5

Table 1.

Stratigraphy of the Horse Hill 1 borehole.

Data for tops of stratigraphic subdivisions (as used by UKOGL; not all expressed using modern formal stratigraphic nomenclature, which is available from https://www.bgs.ac.uk/lexicon) are from the online well log (https://ukogl.org.uk/map/php/pdf.php?subfolder=wells\tops&filename=3041.pdf), supplemented by values from Pullan and Butler [13], NR indicating ‘not reported’. Measured Depth (MD) is measured below a datum at 66.9 m O.D., below the local ground level of 74.5 m O.D. at the wellhead, at TQ 25254 43600. TVDSS is True Vertical Depth below O.D.; TWT is echo time. Values of interval velocity, V_I, are determined in this study, ND indicating ‘not determined’. Notes: 1. The Durlston Beds (or Durlston Formation) are nowadays regarded as earliest Cretaceous; the rest of the Purbeck Group is Late Jurassic. 2. The Cornbrash Formation is too thin here for its interval velocity to be reliably determined. 3. Interval velocities for the Mercia Mudstone and the Dolomitic Conglomerate combined. 4. TVDSS for the top Devonian estimated given the vertical orientation of the deepest part of the well. 5. The well bottoms (at TD) in Upper Devonian mudstone.

H (km)	VP (km s^−1)	VS (km s^−1)
0.0	2.2	1.2
0.2	2.4	1.4
0.4	2.6	1.5
0.7	2.7	1.5
1.2	3.1	1.8
1.5	3.6	2.0
1.8	4.7	2.7
2.1	5.0	2.8
2.4	5.5	3.1
7.6	6.4	3.7
18.9	7.0	4.1
34.2	8.0	4.6

Table 2.

Velocity Model from Hicks et al. [1].

This velocity model was used by Hicks et al. [1] for earthquake relocation and moment tensor inversion. H denotes the depth to the top of each layer; V_P and V_S denote the P-wave and S-wave velocities. Note that this velocity model is significantly slower that that in Table 1; it results in a two-way time to depth 2326.8 m, corresponding to the top of the Carboniferous Limestone at Horse Hill, of 1.466 s rather than the actual 1.337 s.

Figure 9.

Comparisons of one-way time versus depth for different seismic velocity models. (a) Comparison of the Hicks et al. [1] velocity model with the sonic logs from wells CF1, HH1, and BR1 (from [24], and UKOGL), the velocity model (based on well HH1) adopted for this study, and a representative point from the analysis by Pullan and Butler [13] where a ∼ 600 ms one-way time corresponds to a depth of ∼2100 m. (b) Comparison of the Hicks et al. [1] velocity model with velocity models derived from interval velocities obtained from the moveout analysys of seismic line TWLD90–21 (Fig. 7), for a suite of representative CDPs. These interval velocities were calculated in this study using Dix’s formula (e.g., [52, 53]) from root mean square velocities provided by UKOGL. UKOGL did not have such data for the other seismic lines analysed for this study.

No.	Date	Time; (UTC)	Epicentre; (BNG reference)	Depth; (m)	ΔN; (m)	ΔE; (m)	Δz; (m)	Δα; (°)	NP	NS	Δto; (s)	zDD; (m)	ML	MW; 1	MW; 2	zC; (m)	Strike; (°)	Dip; (°)	Rake; (°)	fC; (Hz)	ro; (m)	Δσ; (MPa)
1	18 Jul 2018	03:59:56	TQ 22005 41393	1990	1397	803	1011	151	15	13	0.04	ND	2.01	2.20	2.03	2.00	282	74	178	6.4	48	0.25
2	18 Jul 2018	13:33:18	TQ 21920 41474	1860	1463	737	1014	145	15	15	0.02	ND	2.54	2.56	2.45	2.20	276	75	169	4.3	25	0.27
3	14 Feb 2019	07:43:33	TQ 22959 41543	2220	297	330	379	98	9	7	0.09	2050	2.47	2.52	2.27	2.80	255	86	173	7.7	142	1.05
4	19 Feb 2019	17:03:57	TQ 22872 41538	2050	220	429	393	106	5	4	0.04	2040	1.98	1.95	1.77	2.20	256	61	-163	11.2	86	0.56
5	27 Feb 2019	03:42:21	TQ 22622 41517	2110	286	352	316	98	14	11	0.13	2300	3.18	3.25	2.87	3.60	260	78	178	5.8	169	2.62
6	04 May 2019	00:19:19	TQ 22796 41516	2190	143	165	294	94	13	10	0.11	2440	2.35	2.31	2.17	2.40	255	85	167	16.9	64	8.06

Table 3.

Source parameters for Newdigate earthquakes with focal mechanisms.

Data from supplementary Table S2 of Hicks et al. [1]. The events are numbered to match Figure 1. For events 1 and 2 only conventional hypocentral locations were determined, which yielded the epicentral co-ordinates and focal depths. Double difference focal depths (z_DD) were not determined (ND). For the other events, the epicentral co-ordinates and z_DD are determined from the double difference location procedure and the ‘Depth’ by conventional location. ΔN, ΔE and Δz are the uncertainties in northing, easting, and depth, based for events 3–6 on the double difference solutions. Δα is the maximum gap between ray path azimuths to seismograph stations that recorded each event, N_P and N_S being the numbers of P- and S-wave records. Δt_o is the rms residual in origin time. M_L and M_W are local magnitude and moment magnitude. M_W 1 and the centroid depth z_C are determined from the moment tensor; M_W 2 is from P-wave spectra. Strike, Dip and Rake are for the nodal plane of the focal mechanism that is regarded as the fault plane, being subparallel to the Newdigate Fault (Figure 2). Mean corner frequency f_C, source radius r_O, and stress drop Δσ are determined from seismogram spectra.

3. Seismicity and its correlation with well activities

As already noted, multiple publications have already documented the 2018–2019 Newdigate ‘earthquake swarm’, notably those by Baptie and Luckett [54], Verdon et al. [6], and Hicks et al. [1]. Baptie and Luckett [54] presented a preliminary analysis of 14 earthquakes between 1 April and 18 August 2018; their results informed the OGA workshop. The more extensive analysis by Hicks et al. [1] will now be appraised. These latter authors determined hypocentres and other source parameters for 168 earthquakes between 1 April 2018 and 28 June 2019, some with local magnitude M_L < −1, their location patterns and timeline being depicted in Figures 1,2 and 4 and summarised in Table 3, with Table 4 listing events that post-date their study. The first nine earthquakes up to 10 July 2018 (including one of the largest, with M_L 3.02, on 5 July) were located before any local seismograph stations were operational, using only data from regional stations. Hicks et al. [1] explained that due to the limited available data these events were located by assigning each a fixed focal depth. The resulting reported depths vary between 2.33 and 3.08 km (see Table S2 of [1]), it being unclear on what basis different depths were assumed for different events. The next sixteen events, until 11 July, were located conventionally but including data from local stations. For the rest of the events, both ‘double difference’ relocations (after [55]) and conventional locations were determined, using the velocity structure in Table 2. As detailed in Figures 1 and 2, most of the earthquakes in a zone ∼1.5 km long (E-W) by ∼300 m wide (N-S). The compact width of this zone suggests that many patches of a single strand of the Newdigate fault zone were reactivated.

Date	Time (UTC)	Latitude; (°N)	Longitude; (°W)	BNG	Depth; (km)	ML	Note
9 June 2019	02:43:18.2	51.159	0.237	TQ 23382 41449	2.7	-0.5
9 June 2019	23:00:15.0	51.133	0.295	TQ 19393 38462	3.3	-0.1
6 July 2019	01:03:20.4	51.161	0.242	TQ 23027 41663	2.5	-0.7
6 July 2019	01:03:23.7	51.161	0.242	TQ 23027 41663	2.5	-0.7
6 July 2019	01:03:30.1	51.159	0.241	TQ 23102 41442	2.1	-0.8
6 July 2019	01:03:40.2	51.159	0.241	TQ 23102 41442	2.2	-0.7
6 July 2019	03:57:15.3	51.160	0.239	TQ 23239 41557	2.5	0.1
20 July 2019	22:02:26.0	51.158	0.251	TQ 22405 41315	2.1	-0.6
29 July 2019	03:35:25.5	51.160	0.242	TQ 23029 41552	2.2	-0.1
6 Aug 2019	02:32:00.9	51.157	0.239	TQ 23247 41223	2.2	-0.5
12 Aug 2019	00:46:46.6	51.160	0.241	TQ 23099 41554	2.1	-0.7
12 Aug 2019	00:46:49.2	51.160	0.241	TQ 23099 41554	2.1	-1.4
2 Sep 2019	05:13:04.9	51.160	0.237	TQ 23379 41560	2.0	1.1	1
3 Sep 2019	20:19:13.2	51.161	0.237	TQ 23376 41672	2.0	0.2
6 Sep 2019	07:09:45.5	51.161	0.237	TQ 23376 41672	2.0	1.0
21 Sep 2019	14:43:45.2	51.160	0.237	TQ 23379 41560	2.2	0.6
31 Oct 2019	19:25:16.4	51.160	0.238	TQ 23309 41558	2.0	0.8
29 Apr 2020	00:11:25.6	51.172	0.256	TQ 22019 42863	3.1	0.3
24 May 2020	15:16:56.9	51.157	0.250	TQ 22478 41205	2.4	0.6

Table 4.

Newdigate seismicity since the start of June 2019.

Cataloguing here is complete to 27 August 2020. Data are from https://earthquakes.bgs.ac.uk; these earthquakes have been located using standard BGS procedures, as reported by the International Seismological Centre (http://www.isc.ac.uk). Co-ordinate transformations to the British National Grid, as part of this study, use https://www.bgs.ac.uk/data/webservices/convertForm.cfm

Note: 1. Felt in Newdigate; maximum EMS intensity 2.

Focal mechanisms were determined by Hicks et al. [1] for six events, including the largest, of M_L 3.18 and moment magnitude M_W 3.25, on 27 February 2019, as illustrated in Figure 1 and listed in Table 3. All six events have a nodal plane striking roughly east–west and dipping steeply north. As already noted, this plane is inferred to be the fault plane, indicating predominant right-lateral slip. Available data regarding the state of stress in the Weald Basin are extremely limited; Kingdon et al. [56] and Fellgett et al. [57] provided syntheses of in situ stress data across much of Britain. However, these authors wrote little about the Weald Basin, Fellgett et al. [57] noting that many hydrocarbon wells in this area have yielded stress data but these data had not yet been placed in the public domain. The stress dataset available for the Weald Basin thus remains that presented by Evans and Brereton [58]. As is detailed in the supplement, this limited dataset indicates a NW-SE maximum principal stress and a NE–SW minimum principal stress. The Newdigate focal mechanisms (Figure 1) are consistent with this stress field orientation, given the standard requirement for the maximum principal stress to lie within dilatational quadrants [59].

3.1 Temporal clustering

As detailed by Hicks et al. [1], the Newdigate seismicity between April 2018 and June 2019 involved four ‘bursts’ of activity (Figure 4). The first began at 11:10 on 1 April (M_L 2.66), followed by two events later on the same day, another on 9 April, and a final event on 28 April. The smallest of these events (on 9 April) had M_L 1.28. No local seismograph stations were then in operation; Hicks et al. [1] estimated that the completeness threshold for earthquake detection was circa M_L 2, so many smaller events were undoubtedly missed.

The second ‘burst’ (Figure 4) began at 12:28 on 27 June (M_L 2.52), and included four other events above M_L 2.0 (on 29 June and 5 July, and two on 18 July) including the second largest event overall (M_L 3.02), at 10:53 on 5 July. The installation of local seismograph stations in mid June and early July lowered the completeness threshold for earthquake detection to below M_L 0 [1], resulting in many small events being thereafter detected and enabling use of the aforementioned relative location procedure. After these initial relatively large events this ‘burst’ of earthquakes began to tail off, in terms of both magnitude and frequency of occurrence. The last event with M_L > 0 occurred at 03:21 on 18 August (M_L 0.30), with infrequent smaller events persisting into early 2019.

The third ‘burst’ (Figure 4) began on 14 February 2019 with a relatively large event at 07:43 (M_L 2.47), followed by two other events of M_L ≥ ∼2, at 17:03 on 19 February (M_L 1.98) and at 03:42 on 27 February (M_L 3.18), this being the largest event of the overall sequence. After these initial relatively large events this ‘burst’ also tailed off, although two events with M_L > 0 occurred during April 2019 (on 11 and 22 April; M_L 0.73 and 0.56).

The fourth ‘burst’ (Figure 4) began with a relatively large event (M_L 2.35) at 00:19 on 4 May 2019. As for the preceding ‘bursts’, this seismicity thereafter began to tail off, although events with M_L ∼ 0 persisted until the end of June 2019. Locations by BGS confirm the tailing-off trend through July and August 2019 (Table 4), with a M_L 1.1 event on 2 September, three smaller events later that month, one during October, and none more before the end of 2019 (Table 3).

Overall, this pattern of seismicity, consisting of ‘bursts’ of events, each involving activity tailing off after a peak, with the largest event increasing during successive ‘bursts’, bears a striking resemblance to other earthquake swarms that are inferred to be caused by fluid pressure changes in a fault (e.g., [60]). However, the Newdigate earthquake population is insufficient to permit statistical testing of the patterns expected for this mechanism.

3.2 Correlation of seismicity with well activities

Figure 4(c) indicates how these four ‘bursts’ of earthquakes correlate with activities affecting the Portland reservoir in the HH1 or BRX2Y wells. Production from well BRX2Y resumed in late March 2018: from Hicks et al. [1] ∼4.0 m³ (∼25 barrels) of oil were produced on 23 March followed by ∼1.1, ∼0.9 and ∼1.0 m³ (∼7, ∼6 and ∼6 barrels) on 25–27 March. Reservoir pressure during this and subsequent production has not been reported, but from standard theory (e.g., [61, 62]) one expects it to have again decreased. This start of production occurred nine days before the first Newdigate earthquake on 1 April 2018. Furthermore, as is detailed in Figure 4 and in the online supplement, other brief ‘pulses’ of production occurred from well BRX2Y in June, respectively 20, 19, 16 and 6 days before the start of the second ‘burst’ of seismicity on 27 June.

Although the activities that were planned in the HH1 well in 2018 have been disclosed [63], most of the actual activities that took place, and any associated pressure variations in the well, have not been, other than in the very general terms reported by Hicks et al. [1]. An attempt is made in the supplement to piece together the sequence of events, based on fragments of information available. It is thus evident that before 4 July 2018, the Portland reservoir was reported as isolated from the surface by a removable bridge plug. Claims have been made that the reservoir might have been influenced before this date by surface activities and by activities in the shallow part of the well [2]; if so, this would imply that the bridge plug had failed.

It is evident from Figure 4 that production ceased from well BRX2Y in October 2018; production at HH1 switched from the Portland reservoir to the Kimmeridgian reservoirs around the same time. Around this time, seismicity at Newdigate tailed off significantly. The conceptual model in Figure 5 provides a natural explanation for such a variation.

The seismicity then re-intensified as the third ‘burst’, recognised by Hicks et al. [1], starting on 14 February 2019, which followed the resumption on 11 February 2019 of production from well HH1, now at rates of up to 220 bopd, from the Portland reservoir. As is detailed in the online supplement, production from this reservoir continued until late June 2019, after which it switched back to the Kimmeridgian reservoir, then during December 2019 to the newly-completed horizontal lateral, off well Horse Hill-2 (designated HH2Z), in the Portland reservoir. Seismicity at Newdigate remained significant during this phase of production from the Portland reservoir at HH1. However, production was not continuous; Hicks et al. [1] reported shutdowns during 9–12 April and 4–10 May, the latter corresponding to the start of the fourth ‘burst’ of seismicity as recognised by these authors. Seismicity subsequently tailed off following the end of production at HH1 from the Portland reservoir in late June 2019 and the switch to production from the Kimmeridgian reservoir in early July (Figure 4 and Table 4). Furthermore, seismicity did not resume during the initial flow testing of well HH2Z in December 2019, even though the production rates from the Portland reservoir were much greater, up to 1087 barrels of fluid per day, than they had been from well HH1 (see supplement).

Overall, the correlation between phases of production from the Portland reservoir, from well HH1 or well BRX2Y or both, and ‘bursts’ of seismicity has been compelling (Figure 4). Hicks et al. [1] did not recognise this pattern, apparently because they did not differentiate between the Portland and Kimmeridgian reservoirs as sources of production from well HH1, as is now done (based on details in the supplement). There are particularly clear patterns for the first and third ‘bursts’: the first began 9 days after the resumption of production from well BRX2Y in March 2018; the third began 3 days after the resumption of production from well HH1 in February 2019. However, the patterns are less clear for the other two ‘bursts’ of seismicity, nor has the flow testing of well HH2Z, starting in December 2019, been associated with any significant seismicity.

4. Conceptual geomechanical model

The conceptual geomechanical model already summarised (Figure 5), which might account for seismicity beneath Newdigate, caused by pressure decreases in the Portland reservoir resulting from production (or other activities) from the HH1 or BRX2Y wells, will now be developed quantitatively. The basis of this model (Figure 5) is as follows. The Upper Portland Sandstone reservoir adjoining these wells is assumed to make a subhorizontal hydraulic connection with the seismogenic Newdigate fault zone via a permeable fabric formed in calcite ‘beef’ in the stratigraphically adjacent Lower Portland Sandstone. The seismogenic fault is assumed highly permeable and to provide a downward hydraulic connection to the rocks beneath the Jurassic succession. These rocks are assumed to include the dolomitic conglomerate and Dinantian limestone, as in the HH1 well (Table 1), which are themselves permeable. It is further assumed that the Newdigate seismicity has occurred at locations where these permeable lithologies are in contact across this fault. Pressure reduction in the Portland reservoir will thus be communicated via the ‘beef’ fabric, reducing the fluid pressure in this fault, which will cause flow from within the adjoining permeable lithologies into the fault. The associated reduction in fluid volume within these lithologies will cause them to compact. This will result in surfaces that were previously in contact across this fault to separate slightly, reducing the normal stress across the fault. This will ‘unclamp’ the fault (as in Figure 5(b)), moving it closer to the Mohr-Coulomb failure condition. The fault is itself assumed to be ‘critically stressed’, already near this failure condition, potentially enabling relatively small changes in the state of stress to cause coseismic slip (cf. [64, 65]).

Regarding the assumptions thus made, the presence of calcite ‘beef’ within the Portland Group sediments and its permeability have already been discussed. The permeability of faults is a major issue in Earth science (e.g., [66, 67, 68, 69, 70]). There is no information regarding the permeability of any strand of the Newdigate fault zone; however, although counterexamples exist (e.g., faults made impermeable by cemented fault gouge [71]), the view that faults are generally permeable, especially when critically stressed (e.g., [72]) is widely accepted, as is the precautionary principle that faults are assumed permeable, in the absence of contrary evidence, when assessing the possibility of subsurface fluid migration (e.g., [73, 74]). In the present study area, faults with offsets of tens of metres, where the permeable Portland Sandstone is juxtaposed against the impermeable Purbeck Anhydrite, act as seals for oil reservoirs (e.g., [26, 27]) and are obviously impermeable. However, low-offset faults with the permeable Portland Sandstone juxtaposed on both sides can be expected to be permeable. The question of the continuity of the Dinantian limestone from the HH1 area to the vicinity of the seismogenic strand of the Newdigate fault zone has already been discussed. The uncertainty regarding the state of stress in the Weald Basin is considered in the online supplement. As will become clear below, if the differential stress here is anything like as high as it is the Preese Hall area (after [11]), then any fault with the orientation of that which slipped will be very close to the Mohr-Coulomb failure condition.

The model thus includes three elements: pressure drawdown caused by production in wells, and its communication via the ‘beef’ fabric; compaction of the permeable rocks alongside the Newdigate fault caused by fluid withdrawal accompanying depressurization; and the associated Mohr-Coulomb failure analysis.

4.1 Pressure drawdown accompanying production

In general variations in fluid pressure, P, in a porous medium, are solutions to the diffusion equation

$$\frac{\mathrm{\partial P}}{\mathrm{\partial t}}=D{\nabla }^{2}P,$$E1

where t is time, ∇² is the Laplacian operator, and D is the hydraulic diffusivity of the medium (e.g., [75]). The value of D depends on properties of the medium and fluid and on solution details, such as boundary conditions for pressure and strain (e.g., [76]). For pressure diffusion,

$$D\equiv \frac{kM}{\eta }$$E2

(e.g., [77]), where k and η are the permeability of the medium and viscosity of the fluid and M is the Biot modulus of the fluid-rock combination. M can be expressed as

$$\frac{1}{M}\equiv \frac{1}{B_R}-\frac{B_E}{{B_R}^2}+\varphi \left(\frac{1}{B_F}-\frac{1}{B_R}\right),$$E3

where B_R and ϕ are the bulk modulus and porosity of the rock, B_F is the bulk modulus of the fluid, and B_E is the ‘effective’ bulk modulus of the combination, defined as

$$\frac{1}{B_E}\equiv \frac{1-\varphi }{B_R}+\frac{\varphi }{B_W}.$$E4

Costain [78] expressed D as

$$D\equiv \frac{k{B}_E}{\eta },$$E5

which can be compared with Eq. (2); in practice, the difference in choice of elastic modulus (B_E or M) makes little practical difference (see below); from Eqs. (3) and (4),

$$\frac{1}{M}\equiv \frac{1}{B_E}-\frac{B_E}{{B_R}^2}.$$E6

The pressure variations in the calcite ‘beef’ layer hydraulically connected to, and surrounding, the Brockham and Horse Hill Portland reservoirs can be assumed to have circular symmetry; analysis in terms of cylindrical polar co-ordinates is thus required. Thus, if the flow does not vary azimuthally or across the vertical extent h of the ‘beef’, Eq. (1) reduces to

$$\frac{\partial P}{\partial t}=D\left(\frac{{\partial }^{2}P}{\partial r^2}+\frac{1}{r}\frac{\partial P}{\partial r}\right).$$E7

As others (e.g., [79, 80]) have noted, if production at rate Q starts at time t = 0, the variation in P, ΔP, after t = 0 takes the form

$$\mathrm{\Delta }P=\frac{Q\eta }{4\pi k_B{h}_B}{E}_1\left(\frac{r^2}{4D_{B}t}\right),$$E8

where the subscripts B denote values of D, h and k appropriate for the layer of ‘beef’. Here, E₁ is the Exponential Integral Function, defined as

$$E_1\left(x\right)\equiv \underset{x}{\overset{\infty }{\int }}\frac{exp\left(-s\right)}{s}\mathit{ds}.$$E9

E₁ is not supported directly in Microsoft Excel, but using its relationship to other functions (discussed, e.g., by [81]) it can be evaluated indirectly. This is possible because

$$E_1\left(x\right)\equiv {lim}_{\psi \to 0}\mathrm{\Gamma }\left(\psi ,x\right)$$E10

where Γ(ψ, x) is the Upper Incomplete Gamma Function. As Schurman [82] has noted, this means that an Excel formula providing a good approximation to E₁(x) can be written as

$$\mathit{EXP}\left(\mathit{\text{GAMMALN}}\left(B1\right)\right)\ast \left(1-\mathit{\text{GAMMA}}.\mathit{\text{DIST}}\right.\left(A1,B1,1,\mathit{\text{TRUE}}\right)$$E11

where cell A1 contains x and cell B1 contains a small positive number (say, 10⁻⁸), representing ψ. Values of E₁ thus calculated were checked against the tables by Harris [83] and Stegun and Zucker [84] and against the power series approximations for E₁(x) for the limit of x≪1,

$$E_1\left(x\right)\approx -\gamma -ln\left(x\right)+x-\frac{x^2}{4}+\cdots$$E12

(e.g., [79]), where γ = 0.5772156649… is Euler’s Constant, and for the limit x≫1,

$$E_1\left(x\right)\approx \frac{exp\left(-x\right)}{x}\left(1-\frac{1}{x}+\cdots \right)$$E13

(e.g., [85]), and were found to be accurate to four significant figures or better. This method for evaluating E₁(x) for all values of x is, thus, more accurate in general than the overall approximation formula proposed by Barry et al. [85].

In the near-wellbore volume the pressure variation will be in the Portland sandstone, not in the ‘beef’. Thus, from Eq. (8), at time t the variation at the well rim, of radius r_W, is ΔP_W where

$$\mathrm{\Delta }{P}_W=\frac{Q\eta }{4\pi k_P{h}_P}{E}_1\left(\frac{{r_W}^2}{4D_{P}t}\right),$$E14

where the subscripts P denote values of D, h and k appropriate for the Portland sandstone. Given the high value of D_P for the Portland sandstone, for most of the durations of the production pulses at BRX2Y and HH1, r_W²/(4 D_P t)≪1. One may thus approximate E₁(x) using Eq. (12). Indeed, x will be so small that only the -ln(x) term need be considered. The resulting logarithmic dependence of ΔP_W on t means that ΔP_W remains roughly constant, as observed during the HH1 well test (see supplement).

Using the same general approach, a brief episode of production at rate Q starting at time t = 0 and ending at time Δt causes a pressure variation given by

$$\mathrm{\Delta }P=\frac{Q\eta }{4\pi k_{B}h}\left(E_1\left(\frac{r^2}{4D_{B}t}\right)-E_1\left(\frac{r^2}{4D_B\left(t-\mathrm{\Delta }t\right)}\right)\right).$$E15

Using the definition of E₁ in Eq. (9), for t≫Δt, Eq. (8) can be approximated as

$$\mathrm{\Delta }P\approx \frac{Q\eta \mathrm{\Delta }t}{4\pi k_{B}ht}exp\left(\frac{-r^2}{4D_{B}t}\right).$$E16

This equation can be differentiated with respect to t; by setting the resulting partial derivative to zero one may solve for the time delay t_D for the maximum pressure variation, thus:

$$t_D=\frac{r^2}{4D_B}.$$E17

The maximum pressure variation ΔP_M at distance r and time t_D is thus

$$\mathrm{\Delta }{P}_M\left(r\right)=\frac{Q\eta D_B\mathrm{\Delta }t}{\pi k_{B}he{r}^2},$$E18

indicating that ΔP_M varies inversely with r².

4.2 Compaction alongside the Newdigate fault

The Newdigate fault is envisaged as extremely permeable, such that a pressure variation ΔP applied to any point of it by via the layer of ‘beef’ is transmitted downward, with no significant time delay, to depths where it transects the permeable Dinantian limestone, the presumed seismogenic layer. This model fault is vertical, the permeable seismogenic layer being assigned thickness H_D, hydraulic diffusivity D_D, permeability k_D, and porosity ϕ_D. Depressurization at the point where the ‘beef’ layer intersects this fault is thus inferred to cause a reduction in groundwater pressure ΔP_O at each point below on the fault within the permeable seismogenic layer. The resulting groundwater pressure variation in the Dinantian limestone will be governed by Eq. (1). However, if this variation is independent of vertical position and position parallel to the fault, the one-dimensional variant

$$\frac{\partial P}{\partial t}=D\frac{{\partial }^{2}P}{\partial x^2},$$E19

will require solution, where t is time and x is distance from the fault.

As a solution to Eq. (19), the drawdown δP in groundwater pressure within such a permeable layer is given (e.g., [78]) by

$$\mathit{\delta P}=\mathrm{\Delta }{P}_O\mathit{erfc}\left(\frac{x}{2\sqrt{\left(Dt\right)}}\right)$$E20

erfc() denoting the Complementary Gaussian Error Function. For z > 0, erfc(z) decreases as z increases, reaching ∼0.0047 when z = 2. As Detournay and Cheng [77] noted, this condition indicates an effective outer limit to significant pressure variations, at distance x_M from the model fault, where

$$x_M=4\sqrt{\left(Dt\right)}.$$E21

As time progresses, as a result of continuing production from a well that is hydraulically connected to the model fault, an ever-widening volume of rock, in the x direction perpendicular to the model fault, will thus become depressurized, water previously stored within this volume being released into the fault. Figure 10 illustrates this effect for D = 1 m² s⁻¹.

Figure 10.

Graphs of the predicted variation in pressure δP / ΔP_O with distance x from the model seismogenic fault, calculated using Eq. (20) for hydraulic diffusivity D = 1 m² s⁻¹, representing Dinantian limestone. (a) after times t of 1 hour, 12 hours, 3 days and 2 weeks. (b) after times t of 2 months, 1 year, 5 years, and 20 years.

In general, poroelastic strain responses to changes in fluid pressure can be highly complex (e.g., [77, 86, 87, 88, 89]). Segall [87] noted that in the limit where Δσ_kk = 0, a reduction in fluid pressure by δP will cause a contractional strain ε = α δP / B_E where α is Biot’s coefficient, defined as

$$\alpha \equiv 1-\frac{B_E}{B_R}.$$E22

Contraction will occur in the vertical z direction as well as in the x direction. Partitioning the contractional strain as Δε_xx = γε and Δε_zz = (1 - γ)ε,

$$\mathrm{\Delta }{\epsilon }_{\mathit{xx}}=\frac{\gamma \alpha \mathit{\delta P}}{B_E}.$$E23

Many studies of poroelastic deformation have treated petroleum or geothermal reservoirs as inclusions embedded in surrounding rocks and have treated faults as idealised planes within the inclusion or its surroundings (e.g., [86, 88, 90, 91]). Because the present study aims to explore the effect of fault asperities, the fault cannot be treated in this idealised manner. The fault is instead envisaged as a vertical ‘cut’ in the poroelastic medium, which has (distant) boundaries in both directions perpendicular to the fault plane. In this configuration, with the outer ends of the blocks on either side of the fault fixed or ‘pinned’, depressurization of pore fluid will cause their inner ends, facing each other across the fault, to separate slightly, by distance Δx, as depicted schematically in Figure 5(b). In contrast, if the rocks on either side of a vertical fault at the mid-point of a continuum model (such as that by [88]) were depressurised, these rocks would move towards the fault, the opposite sense of motion to what is required to provide the ability to make an assessment of the effect of asperities on the fault. It is for this essential reason that new theory is derived here for the pressure and stress response in the Dinantian limestone, rather than using existing published theory.

Subject to the above model definition, Δx can be estimated as

$$\mathrm{\Delta }x=2{\int }_{x=0}^{x\to \infty }\mathrm{\Delta }{\epsilon }_{\mathit{xx}}\mathit{dx},$$E24

the factor of 2 taking into account that the rocks on both sides of the fault will move away from it. Evaluation of Δx requires the integral of erfc() (cf. Eq. (20)). From Abramowicz and Stegun [92], p. 299 and Weisstein [93],

$$E\left(x\right)\equiv {\int }_{z=0}^{z=x}\mathit{erfc}\left(z\right)\mathit{dz}=x\mathit{erfc}\left(x\right)+\frac{\left(1-exp\left(-x^2\right)\right)}{\sqrt{\pi }}$$E25

so

$$E\left(\infty \right)\equiv {\int }_{z=0}^{z\to \infty }\mathit{erfc}\left(z\right)\mathit{dz}=\frac{1}{\sqrt{\pi }}$$E26

Using Eqs. (20), (23) and (26), one obtains.

$$\mathrm{\Delta }x=\frac{4\gamma \alpha \mathrm{\Delta }{P}_O}{B_E}\sqrt{\left(\frac{\mathit{DT}}{\pi }\right)},$$E27

this quantity being positive for a reduction in pressure by ΔP_O.

Costain [78] also showed that if, rather than persisting indefinitely, the pressure change ΔP_O imposed at x = 0 persists for durationδt, the resulting pressure variation δP is given by

$$\mathit{\delta P}\left(x,t\right)=\mathrm{\Delta }{P}_O\left(\mathit{erfc}\left(\frac{x}{2\sqrt{\left(Dt\right)}}\right)-\mathit{erfc}\left(\frac{x}{2\sqrt{\left(D\left(t-\mathit{\delta t}\right)\right)}}\right)\right).$$E28

At times t≫Δt, the pressure variation δP is given to a good approximation by

$$\mathit{\delta P}\left(x,t\right)=\frac{\mathrm{\Delta }{P}_{O}Dx\mathit{\delta t}exp\left(-x^2/\left(4Dt\right)\right)}{2\sqrt{\pi {\left(Dt\right)}^3}}.$$E29

The maximum pressure variation at distance x occurs after a time delay t_D given by

$$t_D=\frac{x^2}{6D}.$$E30

It follows that the maximum pressure variation at distance x and time t_D is given by δP_M where

$$\delta P_M\left(x,t\right)=\frac{3\sqrt{6}\mathrm{\Delta }{P}_{O}D\mathit{\delta t}}{\sqrt{\left(\pi e^3\right)}{x}^2}.$$E31

These results, in Eqs. (30) and (31), can be compared with those for the radially symmetric case in Eqs. (17) and (18). In both cases, t_D is proportional to the square of distance and inversely proportional to D, but differs by a numerical factor. Alternative empirical analysis by Hettema et al. [94], based on ‘rules of thumb’ rather than derivation from first principles, predicts a value for t_D that likewise differs by a numerical factor, but does not predict pressure variations. The pressure variation varies inversely with the square of distance for both the radially symmetric case and for the one-dimensional case (cf. Eqs. (18) and (31)).

Δx can thus be calculated using the exact formula for δP (Eq. (28)) as

$$\mathrm{\Delta }x=\frac{4\gamma \alpha \mathrm{\Delta }{P}_O}{B_E}\sqrt{\left(\left(\frac{Dt}{\pi }\right)-\left(\frac{D\left(t-\mathit{\delta t}\right)}{\pi }\right)\right)},$$E32

and using the approximate formula (Eq. (29)) as

$$\mathrm{\Delta }x=\frac{\mathit{\gamma \alpha }\mathrm{\Delta }{P}_O\mathit{\delta t}}{B_E}\sqrt{\left(\frac{D}{\pi t}\right)}.$$E33

Segall [87] reported that, for the Δσ_kk = 0 boundary condition applicable for this analysis,

$$D\equiv \frac{k}{\eta }\frac{2\mu \left(1-\nu \right)}{\left(1-2\nu \right)}\frac{B\left(1+\nu \right)}{3\alpha \left(1-\nu \right)-2B{\alpha }^2\left(1-2\nu \right)},$$E34

where μ is the shear modulus of the rock and B is its Skempton coefficient.

One may likewise quantify the loss of volume ΔV as a result of the vertical compaction Δz of the Dinantian limestone. By analogy with Eq. (24), an upper bound to Δz can be estimated as

$$\mathrm{\Delta }z={\int }_{x=0}^{x=H}\frac{\left(1-\gamma \right)\alpha \mathrm{\Delta }{P}_O}{B_E}\mathit{dz},$$E35

where H is the thickness of the Dinantian limestone. If the volume of limestone thus affected has dimensions L_y parallel to and L_x perpendicular to the model fault, so ΔV = L_x × L_y × Δz or

$$\mathrm{\Delta }V=\frac{\left(1-\gamma \right)\alpha \mathrm{\Delta }{P}_{O}H{L}_x{L}_y}{B_E}.$$E36

4.3 Mohr-Coulomb failure analysis

The tendency for coseismic slip on the seismogenic fault is analysed using the standard Mohr-Coulomb approach. The Mohr-Coulomb failure parameter Φ:

$$\mathrm{\Phi }=\tau -c\left({\sigma }_N-P\right),$$E37

will thus be evaluated where σ_N, τ and c are the resolved normal stress, shear stress and coefficient of friction on the fault plane, and P is the fluid pressure in the fault. Φ = 0 marks this condition, with Φ < 0 indicating frictional stability. This analysis can also be visualised using the standard Mohr circle construction, as a graph of τagainst effective normal stress σ_N’, defined as σ_N − P (see below).

From Eq. (37), other factors remaining constant, a decrease in P will ‘clamp’ a fault, making it more stable, and an increase in P will ‘unclamp’ a fault, potentially resulting in seismicity. The latter change is the accepted mechanism for the widespread occurrence in recent years of seismicity caused by fluid injection (e.g., [95, 96, 97, 98]). On this basis, one might conclude that a decrease in the groundwater pressure within the Newdigate fault cannot be the cause of the local seismicity.

However, it is generally accepted that the mechanics of faults, notably whether they are stable or can slip seismically, are determined by the properties of strong patches – asperities – where the opposing fault surfaces are in frictional contact (e.g., [99]). A fault surface consisting, on a microstructural scale, of a fractal size distribution of asperities with a small proportion of the fault surface in contact, in proportion to the normal stress applied to the fault, can mimic the effect, on a macroscopic scale, of a constant coefficient of friction (e.g., [100, 101]). Brown and Scholz [102] showed that natural rock surfaces follow fractal scaling for surface features of height up to ∼0.1 m. Laboratory simulations of faulting often include asperities on a microstructural scale (e.g., [103, 104]). Most recently, the view has gained ground that the physics of co-seismic faulting is likewise governed by processes on a microstructural scale (e.g., [105, 106]). For example, McDermott et al. [106] deduced that asperities can be patches of fault with areas of no more than a few square metres, their properties being determined by mineral grains with dimensions of microns. Because these strong patches with fault surfaces in contact occupy only a small proportion of a fault surface, they act as stress concentrations. For example, in the laboratory experiments by Selvadurai and Glaser [104], millimetre-sized asperities with micron-sized heights occupy a very small proportion of the fault area; in one experimental run, a decrease in the mean normal stress across the fault area by ∼0.3 MPa caused decreases in the normal stress affecting individual asperities by ∼10 MPa.

Figure 5(b) illustrates how a small increase in separation of fault surfaces, Δx, can destabilise a fault through its effect on contact between asperities. Moving from configuration (i) to configuration (ii), two of the three asperities depicted will no longer contribute to fault stability. The third one will experience a significantly reduced normal stress, as a result of the increased separation of the fault surfaces. This will reduce the maximum shear stress that this asperity can sustain, in accordance with Eq. (37), whereas the shear stress that it is required to sustain to keep the fault stable will increase because the other asperities no longer contribute. Overall, it can thus be seen how a small increase in separation of fault surfaces might bring a fault significantly closer to the condition for slip, and might indeed result in coseismic slip.

As others (e.g., [107, 108]) have noted, a general calculation of this ‘unclamping’ effect for a fault with a general size-distribution of asperities would be very complex; this is thus not attempted here. A simplified calculation is instead presented, assuming that a patch of fault is prevented from slipping by a single asperity. For this patch, of area A, the normal stress and shear stress areσ_N and τ; the single asperity has cross-sectional area a, uncompressed height b, and Young’s modulus E. The effect of σ_N compresses this model asperity to height d (Figure 11). The tip of this asperity will act as a stress concentration, with normal stress and shear stress(A / a) σ_N and (A / a) τ and coefficient of friction c. The areas of fault where the wall rocks are not in contact initially contain fluid with pressure P_O (Figure 11(a)), the fault being stable with Φ = -ΔΦ. The reduction in fluid pressure to P_O -ΔP_O is followed by poroelastic separation of the fault wall rocks by distance Δx, which brings the fault to the condition for slip (Φ = 0). One may thus state versions of Eq. (37) at the tip of the model asperity for these ‘before’ and ‘after’ conditions:

Figure 11.

Schematic model asperity used to calculate effects of poroelastic fault unclamping. (a) Initial state, with fluid pressure within the fault P_O and the asperity (of uncompressed height b and cross-sectional area a) compressed to height d by the normal stress across the fault. (b) Modified state, with fluid pressure within the fault reduced to P_O-ΔP_O and the asperity compressed to height d-Δx/2 as a result of the separation Δx between the wall rocks on both sides of the fault caused by their poroelastic compaction.

$$\tau \frac{A}{a}-c\left(E\frac{\left(b-d\right)}{b}-P_O\right)=-\mathrm{\Delta \Phi },$$E38

and

$$\tau \frac{A}{a}-c\left(E\frac{\left(b-d-\mathrm{\Delta }x/2\right)}{b}-P_O+\mathrm{\Delta }{P}_O\right)=0.$$E39

Subtraction of Eq. (38) from Eq. (39) gives

$$\mathrm{\Delta \Phi }=-c\left(\mathrm{\Delta }{P}_O-\mathrm{\Delta }{\sigma }_N\right)=-c\left(\mathrm{\Delta }{P}_O-\frac{E\mathrm{\Delta }x}{2b}\right),$$E40

or, substituting Δx from Eq. (27),

$$\mathrm{\Delta \Phi }=c\mathrm{\Delta }{P}_O\left(\frac{E\alpha }{b{B}_E}\sqrt{\left(\frac{Dt}{\pi }\right)}-1\right).$$E41

One may also combine Eqs. (27) and (40) by eliminating ΔP_O, to obtain

$$\mathrm{\Delta }x=\frac{2\mathrm{\Delta \Phi }}{c\left(\frac{E}{b}-\frac{B_E}{2\gamma \alpha }\sqrt{\left(\frac{\pi }{Dt}\right)}\right)}.$$E42

In the limit of high D t / b, at time t≫t_s, where

$$t_s=\frac{\pi b^2{B_E}^2}{4{\gamma }^2{\alpha }^2{E}^{2}D},$$E43

this equation simplifies to

$$\mathrm{\Delta }x\approx \frac{2b\mathrm{\Delta \Phi }}{cE},$$E44

indicating that this poroelastic unclamping effect requires Δx / b to be comparable to ΔΦ / E. Eq. (44) may also be rearranged, recalling from Eq. (40) that Δσ_N = E Δx / (2 b), to give

$$\mathrm{\Delta \Phi }=c\mathrm{\Delta }{\sigma }_N.$$E45

The geomechanical consequences of this model can now be illustrated using the Mohr circle construction (Figure 12), for a model fault with c = 0.6. A model stress field is adopted, similar to that deduced by Westaway [11] at 2400 m depth for the Preese Hall case study, with σ_H = 64.488 MPa, σ_V = 54.300 MPa, and σ_h = 37.880 MPa, with hydrostatic groundwater pressure P = 23.544 MPa. As Figure 12(a) indicates, the state of stress on a vertical model fault (characterised by effective normal stress of magnitude σ_N’ = σ_N-P = 22.332 MPa and shear stress τ = 12.199 MPa), with normal vector oriented at 57° to the maximum principal stress, plots below the Coulomb failure envelope, indicating stability, the differential stress Δσ being 26.608 MPa. If P within this fault decreases by 10 kPa to 23.534 MPa and this causes, via the poroelastic mechanism described above, in a reduction in the magnitude of σ_N’ by 2 MPa to 20.332 MPa, the resulting state of stress is depicted in Figure 12(b). This condition, with τ still 12.199 MPa, now plots on the Coulomb failure envelope, indicating instability. The fault normal vector is now oriented at 60° to the maximum principal stress, reflecting the slight rotation of the stress field in the vicinity of the fault, caused by the poroelastic reduction in σ_N, which accompanies an increase in Δσ to 28.260 MPa. The state of stress on the model fault has thus effectively shifted left by 2 MPa on the Mohr-Coulomb plot, moving towards the failure envelope by 1.2 MPa, these adjustments being interrelated via Eq. (45) given c = 0.6. This poroelastic adjustment to the stress field thus involves reducing the magnitude of σ_N’ keeping τ constant, rotating the near-fault stress field and increasing Δσ. It is different from what occurs during ‘fracking’ of impermeable rocks, where an increase in P causes a Mohr circle of constant diameter (indicating constant Δσ) to shift leftward until part of its circumference touches the failure envelope.

Figure 12.

Mohr circle diagrams representing a model stress field at 2400 m depth, to indicate the sense of changes to the state of stress envisaged as causing the Newdigate seismicity. (a) for a model stress field with σ_HO = 64.488 MPa, σ_VO = 54.300 MPa, and σ_hO = 37.880 MPa, with σ_H oriented at 57° to the normal direction to the fault. Hydrostatic groundwater pressure P = 23.544 MPa causes σ_H^’_O = 40.944 MPa, σ_V^’_O = 30.756 MPa, and σ_h^’_O = 14.336 MPa, resulting in σ_L^’_O = 28.723 MPa and σ_M^’_O = 27.706 MPa. The resolved shear stress and normal stress on the model fault, τ = 12.199 and σ_N’_O = 22.332 MPa, plot below the coulomb failure line for c = 0.6, with (from Eq. (37)) Φ = -1.2 MPa, indicating that the fault is stable under these conditions. (b) for revised conditions consistent with the set of processes in Fig. 5. Groundwater pressure adjusts by ΔP = -10 kPa, to P_f = 23.534 MPa, and the principal stresses adjust toσ_H = 65.130 MPa and σ_h = 36.870 MPa keeping σ_V = σ_VO = 54.300 MPa, with σ_H now oriented at 60° to the fault normal. As a result, σ_H^’ = 41.596 MPa, σ_V^’ = 30.766 MPa, and σ_h^’ = 13.336 MPa, resulting in σ_L^’ = 28.566 MPa and σ_M^’ = 27.466 MPa. The resolved shear stress and normal stress on the fault, τ = 12.199 and σ_N’ = 20.332 MPa, now plot on the coulomb failure line for c = 0.6, with Φ = 0, indicating that the fault is now frictionally unstable. The calculated increase in Φ by 1.2 MPa, for Δσ_N = 2 MPa with c = 0.6, is consistent with Eq. (45).

4.4 Estimation of model parameters

Application of the above theory to the Newdigate case study requires determination of many model parameters, representing properties of key lithologies (Portland Sandstone, calcite ‘beef’, and Dinantian limestone) and of the Newdigate fault. To facilitate this, it is noted that the elastic moduli that appear in the foregoing analysis are interrelated via standard formulas, such as

$$E\equiv 3B_R\left(1-2\nu \right)$$E46

and

$$\mu \equiv 3B_R\frac{1-2\nu }{2\left(1-\nu \right)}.$$E47

Hydraulic conductivity K and permeability k are also interrelated thus:

$$K\equiv \frac{k{\rho }_{W}g}{\eta }$$E48

where g is the acceleration due to gravity. Different formulas for hydraulic diffusivity D, subject to different boundary conditions, have already been noted. However, the considerable uncertainty in model parameters for lithologies in the present study region makes such distinctions moot; D will, therefore, be estimated using Eq. (5).

For the Upper Portland Sandstone at Brockham, Angus [27] reported a 3 m thick layer with ϕ0.25 and k 200 mD. Lee [109] reported B_R 13 GPa as typical for dry sandstone with ϕ0.25. Taking B_W 2.15 GPa for water (from [110]), using Eq. (4), B_E for Portland Sandstone with pore space occupied by water is ∼6 GPa. At Horse Hill, Xodus [26] reported a 35 foot or 11 m section in Portland Sandstone with permeability up to 20 mD. The water in contact with the Portland reservoirs at both sites, at ∼600 m depth, is at ∼25 °C (cf. [111]), for which η = 0.9 mPa s [112], with ρ_w 1000 kg m⁻³ and g 9.81 m s⁻². Under these conditions the oil at Brockham has η = 11 cP or 11 mPa s [27]. Using Eq. (5), with k 200 mD, B_E ∼ 6 GPa, and η 0.9 mPa s, D for Portland Sandstone is ∼1.3 m² s⁻¹, which can be rounded to 1 m² s⁻¹.

Many workers (e.g., [42, 113, 114]) have investigated the aperture, or width, of bedding-parallel fractures (typically filled with ‘beef’) in shale, as a guide to its hydraulic properties. In a study spanning several shale provinces, Wang [113] found fractures with width between 15 μm and 87 mm. Many of the wider ones could be seen to have opened by multiple increments, each adding a few tens of microns of width, prior to cementation due to growth of calcite. Permeability and fracture aperture can be interrelated by comparing the Darcy equation for laminar fluid flow, Q = (k A / η) dP/dx, and the Poiseulle equation for laminar flow between parallel boundaries, Q = (D W² / (12 η)) dP/dx, which is a solution to the more general Navier–Stokes equation for fluid flow (e.g., [115]). Here Q is the volume flow rate, η the viscosity of the fluid, dP/dx the pressure gradient in the direction of flow, k the permeability of the medium, A the cross-sectional area of the flow, and W and D the width of the channel and its length in the direction perpendicular to the flow. Combining these two formulae, equating A to D × W, gives k ≡ W²/12. This formula gives the permeability equivalent to W = 15 μm as ∼20 D (∼2 × 10⁻¹¹ m²). Overall, it is inferred that the ∼900 mD value, from Carey et al. [50], might be applicable to calcite ‘beef’ in the present study area. Identifying a suitable representative value for B_R, the bulk modulus for calcite ‘beef’, is problematic, because of its strongly anisotropic character. ‘Beef’ is abundant within mudstones of the Neuquen Basin of Argentina (e.g., [116]). Sosa Massaro et al. [117] estimated a representative vertical Young’s modulus and Poisson’s ratio for this lithology as ∼15 GPa and ∼ 0.25; using Eq. (46) these parameters yield a bulk modulus of ∼10 GPa. Using Eq. (5), with k 900 mD, K_B ∼ 10 GPa, and η 0.9 mPa s, D for calcite ‘beef’ can be estimated - subject to considerable uncertainty - as ∼10 m² s⁻¹. The calcite ‘beef’, reported in the BGS borehole viewer online documentation as ‘veins of secondary calcite’, in the CF1 borehole occurred at depths of 2327–2329 feet or 709.3–709.9 m TVD, thus ∼626 m TVDSS. Although this interval was cored, the core was not analysed for permeability; however, core between 2300 and 2301 feet TVD yielded k 1650 mD. The top Portland in this borehole is at 1753 feet or ∼ 534 m TVDSS, the estimated base of the oil reservoir at ∼566 m TVDSS (Figure 7); this ‘beef’ layer is thus ∼100 m below the top Portland. Based on this information this layer is assigned a nominal thickness h_B of 1 m for the purpose of modelling.

Carbonate rocks such as the Dinantian limestone are likely to be complex, being fractured, so water storage within them will be in part by opening of fractures and in part by opening of pore space. Dinantian limestone typically has low matrix porosity (e.g., [118]), its ability to store and transmit groundwater being largely via fractures. Parameter values for Dinantian limestone include B_R = 50 GPa and ϕ = 0.04, along with Young’s modulus E_R = 75 GPa and Poisson’s ratio ν = 0.25, from Bell [119]. With this set of values, B_E is ∼27 GPa, and α = 1–27/50 = 0.46 (Eq. (22)). Poisson’s ratio ν ranges between 0.19 and 0.31 [119] so 0.25 is adopted, for which μ ≡ B (Eq. (47)). No attempt is made here to determine the value of γ; a value of 0.5 will be assumed, consistent with depressurization causing closure of both vertical and horizontal fractures to an equivalent extent. Skempton’s coefficient B has been reported as ∼0.4 for many limestone samples (e.g., [120, 121, 122]). Bell [119] noted a range of values of K for laboratory samples of Dinantian limestone, ranging from 0.07 × 10⁻⁹ m s⁻¹ to 0.3 × 10⁻⁹ m s⁻¹. Lewis et al. [123] reported much higher values ranging from ∼10⁻⁶ m s⁻¹ to ∼10⁻² m s⁻¹ in karstified regions, which correspond (using Eq. (48)) to k ≥ 100 mD. Using the latter set of values, Eq. (2) gives values for D ranging upward from ∼3 m² s⁻¹. For comparison, Shepley [124] determined an upper bound to D for Dinantian limestone in the Peak District of central England by modelling the hydrology of Meerbrook Sough, a disused mine drainage adit that drains a ∼ 40 km² area. His analysis reported an upper bound of 50,000 m² day⁻¹ or ∼ 0.6 m² s⁻¹. However, this analysis did not reproduce the observed magnitude of seasonal fluctuations in flow, favouring a higher value of D. Overall, it is inferred that that D ∼ 1 m² s⁻¹ is appropriate for karstified Dinantian limestone, subject to considerable uncertainty. For comparison, Hornbach et al. [125] deduced that a poroelastic pressure pulse resulting from large-scale injection of waste water propagated for up to ∼40 km through the Ellenburger Formation, a karstified limestone of Ordovician age, in ∼6 years, resulting in earthquakes in the vicinity of Dallas, Texas. Using Eq. (30) gives an upper bound to D for the Ellenburger Formation of ∼1.4 m² s⁻¹, in reasonable agreement. Zhang et al. [126] reported a nominal value of 1 m² s⁻¹ for this karstified Ordovician limestone.

For an ensemble of faults in different lithologies, Brodsky et al. [107] deduced that the typical asperity height b in length L of a fault scales as

$$b=b_O{L}^{\zeta },$$E49

where ζ = ∼0.6 and b_O = ∼10⁻³ m^0.4. For example, with L = 100 m, Eq. (49) gives b ∼ 16 mm.

To investigate the area of the patch of fault that slipped in each earthquake, and to thus quantify the associated asperity height, seismic moment M_O was first determined using the standard formula

$${log}_{10}\left(M_o/Nm\right)=9.05+1.5M_W$$E50

[127], where M_W is moment magnitude. Next, the radius a of the equivalent circular seismic source was determined from M_o, assuming a nominal coseismic stress drop δσ:

$$M_o=\frac{16\mathit{\delta \sigma }\left(1-\nu \right)a^3}{3\left(2-\nu \right)}$$E51

(e.g., [128]). The corresponding diameter 2a is equated to fault length L to characterise asperity height. One thus obtains

$$L=2a={\left(\frac{3\left(2-\nu \right)M_o}{2\mathit{\delta \sigma }\left(1-\nu \right)}\right)}^{1/3}.$$E52

4.5 Application of model: Horse Hill

The first topic to be assessed is the possibility that the resumption in February 2019 of production from the Portland reservoir in well HH1 caused the third ‘burst’ of Newdigate seismicity. This phase of production is inferred to have begun at 08:00 GMT on 11 February 2019 and caused earthquakes starting with the event (M_W 2.5; M_O 6.8 × 10¹² N m) at 07:43 GMT on 14 February 2019 (located at TQ 22959 41543, after [1], at a probable depth of ∼2400 m), event H in Figure 7. The time delay in this instance was thus 3 days or ∼ 72 hours, for a separation of ∼3 km (Figures 1,7), the HH1 well bottom being ∼2.85 km from nearest point where the top Portland sandstone is transected by the Newdigate Fault and ∼ 3.0 km from an updip projection of the hypocentre to this cutoff. This production took place at ∼220 barrels per day or ∼ 0.40 l s⁻¹. The pressure drawdown during this phase has not been reported. However, as detailed in the online supplement, during the flow testing of the Portland reservoir in well HH1 in July–August 2018, the developer reported production at a sustained rate of ∼190 barrels per day or ∼ 0.35 l s⁻¹ (∼140–160 barrels per day being of oil), with stable bottom hole pressures ∼1.4 MPa below the initial reservoir pressure of ∼6.3 MPa; scaling in proportion gives a ∼ 1.6 MPa pressure drawdown at the HH1 well bottom in February 2019 (i.e., ∼1.4 MPa × 220/190). Using Eq. (52), with ν = 0.25 and δσ = 10 MPa, gives L = 134 m; using Eq. (49) gives b = 18.9 mm.

Figure 13(a) shows the predicted pressure variation in the ‘beef’ layer that is inferred to hydraulically connect the Horse Hill Portland reservoir and the Newdigate fault. With the chosen parameter values, the resulting pressure decrease at 3 km distance is calculated using Eq. (8) as ∼6.5 kPa after 2.5 days and ∼ 8.7 kPa after 3 days. Taking h_P = 11 m, η = 0.9 mPa s, D = 1.3 m² s⁻¹, and k_P = 35 mD (somewhat higher than the 20 mD reported by Xodus [26], with r_W = 88.9 mm (appropriate for 7-inch diameter casing), using Eq. (14) the pressure decline at the well bottom, 2–4 weeks after the start of production, would be ∼1.6 MPa, as expected. If η were adjusted to 11 mPa s, to reflect oil rather than water, k_P ∼ 430 mD would be required to match this pressure variation.

Figure 13.

Modelling of hydraulic consequences of the phase of production from well HH1 starting in February 2019. (a) Graphs of -ΔP versus radial distance r within the layer of calcite ‘beef’ at ∼600 m depth, which is inferred to connect the Horse Hill oil reservoir with the Newdigate fault, at different times t after the start of production. Calculations use Eq. (8) with Q = 0.4 l s⁻¹, η 0.9 mPa s, h_B 1 m, and D_B 10 m² s⁻¹. Dashed lines indicate, for r = 3 km, ΔP = -6.5 kPa after t = 2.5 days and ΔP = -8.7 kPa after t = 3 days. (b) Graphs of -δP/ΔP_O versus distance x perpendicular to the Newdigate fault, at different times after the pressure variation in (a) reached this fault. Calculations use Eq. (20) with D_D 1 m² s⁻¹, for Dinantian limestone. (c) Graph of Δx versus time for the pressure variations depicted in (b). Calculations use Eq. (27) with ΔP_O 8.7 kPa (from (a)), and B_E 27 GPa, α = 0.46, and D again 1 m² s⁻¹, for Dinantian limestone. Dashed line indicates Δx = 0.00504 mm, which arises after 0.25 hours or 15 minutes. (d) Graph of ΔΦ for the patch of the Newdigate fault that slipped in the 14 February 2019 earthquake versus time for the variations in Δx depicted in (c). Calculations use Eq. (44), with the same parameters as for (c) plus b 18.9 mm, c 0.6, and E 75 GPa for the model asperity on the seismogenic patch of fault. Dashed line indicates ΔΦ = 6 MPa, which arises after 0.25 hours or 15 minutes.

Figure 13(b) shows the variation in fluid pressure inside the Dinantian limestone alongside the seismogenic part of the Newdigate fault, as a result of the imposed pressure variation in the ‘beef’. Figure 13(c) shows the corresponding variation in Δx, the poroelastic separation of the opposing faces of this fault, and Figure 13(d) shows the corresponding variation in ΔΦ. The upper limit to ΔΦ is taken as cδσ, or 6 MPa, the coseismic stress drop, δσ = 10 MPa, also being the increase in shear stress during a complete earthquake cycle; this value of ΔΦ corresponds (using Eq. (44)) to Δx = 0.005 mm. These conditions are taken as indicating the condition for an earthquake to occur on the specified patch of fault. It is evident that with the specified combination of parameter values, these changes will occur very rapidly, in just 15 minutes after the idealised step onset of the pressure variation in the fault. In reality, Figure 13(a) indicates a gradual onset of this pressure variation, rather than an abrupt step, so the calculated changes Δx and ΔΦ will occur gradually, leading to an earthquake at the calculated time, roughly three days after the start of production from the well. This modelling demonstrates that, with the parameter values adopted, production from the Portland reservoir by well HH1 has been readily able to cause seismicity on the Newdigate fault.

4.6 Application of model: Brockham

The second topic assessed will be the possibility that the ∼4 m³ of production from well BRX2Y on 23 March 2018, inferred to have started at 08:00 BST, caused earthquakes starting with the event (M_W ∼ 2.8; M_O ∼ 1.3 × 10¹³ N m) at 12:10 BST on 1 April 2018 (located at TQ 21992 41976, after [1], also at a probable depth of ∼2400 m), event B in Figure 7, with hypocentre close to the western or WNW limit of the seismogenic part of the Newdigate fault. The time delay in this instance was 9 days ∼4 hours or ∼ 220 hours. In a straight line, the BRX2Y well bottom is ∼7.2 km NNW of the point, updip from this hypocentre, on the top Portland cutoff of the Newdigate fault (between circa TQ 18440 48310 and circa TQ 22100 42200). However, as already noted, it is possible that the high throw on the Brockham Fault along this direct path blocks any high-permeability connection, in which case the actual path length between these end points, through localities farther west with low fault offsets (Figure 8), and around the western end of the Leigh Fault (Figure 1), might be >1 km longer. Using Eq. (52), with ν = 0.25 and δσ = 10 MPa, M_O ∼ 1.3 × 10¹³ N m for the 1 April 2018 earthquake gives L = 166 m; using Eq. (49) gives b = 21.5 mm.

Figure 14(a) shows the predicted variation in fluid pressure in the ‘beef’ layer that is inferred to hydraulically connect the Brockham Portland oil reservoir and the Newdigate fault, calculated using Eq. (15). With the chosen set of parameter values, following 4 m³ of production the resulting pressure decrease at r = 8 km distance is predicted to be ∼50 Pa after 9 days for D_B = 10 m² s⁻¹. Assuming h_P = 3 m, η = 0.9 mPa s, D_P = 1.3 m² s⁻¹, and k_P = 200 mD, with r_W = 88.9 mm (again appropriate for 7-inch diameter casing), using Eq. (14) the pressure decline at the well bottom at the end of production would be ∼530 kPa if the production took 4 hours or ∼ 1.0 MPa if completed in 2 hours. Within an hour of this idealised ∼50 Pa step pressure variation affecting the Dinantian limestone, Δx would be ∼0.06 μm (Figure 14(b)). This would be sufficient to unclamp the Newdigate fault by ∼60 kPa or ∼ 1% of δσ (Figure 14(c)). With D_B adjusted to 20 m² s⁻¹, the pressure decrease at r = 8 km would peak, circa 9 days after this pulse of production, at ∼150 Pa (Figure 14(d)). Within three quarters of an hour of this idealised ∼150 Pa step pressure variation affecting the Dinantian limestone, Δx would be ∼0.14 μm (Figure 14(e)). This would be sufficient to unclamp the Newdigate fault by ∼150 kPa or ∼ 2.5% of δσ (Figure 14(f)). Evidently, for such a small change in the state of stress to have caused the observed seismicity, the patch of this fault that slipped on 1 April 2018 must have already been within a very small proportion of δσ of the Mohr-Coulomb condition for slip.

Figure 14.

Modelling of hydraulic consequences of the pulse of production from well BRX2Y on 23 march 2018. (a) Graphs of -ΔP versus radial distance r within the layer of calcite ‘beef’ at ∼600 m depth, which is inferred to connect the Brockham oil reservoir with the Newdigate fault, at different times t after the start of the pulse of production from well BRX2Y on 23 march 2018, time t being measured from the start of production. Calculations using Eq. (15) assume 4 m³ produced volume, calculated as Q 1.39 × 10⁻⁴ m³ s⁻¹ for Δt 8 hours, η 0.9 mPa s, h_B 1 m, k_B 900 mD, and D_B 10 m² s⁻¹. Dashed line indicates ΔP = -51 Pa for r = 8 km after t = 9 days. (b) Graph of Δx versus time following a step pressure reduction within the Newdigate fault. Calculations use Eq. (27) with ΔP_O = 51 Pa (cf. (a)), with K_D 27 GPa, and D_D again 1 m² s⁻¹, for Dinantian limestone. Dashed line indicates Δx = 0.057 μm, which arises after 0.95 hours or 57 minutes. (c) Graph of ΔΦ for the patch of the Newdigate fault that slipped in the 1 April 2018 earthquake versus time for the variations in Δx depicted in (b). Calculations use Eq. (44), with the same parameters as for (b) plus b 21.5 mm, c 0.6, and E 75 GPa for the model asperity on the seismogenic patch of fault. Dashed line indicates ΔΦ = 60 kPa, which arises after 0.95 hours or 57 minutes. (d) Graphs as for (a), except D_B is now 20 m² s⁻¹. Dashed line now indicates Δ_P = -146 Pa for r = 8 km after t = 9 days. (e) Graph as for (b), except D_B is now 20 m² s⁻¹. Dashed line now indicates Δx = 0.143 μm, which arises after 0.72 hours or 43 minutes. (f) Graph as for (c), except D_B is now 20 m² s⁻¹. Dashed line indicates ΔΦ = 150 kPa, which arises after 0.72 hours or 43 minutes.

4.7 Application of model: interference between the wells

The final topic assessed will be pressure interference between the wells. Like the path from the Brockham well to the Newdigate fault, the most permeable connection, through ‘beef’, between the two wells will exceed the straight line distance; separation r = 10 km is adopted. Again calculated using Eq. (29), with D_B = 10 m² s⁻¹ the 4 m³ of production from well BRX2Y would cause a maximum pressure decline in the vicinity of well HH1 of ∼50 Pa after a ∼ 1 month delay (Figure 15(a)). With D_B increased to 20 m² s⁻¹, this production would cause a maximum pressure decline of ∼90 Pa after a ∼ 15 day delay (Figure 15(b)). Variations in HH1 bottom hole pressure of this order, developing and dissipating on timescales of weeks or months in response to attempts at production from well BRX2Y, would be extremely difficult, if not impossible, to recognise given the >1 MPa pressure variations caused by the production from this well. This is consistent with the statement by the OGA [5] that no pressure variation at HH1 was detectable in response to the pulses of production from BRX2Y.

Figure 15.

Modelling of pressure interference of production from well BRX2Y on well HH1. (a) Graphs of pressure variations -Δ_P in ‘beef’ adjoining well HH1 caused by the pulse of production from well BRX2Y on 23 March 2018, time t being measured from the start of production. Calculations using Eq. (15) assume 4 m³ produced volume, calculated as Q 1.39 × 10⁻⁴ m³ s⁻¹ for Δt 8 hours, η 0.9 mPa s, h_B 1 m, k_B 900 mD, and D_B 10 m² s⁻¹. Horizontal dashed line indicates Δ_P = -47 Pa for r = 10 km after t = 35 days. (b) Graphs as for (a), except D_B is now 20 m² s⁻¹. Horizontal dashed line now indicates Δ_P = -93 Pa for r = 10 km after t = 35 days. (c) Graphs illustrating the pressure drawdown caused by 20 years of production from well BRX2Y. Calculations using Eq. (8) assume Q 6 × 10⁻⁵ m³ s⁻¹, η 0.9 mPa s, and h_B 1 m, for variable k_B and D_B. For D_B = 10 m² s⁻¹, k_B = 900 mD; for other values of D_B, k_B is adjusted in proportion (cf. Eq. (5)). Horizontal dashed line now indicates ΔP = -23.7 kPa for r = 10 km with D_B = 10 m² s⁻¹. Using Eq. (8), with k_P = 20 mD, h_P = 2 m, r_w = 0.0889 m, and D_P = 1 m² s⁻¹, gives a predicted bottom-hole pressure decline at well BRX2Y of ∼2.9 MPa, roughly as observed.

A second test is possible, given that ∼20 years of production at Brockham (followed by two years of shut-in, during 2016–2018) resulted in a bottom-hole pressure of ∼3.4 MPa, roughly 3 MPa below hydrostatic, whereas the initial bottom-hole pressure in well HH1 was ∼6.3 MPa, roughly hydrostatic. The pressure drawdown at Brockham evidently had no significant effect on the pressure at HH1. It might thus be inferred that the two wells are not hydraulically connected, in contrast with the arguments in the present study. To explore this issue, Eq. (8) is used to calculate the effect of twenty years of production at the steady rate of 6 × 10⁻⁵ m³ s⁻¹ to obtain the 36,900 m³ produced (Figure 15(c)). Taking into account the dependence of D_B on permeability, it is evident that for D_B ∼ 10–20 m² s⁻¹, as envisaged in the present study, the predicted pressure decline at HH1 was only a few tens of kilopascals, thus not significant in relation to the measured pressure. The large difference in bottom-hole pressure between the two wells in 2018 is thus not evidence against the proposed model.

The production test from well HH1 that started on 11 February lasted until late June 2019, some 140 days or 20 weeks, albeit with some intermittency (see supplement). Using Eq. (8), the resulting pressure decline is depicted in Figure 16(a) for D_B = 10 m² s⁻¹ and in Figure 16(b) for D_B = 20 m² s⁻¹. At r = 10 km, in the vicinity of well BRX2Y, the maximum pressure decline is ∼40 kPa (Figure 16(a)) or ∼ 60 kPa (Figure 16(b)). The production from well HH1 thus influenced bottom hole pressure in well BRX2Y by three orders of magnitude more than for the effect of well BRX2Y on well HH1. The OGA [6] made no mention of any pressure data for well BRX2Y that might be used to test this deduction.

Figure 16.

Modelling of pressure interference of production from well HH1 on well BRX2Y. (a) Graphs of pressure variations -Δ_P in ‘beef’ adjoining well BRX2Y caused by the phase of production from well HH1 starting in February 2019, time t being measured from the start of production. Calculations using Eq. (8) assume Q 4 × 10⁻³ m³ s⁻¹, η 0.9 mPa s, h_B 1 m, k_B 900 mD, and D_B 10 m² s⁻¹. Horizontal dashed line indicates ΔP = -37.8 kPa for r = 10 km after t = 140 days. (b) Graphs as for (a), except D_B is now 20 m² s⁻¹. Horizontal dashed line now indicates ΔP = -56.6 kPa for r = 10 km after t = 140 days.

5. Discussion

The proposed mechanism for the Newdigate seismicity depends on a pressure drop within the Dinantian limestone alongside the seismogenic strand of the Newdigate fault zone, as a result of depressurization of the water within this fault, caused by oil production from neighbouring wells (Figure 5). For the production from BRX2Y to have caused seismicity by this mechanism, the seismogenic fault must already have been extremely close (maybe within ∼60 kPa; see above) to the Mohr-Coulomb failure condition. It can be inferred that the same mechanism, operating during the previous production from this well, contributed to creating this state of stress by progressively depressurizing the Dinantian limestone. To test this possibility, one may use Eq. (36), noting the 2.7 MPa depressurization of the reservoir over ∼20 years and estimating from the previous analyses (e.g., Figure 15(d)) that the resulting value of ΔP (and, thus, δP) within this limestone would be ∼10 kPa. With B_E = 27 GPa, α = 0.46, H = 70 m, andΔV = 36,900 m³ to balance the production, depressurization of a ∼ 6000 km² area would be indicated; if roughly equidimensional, this would have a radius of ∼40 km. However, in reality, it is to be expected that such depressurization would be largely cancelled by recharge of water into the Dinantian limestone from other directions, which is not incorporated into the model. For example, if after this cumulative production, x_M were ∼ 6 km and ΔP_O were ∼ 1 kPa, then substituting Eq. (21) into Eq. (41), and taking b ∼ 0.1 m (obtained for L = 1.5 km, the length of the seismogenic part of the Newdigate Fault, using Eq. (49)), ΔΦ would be ∼6 MPa. Notwithstanding the approximations made in the model, it is thus indeed plausible that the cumulative production at Brockham brought this fault to the condition for shear failure, assuming that it was already critically stressed before this production began.

In principle, testing of the proposed mechanism is possible, given the predicted vertical compaction in the Dinantian limestone (Eq. (35)). Such compaction will cause subsidence of the Earth’s surface, and so is in principle observable using multiple techniques, including interferometric synthetic-aperture radar (InSAR) and repeated gravity and GPS measurements. A combined dataset of this type has been analysed for a region of southeast England, including the northern Weald Basin, by Aldiss et al. [129]. At the October 2018 workshop attention was also drawn to an InSAR-derived surface deformation map of Britain by GVL [130], spanning October 2015 to October 2017. However, the predicted subsidence, resulting from the two decades of production at Brockham, will be only a small fraction of 1 mm (∼36,900 m³ / ∼200 km² ≈ 0.2 mm), even if none of the fluid withdrawal from this limestone were recharged. The Aldiss et al. [129] analysis revealed vertical crustal motions at ∼1 mm a⁻¹, caused by processes such as extraction from or recharge of shallow groundwater reservoirs. Such rates make it impossible to resolve the much smaller effect expected from compaction of the Dinantian limestone at Newdigate.

Much has been made by participants in the OGA [5] workshop regarding the extent to which the Newdigate earthquake ‘swarm’ might fit the standard criteria identified by Davis and Frohlich [131] for establishing whether instances of seismicity are anthropogenic (e.g., [54]). UKOG [8] have argued that this set of criteria is inapplicable as they relate to seismicity caused by fluid injection, which is not the causal mechanism in this case. However, familiarity with the literature in this field (e.g., [132]) indicates that these criteria are widely used irrespective of the geomechanical cause of any particular anthropogenic earthquake. Verdon et al. [6] proposed a different approach to assessing anthropogenic seismicity. This approach appears problematic, since it replaces the objective (yes / no) criteria recommended by Davis and Frohlich [131] with subjective numerical scores. The development of a conceptual geomechanical model for the Newdigate earthquakes supersedes the other Davis and Frohlich [131] criteria; nonetheless, an appraisal of this seismicity in terms of these criteria is included in the online supplement.

Although fluid injection is nowadays recognised as a widespread cause of induced seismicity, it is worth noting that reductions in fluid pressure caused by fluid extraction has been linked to seismicity for far longer. The first instance recognised is by Pratt and Johnson [133], for earthquakes accompanying oil production near Houston, Texas. Other case studies subsequently recognised include those by Calloi et al. [134], Kovach [135], Rothé and Lui [136], Simpson and Leith [137], Pennington et al. [138], Wetmiller [139], Grasso and Wittlinger [140], Doser et al. [141], Ottemöller et al. [142], Dahm et al. [143], and Hornbach et al. [144], whereas works discussing physical mechanisms for such seismicity, including compaction of reservoir rocks, include those by Yerkes and Castle [145], Simpson et al. [146], Segall [147], Segall and Fitzgerald [88], Ottemöller et al. [142], and Soltanzadeh and Hawkes [90]. Some of the above works (e.g., [138]) have recognised the significance of processes required in the present conceptual model (e.g., compaction of limestone and failure of asperities), and others (e.g., [148, 149, 150]) have recognised that highly permeable connections can cause seismicity at significant distance from the source of the causative change in fluid pressure. However, no previous case study known to the present author has proposed a geometry between the source of depressurization and the seismogenic fault that resembles the present conceptual model (Figure 5).

Hicks et al. [1] considered and rejected the possibility that the Newdigate seismicity was caused by compaction, as a result of depressurization caused by oil production, from one of the oil reservoirs in the area. They reached this conclusion on the basis of the strike-slip focal mechanisms (Figure 1), because in their view compaction would be expected to cause dip-slip earthquakes; they thus argued that these strike-slip earthquakes must be natural. They cited in support of this conclusion work on the Groningen seismicity in the Netherlands, where gas field depletion has caused many earthquakes, almost all with normal-faulting focal mechanisms (e.g., [151]). Compaction was initially thought to be the cause of the Groningen seismicity [152]. However, these events were later reinterpreted as caused by the combined effects of compaction and poro-elastic changes to the state of stress on faults [153], the normal faulting focal mechanisms thus reflecting the extensional stress state in the Netherlands. There is therefore no contradiction between this work and the occurrence at Newdigate of strike-slip earthquakes, the focal mechanism orientation again consistent with the local state of stress. Differences between the Newdigate and Groningen case studies concern, for the latter, the much larger scale (reservoir area ∼ 900 km²; surface subsidence ∼30 cm; cumulative production ∼10¹⁰ m³ at reservoir pressure; [154]), and the lower hydraulic diffusivity and elastic moduli (D 0.38 m² s⁻¹ and E ∼ 8–25 GPa; [153, 155]). Regarding the geomechanics, the principal differences are that the analysis for Groningen has neglected consideration of fault asperities and has assumed that the elastic moduli of the fault wall rocks vary over time following the start of compaction, relaxing from short timescale to long timescale values (this effect being assumed to occur over a characteristic timescale, specified as 7.3 years for no clear reason other than to make the model work). In the present study, the time-dependence of the response develops as a result of the time required for poroelastic compaction in the Dinantian limestone to create fault opening Δx (Eq. (24)) comparable to the typical fault asperity height b. This seems a more physically realistic approach, avoiding any arbitrary timescale parameter and in keeping with modern ideas (noted above) that coseismic slip on faults is governed by interactions of asperities. Nonetheless, as is evident, the present analysis makes many simplifying assumptions, not least regarding the geometry of the fluid flow and its variations, for example: recharge of the Dinantian limestone is neglected; and smooth variations in the pressure of the fluid being drawn from the Newdigate Fault (as in Figure 13(a)) are approximated as step changes (as in Figure 13(b)). Furthermore, as Segall [87] noted, the choice of boundary conditions for analysis of the combined pressure and strain response in the Dinantian limestone is one end member of a range of possibilities. The present analysis is essentially a ‘proof-of-concept’, to demonstrate that it is plausible that oil production caused the Newdigate seismicity and to shed light on the combination of hydraulic properties that makes this possible; it does not, of course, prove that the production caused this seismicity.

The remainder of this discussion will concentrate on geomechanical issues. Following each of the Newdigate earthquakes, the spatially averaged shear stress on the patch of fault that slipped will reduce by a value equal to the coseismic stress drop Δσ, moving the state of stress away from the Coulomb failure condition. The decline in seismicity in late 2019 suggests that it had reached some form of limiting state; the conditions governing this decline will now be assessed. To determine the area of the patch of fault that slipped in each earthquake, seismic moment M_O was first calculated from M_W using Eq. (50). For most of the Newdigate earthquakes, M_W values are unavailable from Hicks et al. [1], only M_L has been determined. Nonetheless, Deichmann [156] has shown that M_L can serve as an equivalent proxy of earthquake ‘size’ as M_W, provided it is appropriately calibrated. Hicks et al. [1] used the Luckett et al. [157] M_L calibration, which in the UK has superseded the more familiar Ottemöller and Sergeant [158] version. No calibration against M_W of the Luckett et al. [157] M_L formula has been reported; nonetheless, for the Newdigate events for which both M_W and M_L have been determined, both these measures are in close agreement, so M_L is used as here a proxy for M_W. Next, the radius a of the equivalent circular seismic source is determined from M_o, using Eq. (51). The source area A = π a² and the mean slip u = M_o / (μ A) are then calculated. For an ideal circular earthquake source, u = 16 (1 – ν) a δσ / (3 (1 – ν) πμ) and the maximum slip u’ is 8 (1 – ν) a δσ / ((1 – ν) πμ) (e.g., [128]), so u’ = 3 u / 2. Coseismic stress drop δσ is set as 10 MPa (a plausible upper bound), Poisson’s ratio ν is set as 0.25, and shear modulus μ ≡ B_R for ν = 0.25 (Eq. (47)); as before B_R = 50 GPa, from Bell [119]. This task was carried out for the complete Hicks et al. [1] earthquake dataset, plus the additional events listed in Table 4. The cumulative seismic moment thus obtained was ∼1.5 × 10¹⁴ N m, equivalent to a single event of M_W ∼ 3.4, with cumulative area of fault rupture ∼2.1 × 10⁵ m² and maximum coseismic slip in the largest earthquake ∼3.4 cm.

Taking 1.5 km as the length of the seismogenic fault, from Figures 1 and 7, and 70 m as the thickness of the Dinantian limestone, the area of fault in this lithology is ∼10⁵ m². Calculated on this basis, the total area of coseismic ruptures exceeds the area of the fault, and would be even greater if lower δσ were adopted. Thus, either patches of fault slipped more than once, or that the seismicity propagated into the overlying and/or underlying lithologies, although the compact hypocentre ‘cloud’ (Figure 1) indicates no clear propagation in any direction. Nonetheless, the calculations indicate that the eight largest earthquakes (with M_W ≥ 1.96) have source diameters larger than the estimated 70 m thickness of the limestone; evidently, these events either ruptured outside this layer or ruptured patches of fault that are elongated horizontally. Assuming the latter explanation, and that the overall earthquake population was distributed to produce constant overall coseismic slip across the seismogenic fault, this amount is determined as was ∼1.5 × 10¹⁴ N m / (50 GPa × 1500 m × 70 m) or ∼2.9 cm, roughly equal to the maximum slip in the largest individual earthquake (∼3.4 cm for M_W 3.25) (the value would be ∼3.3 cm, in better agreement, if the thickness of this lithology were taken as 60 m). It is thus possible that the earthquake swarm was indeed ‘self-limiting’, and that once the full extent of the seismogenic fault had slipped by this distance, the fault was effectively ‘de-stressed’ (δσ = 10 MPa having reduced ΔΦ by ∼6 MPa) and the activity died out, consistent with its observed decline and near cessation in late 2019 (Table 4). This seismogenic patch of fault is bounded to the east by the eastern end of the Newdigate Fault, but has no obvious boundary to the west, although it is noted that at its limit the downthrow of the top Portland sandstone is ≤∼60 m (Figure 7), comparable to the thickness of the Dinantian limestone; the state of stress must be different, for some reason yet to be resolved, farther west along this fault. Further analysis of this aspect is evidently warranted, given the possibility that significant seismicity might resume (two very small earthquakes having occurred in the spring of 2020; Table 4), as a result of pressure changes and fluid withdrawal arising from the planned increase in production at the Horse Hill site. A simple approach to mitigation, suggested by the present study, would be to ensure that oil production is balanced by reinjection to minimise the spatial extent of subsurface pressure changes. Suitably targeted reinjection could ensure that pressure decreases are confined to the Portland reservoir and do not propagate to the ‘beef’ and thence to the Newdigate Fault or other faults.

As already noted, the proposed physical mechanism, whereby a decrease in the fluid pressure within a fault can destabilise the fault, is the opposite of what might be termed the ‘standard’ effect, familiar from many studies of ‘fracking’: the unclamping that will occur when the fluid pressure within a fault increases and reduces the effective normal stress. The possibility of these two opposite responses to a fluid pressure change in a fault can be seen from inspection of Eq. (41). In the limit of D = 0 this adjusts to the conventional expression ΔΦ = −c ΔP_O, whereby a pressure decrease by ΔP_O would cause fault ‘clamping’ and a pressure increase would be necessary for ‘unclamping’, but if D is large enough, a pressure decrease by ΔP_O can outweigh this effect, causing fault ‘unclamping’ and coseismic slip. The ‘standard’ effect of an increase in fluid pressure causing fault unclamping is to be expected if the fault is in impermeable rocks, where the fluid pressure only acts within the fault and not within the adjoining rock volume (e.g., [11, 159]). In general, for faults within permeable rocks, one can expect these two effects to counteract each other; whether the poroelastic effect will predominate or not will depend on the conditions in each case.

In this context, it is noteworthy that much of the seismicity associated with fluid injection in the USA occurs as a result of pressure increases in faults in impermeable basement rocks, rather than in the permeable rocks into which the injection takes place. A significant case study illustrating this issue is provided by the attempt to control seismicity in the Rangely, Colorado, oilfield by varying the reservoir pressure [160]. In this instance earthquakes occurred on patches of a fault where it transects both the reservoir and deeper crustal basement, their frequency of occurrence in both locations increasing with fluid pressure. The fault unclamping effect of increasing fluid pressure was the same in the reservoir rocks as in the impermeable basement rocks. In this case, evidently, the direct effect of fluid pressure on the Mohr-Coulomb failure condition for the patch of fault in the reservoir rocks outweighed the poroelastic effect (cf. Eq. (40)), as might be expected given the relatively low permeability (∼1 mD; [160]) of the reservoir rocks, orders-of-magnitude smaller than is expected for Dinantian limestone (see above). A second case study highlighting the complexity of this issue was documented by Hornbach et al. [144] (see also, e.g., [126, 161]), at Azle near Fort Worth, Texas, where earthquake activity began in 2013 in a locality that had experienced both injection (of industrial wastewater) and production (of brine, oil and natural gas). The injection was initially suspected as the cause, on account of its very large volume, but Hornbach et al. [144] deduced that the pressure reduction caused by oil and gas production was the most important individual factor. In other instances, for example that discussed by Justinic et al. [162], authors have emphasised the proximity of hypocentres to injection points to highlight the anthropogenic cause, when many hypocentres are in fact rather deeper and indicate earthquakes within the underlying impermeable basement. Hincks et al. [163] have noted that fluid injection into faults or fractures in basement or near the contact with basement at the base of permeable sediments is statistically much more likely to result in seismicity than injection well above basement. Consideration of poroelasticity provides a natural explanation for this empirical observation.

6. Conclusions

The seismicity at Newdigate, Surrey, during 2018–2019, has been reassessed, amending aspects of the Hicks et al. [1] analysis. First-order correction of their seismic velocity model, which was too slow for the local stratigraphy, adjusts the hypocentres ∼400 m deeper than previously thought, to depths of ∼2400 m, placing them within the Palaeozoic ‘basement’ beneath the Weald Basin rather than within its Jurassic sedimentary sequence. These earthquakes involved mainly right-lateral slip on a steeply north dipping fault, part of the Newdigate fault zone (Figure 2).

Oil was produced during 2018–2019 in this vicinity from the Upper Portland Sandstone by the Brockham-X2Y and Horse Hill-1 wells. The correlation between phases of production from this reservoir and ‘bursts’ of earthquake activity (Figure 4) warrants consideration of potential geomechanical mechanisms. A conceptual model that can account for this causal connection is indicated schematically in Figure 5. It is thus suggested that the seismicity occurred within a thin (estimated ∼70 m thick) layer of permeable Dinantian limestone, hydraulically connected to the Portland reservoir via permeable strands of the Newdigate fault zone and by the highly permeable calcite ‘beef’ fabric within the Portland sandstone. It is hypothesised that past oil production at Brockham depressurized the Portland reservoir around this well and drew groundwater from the Dinantian limestone, causing it to compact and ‘unclamp’ the seismogenic fault but not sufficiently to reach the Mohr-Coulomb failure criterion to initiate coseismic slip. The resumption of production at Brockham in March 2018 caused a negative pressure pulse to propagate through the hydraulic connection to the Dinantian limestone, which, it is suggested, reached the failure threshold, initiating the first ‘burst’ of Newdigate seismicity in April 2018. Likewise, negative pressure pulse following resumption of production from the Portland reservoir at Horse Hill in February 2019 initiated a subsequent ‘burst’ of seismicity. This mechanism requires hydraulic diffusivity ∼10–20 m² s⁻¹ in the calcite ‘beef’ and ∼ 1 m² s⁻¹ in the Dinantian limestone; it predicts unclamping of fault patches by many megapascals as a result of the Horse Hill production in February 2019 and by up to ∼0.1 MPa as a result of the Brockham production in March 2018. At other times, the complexity of production patterns (e.g., from both BRX2Y and HH1 in summer 2018) and the absence of pressure data prevent any detailed conclusions being drawn, although the general correlation of seismicity with production from the Portland reservoir (Figure 4) is compelling. The proposed ‘unclamping’ effect requires consideration of the roughness of the seismogenic fault, determined by the height of its asperities and their response to compaction of the adjoining limestone. Such behaviour is particularly significant in this instance because of the high permeability of the Dinantian limestone; in impermeable rocks a reduction in pore pressure would cause fault clamping rather than unclamping. In principle this model is fully testable, but required data, notably the history of pressure variations in the oil wells, is not currently in the public domain. The recognition that this instance of seismicity is arguably caused by human activity, and the role of highly permeable hydraulic connections extending for many kilometres, has significant implications for regulation to mitigate the potential nuisance from future seismicity caused by oil production in the Weald Basin, and may also inform the understanding of anthropogenic seismicity in other settings.

Acknowledgments

Imagery and metadata for seismic lines C79-36, TWLD90-15 and TWLD90-21 were kindly provided by Malcolm Butler from the UK Onshore Geophysical Library / ‘Beneath Britain’ archive. This research did not receive any specific grant from any funding agency.

Declaration of interests

The author declares that he has no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A. Supplementary data

Provided online: https://eprints.gla.ac.uk/227098/