Stratigraphy of the Horse Hill 1 borehole.
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.
- anthropogenic seismicity
- calcite ‘beef’
- Weald Basin
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.
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 . 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 , 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’ (, 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., ).
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.  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.
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  and Gallois and Worssam ; Trueman , DECC , 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  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., ; 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.
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., ; Table 1). The oil–water contact for the Horse Hill reservoir has been inferred as ∼580 m TVDSS , 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.  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., ) (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.
The smaller (∼2 km × ∼1 km) Brockham reservoir (Figure 8) has a permeability of up to ∼200 mD (e.g., ), its base at 582 m TVDSS . 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.  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  fault reconstruction is favoured in the present study, being assessed as more soundly based than either the Hicks et al.  version or the earlier interpretation by Europa  (see supplement). As detailed in the supplement, during prolonged production from Brockham well BRX2Y from 1998 to 2016, amounting to 36,900 m3 , the reservoir pressure had decreased from ∼900 to ∼500 psi. Angus  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.
The geometry of the Newdigate Fault also differs radically between the Xodus  and Hicks et al.  interpretations (Figures 1 and 7). Hicks et al.  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  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  interpretation, thus also casting doubt on other aspects of the Hicks et al.  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 . 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  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., ). 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., ). Busby and Smith  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  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  showed, there is no well control for ∼20 km distance SW of the HH1 well, so no direct evidence either way. Pullan and Butler  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., ) 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. , will be considered.
2.1 Calcite ‘beef’ and its significance
Calcite ‘beef’, first reported in 1826 by Webster , consists of bedding-parallel veins of diagenetic calcite (e.g., [32, 33]). In 1835, Buckland and De la Beche  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., ). 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., ). The conditions for calcite ‘beef’ development include palaeo-temperature ∼ 70–120 °C . 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., ). The idea that the properties of calcite ‘beef’ enable the Newdigate Fault and neighbouring oil reservoirs to be hydraulically connected was suggested by Geosierra , 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., ). In the Howett  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  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., ), 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 , 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., ). 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−13 m2), being by Carey et al.  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.
The study area has been illustrated using the map (Figure 1), and seismic cross-section (Figure 2) from Hicks et al. . 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. ; it is evident that they are from the UKOGL location map. Hicks et al.  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. .
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  (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.  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. ), which merge upwards by ∼0.5 s two way time (TWT).
A significant issue to have emerged from checking the Hicks et al.  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 . The Hicks et al.  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.  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. . Furthermore, as is discussed in the supplement, Hicks et al.  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.  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|
|Lower Tunbridge Wells Sands||234.7||160.2||NR||ND|
|Purbeck Durlston Beds||396.8||322.3||NR||ND||1|
|Purbeck Main Anhydrite||604.7||530.2||0.370||2500|
|Upper Portland Sandstone||622.4||547.7||0.384||2531|
|Lower Portland Sandstone||708.4||632.5||0.451||5011|
|Kimmeridgian Micrite 1||851.3||765.4||0.532||2861|
|Kimmeridgian Micrite 2||939.7||835.5||0.581||2961|
|Main Great Oolite||1682.8||1420.1||0.950||5095|
|H (km)||VP (km s−1)||VS (km s−1)|
|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|
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 , Verdon et al. , and Hicks et al. . Baptie and Luckett  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.  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 ML < −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 ML 3.02, on 5 July) were located before any local seismograph stations were operational, using only data from regional stations. Hicks et al.  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 ), 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 ) 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.
|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|
Focal mechanisms were determined by Hicks et al.  for six events, including the largest, of ML 3.18 and moment magnitude MW 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.  and Fellgett et al.  provided syntheses of in situ stress data across much of Britain. However, these authors wrote little about the Weald Basin, Fellgett et al.  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 . 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 .
3.1 Temporal clustering
As detailed by Hicks et al. , 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 (ML 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 ML 1.28. No local seismograph stations were then in operation; Hicks et al.  estimated that the completeness threshold for earthquake detection was circa ML 2, so many smaller events were undoubtedly missed.
The second ‘burst’ (Figure 4) began at 12:28 on 27 June (ML 2.52), and included four other events above ML 2.0 (on 29 June and 5 July, and two on 18 July) including the second largest event overall (ML 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 ML 0 , 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 ML > 0 occurred at 03:21 on 18 August (ML 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 (ML 2.47), followed by two other events of ML ≥ ∼2, at 17:03 on 19 February (ML 1.98) and at 03:42 on 27 February (ML 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 ML > 0 occurred during April 2019 (on 11 and 22 April; ML 0.73 and 0.56).
The fourth ‘burst’ (Figure 4) began with a relatively large event (ML 2.35) at 00:19 on 4 May 2019. As for the preceding ‘bursts’, this seismicity thereafter began to tail off, although events with ML ∼ 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 ML 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., ). 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.  ∼4.0 m3 (∼25 barrels) of oil were produced on 23 March followed by ∼1.1, ∼0.9 and ∼1.0 m3 (∼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 , 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. . 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 ; 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. , 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.  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.  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 ), the view that faults are generally permeable, especially when critically stressed (e.g., ) 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 ), 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
where t is time, ∇2 is the Laplacian operator, and D is the hydraulic diffusivity of the medium (e.g., ). 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., ). For pressure diffusion,
(e.g., ), 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
where BR and ϕ are the bulk modulus and porosity of the rock, BF is the bulk modulus of the fluid, and BE is the ‘effective’ bulk modulus of the combination, defined as
Costain  expressed D as
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
where the subscripts B denote values of D, h and k appropriate for the layer of ‘beef’. Here, E1 is the Exponential Integral Function, defined as
E1 is not supported directly in Microsoft Excel, but using its relationship to other functions (discussed, e.g., by ) it can be evaluated indirectly. This is possible because
where Γ(ψ, x) is the Upper Incomplete Gamma Function. As Schurman  has noted, this means that an Excel formula providing a good approximation to E1(x) can be written as
where cell A1 contains x and cell B1 contains a small positive number (say, 10−8), representing ψ. Values of E1 thus calculated were checked against the tables by Harris  and Stegun and Zucker  and against the power series approximations for E1(x) for the limit of x≪1,
(e.g., ), where γ = 0.5772156649… is Euler’s Constant, and for the limit x≫1,
(e.g., ), and were found to be accurate to four significant figures or better. This method for evaluating E1(x) for all values of x is, thus, more accurate in general than the overall approximation formula proposed by Barry et al. .
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 rW, is ΔPW where
where the subscripts P denote values of D, h and k appropriate for the Portland sandstone. Given the high value of DP for the Portland sandstone, for most of the durations of the production pulses at BRX2Y and HH1, rW2/(4 DP t)≪1. One may thus approximate E1(x) using Eq. (12). Indeed, x will be so small that only the -ln(x) term need be considered. The resulting logarithmic dependence of ΔPW on t means that ΔPW 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
This equation can be differentiated with respect to t; by setting the resulting partial derivative to zero one may solve for the time delay tD for the maximum pressure variation, thus:
The maximum pressure variation ΔPM at distance r and time tD is thus
indicating that ΔPM varies inversely with r2.
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 HD, hydraulic diffusivity DD, permeability kD, and porosity ϕD. Depressurization at the point where the ‘beef’ layer intersects this fault is thus inferred to cause a reduction in groundwater pressure ΔPO 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
will require solution, where t is time and x is distance from the fault.
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  noted, this condition indicates an effective outer limit to significant pressure variations, at distance xM from the model fault, where
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 m2 s−1.
In general, poroelastic strain responses to changes in fluid pressure can be highly complex (e.g., [77, 86, 87, 88, 89]). Segall  noted that in the limit where Δσkk = 0, a reduction in fluid pressure by δP will cause a contractional strain ε = α δP / BE where α is Biot’s coefficient, defined as
Contraction will occur in the vertical z direction as well as in the x direction. Partitioning the contractional strain as Δεxx = γε and Δεzz = (1 - γ)ε,
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 ) 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
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 , p. 299 and Weisstein ,
this quantity being positive for a reduction in pressure by ΔPO.
Costain  also showed that if, rather than persisting indefinitely, the pressure change ΔPO imposed at x = 0 persists for durationδt, the resulting pressure variation δP is given by
At times t≫Δt, the pressure variation δP is given to a good approximation by
The maximum pressure variation at distance x occurs after a time delay tD given by
It follows that the maximum pressure variation at distance x and time tD is given by δPM where
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, tD 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. , based on ‘rules of thumb’ rather than derivation from first principles, predicts a value for tD 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
and using the approximate formula (Eq. (29)) as
Segall  reported that, for the Δσkk = 0 boundary condition applicable for this analysis,
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
where H is the thickness of the Dinantian limestone. If the volume of limestone thus affected has dimensions Ly parallel to and Lx perpendicular to the model fault, so ΔV = Lx × Ly × Δz or
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 Φ:
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., ). 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  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.  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 , 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 PO (Figure 11(a)), the fault being stable with Φ = -ΔΦ. The reduction in fluid pressure to PO -ΔPO 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:
or, substituting Δx from Eq. (27),
In the limit of high D t / b, at time t≫ts, where
this equation simplifies to
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  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.
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
Hydraulic conductivity K and permeability k are also interrelated thus:
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  reported a 3 m thick layer with ϕ0.25 and k 200 mD. Lee  reported BR 13 GPa as typical for dry sandstone with ϕ0.25. Taking BW 2.15 GPa for water (from ), using Eq. (4), BE for Portland Sandstone with pore space occupied by water is ∼6 GPa. At Horse Hill, Xodus  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. ), for which η = 0.9 mPa s , with ρw 1000 kg m−3 and g 9.81 m s−2. Under these conditions the oil at Brockham has η = 11 cP or 11 mPa s . Using Eq. (5), with k 200 mD, BE ∼ 6 GPa, and η 0.9 mPa s, D for Portland Sandstone is ∼1.3 m2 s−1, which can be rounded to 1 m2 s−1.
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  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 W2 / (12 η)) dP/dx, which is a solution to the more general Navier–Stokes equation for fluid flow (e.g., ). 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 ≡ W2/12. This formula gives the permeability equivalent to W = 15 μm as ∼20 D (∼2 × 10−11 m2). Overall, it is inferred that the ∼900 mD value, from Carey et al. , might be applicable to calcite ‘beef’ in the present study area. Identifying a suitable representative value for BR, 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., ). Sosa Massaro et al.  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, KB ∼ 10 GPa, and η 0.9 mPa s, D for calcite ‘beef’ can be estimated - subject to considerable uncertainty - as ∼10 m2 s−1. 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 hB 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., ), its ability to store and transmit groundwater being largely via fractures. Parameter values for Dinantian limestone include BR = 50 GPa and ϕ = 0.04, along with Young’s modulus ER = 75 GPa and Poisson’s ratio ν = 0.25, from Bell . With this set of values, BE is ∼27 GPa, and α = 1–27/50 = 0.46 (Eq. (22)). Poisson’s ratio ν ranges between 0.19 and 0.31  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  noted a range of values of K for laboratory samples of Dinantian limestone, ranging from 0.07 × 10−9 m s−1 to 0.3 × 10−9 m s−1. Lewis et al.  reported much higher values ranging from ∼10−6 m s−1 to ∼10−2 m s−1 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 m2 s−1. For comparison, Shepley  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 km2 area. His analysis reported an upper bound of 50,000 m2 day−1 or ∼ 0.6 m2 s−1. 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 m2 s−1 is appropriate for karstified Dinantian limestone, subject to considerable uncertainty. For comparison, Hornbach et al.  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 m2 s−1, in reasonable agreement. Zhang et al.  reported a nominal value of 1 m2 s−1 for this karstified Ordovician limestone.
For an ensemble of faults in different lithologies, Brodsky et al.  deduced that the typical asperity height b in length L of a fault scales as
where ζ = ∼0.6 and bO = ∼10−3 m0.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 MO was first determined using the standard formula
, where MW is moment magnitude. Next, the radius a of the equivalent circular seismic source was determined from Mo, assuming a nominal coseismic stress drop δσ:
(e.g., ). The corresponding diameter 2a is equated to fault length L to characterise asperity height. One thus obtains
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 (MW 2.5; MO 6.8 × 1012 N m) at 07:43 GMT on 14 February 2019 (located at TQ 22959 41543, after , 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−1. 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−1 (∼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 hP = 11 m, η = 0.9 mPa s, D = 1.3 m2 s−1, and kP = 35 mD (somewhat higher than the 20 mD reported by Xodus , with rW = 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, kP ∼ 430 mD would be required to match this pressure variation.
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 m3 of production from well BRX2Y on 23 March 2018, inferred to have started at 08:00 BST, caused earthquakes starting with the event (MW ∼ 2.8; MO ∼ 1.3 × 1013 N m) at 12:10 BST on 1 April 2018 (located at TQ 21992 41976, after , 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, MO ∼ 1.3 × 1013 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 m3 of production the resulting pressure decrease at r = 8 km distance is predicted to be ∼50 Pa after 9 days for DB = 10 m2 s−1. Assuming hP = 3 m, η = 0.9 mPa s, DP = 1.3 m2 s−1, and kP = 200 mD, with rW = 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 DB adjusted to 20 m2 s−1, 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.
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 DB = 10 m2 s−1 the 4 m3 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 DB increased to 20 m2 s−1, 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  that no pressure variation at HH1 was detectable in response to the pulses of production from BRX2Y.
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−5 m3 s−1 to obtain the 36,900 m3 produced (Figure 15(c)). Taking into account the dependence of DB on permeability, it is evident that for DB ∼ 10–20 m2 s−1, 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 DB = 10 m2 s−1 and in Figure 16(b) for DB = 20 m2 s−1. 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  made no mention of any pressure data for well BRX2Y that might be used to test this deduction.
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 BE = 27 GPa, α = 0.46, H = 70 m, andΔV = 36,900 m3 to balance the production, depressurization of a ∼ 6000 km2 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, xM were ∼ 6 km and ΔPO 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. . At the October 2018 workshop attention was also drawn to an InSAR-derived surface deformation map of Britain by GVL , 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 m3 / ∼200 km2 ≈ 0.2 mm), even if none of the fluid withdrawal from this limestone were recharged. The Aldiss et al.  analysis revealed vertical crustal motions at ∼1 mm a−1, 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  workshop regarding the extent to which the Newdigate earthquake ‘swarm’ might fit the standard criteria identified by Davis and Frohlich  for establishing whether instances of seismicity are anthropogenic (e.g., ). UKOG  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., ) indicates that these criteria are widely used irrespective of the geomechanical cause of any particular anthropogenic earthquake. Verdon et al.  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  with subjective numerical scores. The development of a conceptual geomechanical model for the Newdigate earthquakes supersedes the other Davis and Frohlich  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 , for earthquakes accompanying oil production near Houston, Texas. Other case studies subsequently recognised include those by Calloi et al. , Kovach , Rothé and Lui , Simpson and Leith , Pennington et al. , Wetmiller , Grasso and Wittlinger , Doser et al. , Ottemöller et al. , Dahm et al. , and Hornbach et al. , whereas works discussing physical mechanisms for such seismicity, including compaction of reservoir rocks, include those by Yerkes and Castle , Simpson et al. , Segall , Segall and Fitzgerald , Ottemöller et al. , and Soltanzadeh and Hawkes . Some of the above works (e.g., ) 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.  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., ). Compaction was initially thought to be the cause of the Groningen seismicity . 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 , 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 km2; surface subsidence ∼30 cm; cumulative production ∼1010 m3 at reservoir pressure; ), and the lower hydraulic diffusivity and elastic moduli (D 0.38 m2 s−1 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  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 MO was first calculated from MW using Eq. (50). For most of the Newdigate earthquakes, MW values are unavailable from Hicks et al. , only ML has been determined. Nonetheless, Deichmann  has shown that ML can serve as an equivalent proxy of earthquake ‘size’ as MW, provided it is appropriately calibrated. Hicks et al.  used the Luckett et al.  ML calibration, which in the UK has superseded the more familiar Ottemöller and Sergeant  version. No calibration against MW of the Luckett et al.  ML formula has been reported; nonetheless, for the Newdigate events for which both MW and ML have been determined, both these measures are in close agreement, so ML is used as here a proxy for MW. Next, the radius a of the equivalent circular seismic source is determined from Mo, using Eq. (51). The source area A = π a2 and the mean slip u = Mo / (μ 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., ), 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 μ ≡ BR for ν = 0.25 (Eq. (47)); as before BR = 50 GPa, from Bell . This task was carried out for the complete Hicks et al.  earthquake dataset, plus the additional events listed in Table 4. The cumulative seismic moment thus obtained was ∼1.5 × 1014 N m, equivalent to a single event of MW ∼ 3.4, with cumulative area of fault rupture ∼2.1 × 105 m2 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 ∼105 m2. 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 MW ≥ 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 × 1014 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 MW 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 ΔPO, whereby a pressure decrease by ΔPO would cause fault ‘clamping’ and a pressure increase would be necessary for ‘unclamping’, but if D is large enough, a pressure decrease by ΔPO 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 . 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; ) 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.  (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.  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. , 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.  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.
The seismicity at Newdigate, Surrey, during 2018–2019, has been reassessed, amending aspects of the Hicks et al.  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 m2 s−1 in the calcite ‘beef’ and ∼ 1 m2 s−1 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.
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/