Seismicity induced by longwall coal mining at the Thoresby Colliery, Nottinghamshire, U.K.

Seismicity induced by longwall coal mining at the Thoresby Colliery, Nottinghamshire, U.K. Abstract The United Kingdom has a long history of deep coal mining, and numerous cases of mining-induced seismicity have been recorded over the past 50 yr. In this study, we examine seismicity induced by longwall mining at one of the United Kingdom’s last deep coal mines, the Thoresby Colliery, Nottinghamshire. After public reports of felt seismicity in late 2013 a local seismic monitoring network was installed at this site, which provided monitoring from February to October 2014. This array recorded 305 seismic events, which form the basis of our analysis. Event locations were found to closely track the position of the mining face within the Deep Soft Seam, with most events occurring up to 300 m ahead of the face position. This indicates that the seismicity is being directly induced by the mining, as opposed to being caused by activation of pre-existing tectonic features by stress transfer. However, we do not observe correlation between the rate of excavation and the rate of seismicity, and only a small portion of the overall deformation is being released as seismic energy. Event magnitudes do not follow the expected Gutenberg–Richter distribution. Instead, the observed magnitude distributions can be reproduced if a truncated power-law distribution is used to simulate the rupture areas. The best-fitting maximum rupture areas correspond to the distances between the Deep Soft Seam and the seams that over- and underlie it, which have both previously been excavated. Our inference is that the presence of a rubble-filled void (or goaf) where these seams have been removed is preventing the growth of larger rupture areas. Source mechanism analysis reveals that most events consist of dip-slip motion along near-vertical planes that strike parallel to the orientation of the mining face. These mechanisms are consistent with the expected deformation that would occur as a longwall panel advances, with the under- and overburdens moving upwards and downwards respectively to fill the void created by mining. This further reinforces our conclusion that the events are directly induced by the mining process. Similar mechanisms have been observed during longwall mining at other sites. Geomechanics, Earthquake source observations, Induced seismicity, Seismic anisotropy 1 INTRODUCTION Seismicity induced by coal mining has been a common occurrence in the United Kingdom (e.g. Redmayne 1988). Indeed, Wilson et al. (2015) estimated that between 20 and 30 per cent of all earthquakes recorded in the United Kingdom between 1970 and 2012 were induced by coal mining. From the late 1980s onwards the rate of coal production has declined significantly, as has the rate of associated earthquakes (Fig. 1). Figure 1. View largeDownload slide Deep mined coal production in the United Kingdom by year (bars) and the number of induced earthquakes per year associated with coal mining (grey line), as categorised by Wilson et al. (2015). The drop in both production and induced seismicity in 1984 is associated with the U.K. miner's strike. Figure 1. View largeDownload slide Deep mined coal production in the United Kingdom by year (bars) and the number of induced earthquakes per year associated with coal mining (grey line), as categorised by Wilson et al. (2015). The drop in both production and induced seismicity in 1984 is associated with the U.K. miner's strike. Nevertheless, seismicity associated with deep coal mining still occurs in the United Kingdom. Between December 2013 and January 2014, the United Kingdom’s national seismometer network detected a series of over 40 earthquakes near to the village of New Ollerton, Nottinghamshire. The largest of these events had a magnitude of ML = 1.7. Given the generally low levels of seismicity in the United Kingdom, the village was dubbed the ‘United Kingdom’s Earthquake Capital’ (Turvill 2014). The area has a history of seismic activity relating to coal mining (e.g. Bishop et al. 1993), and it was soon identified that the events were likely to be associated with longwall coal mining at the nearby Thoresby Colliery, which at the time was one of the few remaining deep coal mining sites in the United Kingdom. In response to the felt earthquakes, a temporary local monitoring network of surface seismometers was deployed between the 5th February and the October 30th, 2014 by the British Geological Survey (BGS). This network recorded a further 300 seismic events. The high quality of the data recorded by the local network permits a detailed study into the nature of seismicity and deformation induced by the longwall mining process. 1.1 Longwall coal mining at Thoresby The Thoresby Colliery opened in 1925. Over the history of the site, at least four different seams have been mined, including the High Hazels, Top Hard, Deep Soft and Parkgate Seams, in order from shallowest to deepest: see Edwards (1967) for a stratigraphic section showing the positions of these and other seams in the region. The Deep Soft Seam was the last to be developed, with work beginning in 2010: this was the only seam being actively mined during the study period. The colliery closed entirely in mid-2015. This was for economic reasons related to the low price of coal, not because of the induced seismicity. The Deep Soft Seam was mined using standard longwall methods: hydraulic jacks are used to support the roof while a shearing device cuts coal from the face. As the face advances, the jacks are moved forward, allowing the roof to collapse into the cavity that is left behind. The collapsed, brecciated roof material filling this void is known as goaf (e.g. Younger 2016). At Thoresby, each longwall panel has dimensions of approximately 300 m width, between 1000 and 3000 m length, and approximately 2.5 m height. 1.2 Seismicity associated with longwall coal mining Seismicity has often been associated with the longwall mining process (e.g. Cook 1976; Gibowicz et al. 1990; Bishop et al. 1993; Stec 2007; Bischoff et al. 2010; Sen et al. 2013). Seismic events associated with coal mining have often been divided into two categories: ‘mining-tectonic’ activity, produced by activation of pre-existing tectonic faults, and ‘mining-induced’ activity, directly associated with the mining excavations (e.g. Stec 2007). Observed magnitudes have typically ranged from 0.5 < ML < 3.5. At some sites event magnitudes have followed the Gutenberg & Richter (1944) distribution (e.g. Bishop et al. 1993; Kwiatek et al. 2011), while in other cases bimodal or other frequency–magnitude distributions have been observed (e.g. Stec 2007; Hudyma et al. 2008; Bischoff et al. 2010). These non-Gutenberg–Richter distributions have been attributed to the presence of characteristic length scales (the dimensions of the mined panels, for example) that provide a control on rupture dimensions and thereby event magnitudes. Analysis of event focal spheres has revealed a variety of source mechanisms in different settings (e.g. Stec 2007; Bischoff et al. 2010; Sen et al. 2013) including: non-double-couple events, indicating a volumetric component of deformation usually associated with the roof collapse process; double-couple events showing a direct relationship to mined panels, with vertical fault planes running parallel to the mining face, on which dip-slip motion occurs; and double-couple events that correspond to regional fault orientations and in situ tectonic stress conditions. In this paper we follow the processes developed in the aforementioned studies to characterize the seismicity induced by mining at the Thorseby Colliery. We begin by locating events, comparing the event locations to the propagation of the mining faces with time, and seismicity rates with the volume of coal extracted from the mine. We investigate the source characteristics of the events, using spectral analysis combined with event frequency–magnitude distributions to assess the length-scales of structures that have generated the observed events. We use shear-wave splitting analysis to image in situ stress orientations at the site, and we calculate focal mechanisms for the events to establish the orientations of fault planes and slip directions generated by the mining process. 2 EVENT DETECTION AND LOCATION 2.1 Monitoring array and event detection The local surface network deployed to monitor seismicity at the Thoresby Colliery comprised of 4 3-component Guralp 3ESP broad-band seismometers (stations NOLA, NOLD, NOLE and NOLF) and 3 vertical-component Geotech Instruments S13J short-period seismometers (NOLB, NOLC and NOLG)). The station positions are shown in Fig. 2. Events were detected using the BGS’s in-house event detection algorithm, which is based on identification of peaks in running short-time/long-time averages (STA/LTA), as described by Allen (1982). A total of 305 events were identified during the deployment of the local monitoring network. Figure 2. View largeDownload slide Map of event hypocentres, with events coloured by occurrence date. Also shown are the positions of the monitoring network (triangles) and the mining panels (brown rectangles). Panels DS-4 and DS-5 were active during the monitoring period, and the coloured bars running across these panels show the forward movement of the mining faces with time. The position of the cross-section A–B (Fig. 5) is marked by the dashed line. Figure 2. View largeDownload slide Map of event hypocentres, with events coloured by occurrence date. Also shown are the positions of the monitoring network (triangles) and the mining panels (brown rectangles). Panels DS-4 and DS-5 were active during the monitoring period, and the coloured bars running across these panels show the forward movement of the mining faces with time. The position of the cross-section A–B (Fig. 5) is marked by the dashed line. P- and S-wave arrival times were repicked manually for every event (e.g. Fig. 3). For most event-station pairs the P-wave arrival was clear and unambiguous, and so could be accurately picked (83 per cent of station-event pairs where a pick could be manually assigned). Stations NOLB, NOLC and NOLG were single, vertical component stations, so S-wave picks were not made for these stations. For smaller events with lower signal-to-noise ratios, clear S-wave arrivals were sometimes difficult to identify, resulting in a lower number of picks (74 per cent of station-event pairs where a pick could be manually assigned). Figure 3. View largeDownload slide Recorded waveforms for a larger event (ML = 1.3). The N (red), E (blue) and Z (green) components for each station are overlain. Stations NOLB, NOLC and NOLG are single (Z) component stations. The P- and S-wave picks are marked by the solid and dashed tick marks. Figure 3. View largeDownload slide Recorded waveforms for a larger event (ML = 1.3). The N (red), E (blue) and Z (green) components for each station are overlain. Stations NOLB, NOLC and NOLG are single (Z) component stations. The P- and S-wave picks are marked by the solid and dashed tick marks. The velocity model used to locate the events is taken from Bishop et al. (1993), and is listed in Table 1. The arrival time picks were inverted for the best-fitting location that minimizes the least-squares residual between modelled and picked arrival times. The search for the best-fitting location was performed using the Neighbourhood Algorithm (Sambridge 1999), and the modelled travel times were calculated using an eikonal solver (Podvin & Lecomte 1991). A map of event hypocentres is shown in Fig. 2, in which the mining panels and the position of the mining face with time are also shown. Table 1. 1D, layered, isotropic velocity model used to locate events. Model is based on that used by Bishop et al. (1993). Layer no.  Depth to layer top (m)  VP (ms−1)  VS (ms−1)  1  0  1900  1280  2  60  2750  1540  3  135  3100  1740  4  275  3500  1970  5  1019  4200  2360  6  1351  5250  2920  7  2751  6000  3370  Layer no.  Depth to layer top (m)  VP (ms−1)  VS (ms−1)  1  0  1900  1280  2  60  2750  1540  3  135  3100  1740  4  275  3500  1970  5  1019  4200  2360  6  1351  5250  2920  7  2751  6000  3370  View Large In Fig. 4 we show histograms of the event location uncertainties laterally and in depth. Note that these uncertainties pertain solely to the residuals between picked and modelled arrival times, and do not account for velocity model uncertainties. The velocity model used is based on limited site-specific data, relying mainly on regional seismic refraction surveys (Bishop et al. 1993). Figure 4. View largeDownload slide Histograms showing the lateral and depth uncertainties for the located events. Figure 4. View largeDownload slide Histograms showing the lateral and depth uncertainties for the located events. A brief sensitivity analysis suggested that velocity model uncertainties of up to 10 per cent may affect depth locations by as much as 150 m, while lateral locations are relatively unaffected. This reflects the geometry of the array, which provides reasonable azimuthal coverage but with surface stations only, such that an uncertain velocity model will primarily affect the event depths. Fig. 5 shows a cross-section of event depths relative to the coal seams. We note that, while it appears that the events are located below the seam depths, given the likely velocity model uncertainties; it is not possible to rule out that these events are actually located at the same depths as the Deep Soft Seam being mined. Figure 5. View largeDownload slide Events depths shown along cross-section A–B (see Fig. 2). The positions of the Top Hard, Deep Soft and Parkgate Seams are also marked. Note that velocity uncertainties mean that the event depths may not be particularly well constrained. Figure 5. View largeDownload slide Events depths shown along cross-section A–B (see Fig. 2). The positions of the Top Hard, Deep Soft and Parkgate Seams are also marked. Note that velocity uncertainties mean that the event depths may not be particularly well constrained. 2.2 Event locations with respect to mining activities The positions of the mining panels, and the progress of the mining face with time, have been provided by the UK Coal Authority in their Mine Abandonment Plans (2017). The position of the mining face with respect to the events can be seen in Fig. 2. It is immediately apparent that the event locations are tracking the position of the face as it moves SE along panel DS-4, before switching to DS-5 and again tracking the mining front to the SE. The monitoring period ceases when the events have propagated approximately half-way along the length of panel DS-5. We investigate the position of events in relation to the mining face in greater detail in Fig. 6, which shows a histogram of event positions relative to the mining face, along an axis parallel to the mining panels. Most events are found to occur ahead of the face, with most events occurring within 300 m of the face. This close correlation between events and the mining face implies that the events are being directly induced by mining activities, as opposed to the activation of pre-existing tectonic features, in which case we would expect the events to align along an activated fault. As per the categorisation described by Stec (2007), we characterise these as mining-induced events. Figure 6. View largeDownload slide Histograms showing the lateral position of each event relative to the mining face at the time of event occurrence, where a positive distance represents events occurring in advance of the face. Figure 6. View largeDownload slide Histograms showing the lateral position of each event relative to the mining face at the time of event occurrence, where a positive distance represents events occurring in advance of the face. However, we also note small cluster of five events that is found at greater depths (>2000 m), to the SW of the DS-4 panel. Four of these five events occurred within a single 7-hr period. Establishing the causality of these events is more difficult. It is possible that these events have been have been triggered by the static transfer of stress changes to greater depths, leading to fault activation. As per the Stec (2007) categorisation, these may be mining-tectonic events. However, it is not possible to rule out that these deeper events may in fact have a natural origin. 3 CORRELATION BETWEEN SEISMICITY AND MINING RATES? In Fig. 7 we show the volume of rock removed from the mine on a weekly basis (ΔV), the number of events per week (NE), and the cumulative seismic moment (∑MO) released per week. The volume of rock removed per week is estimated from the forward progress of the mining face, multiplied by its dimensions (width and height). To further investigate any correlation between the extracted volume and seismicity, in Fig. 8 we cross-plot these parameters. From Fig. 8 it is apparent that there is little immediate correlation between ΔV and NE and ∑MO on a weekly basis. Figure 7. View largeDownload slide Weekly rock volume extracted (black lines) compared with (a) the weekly number of recorded events and (b) the weekly cumulative seismic moment released (grey lines). Figure 7. View largeDownload slide Weekly rock volume extracted (black lines) compared with (a) the weekly number of recorded events and (b) the weekly cumulative seismic moment released (grey lines). Figure 8. View largeDownload slide Cross-plots examining potential correlation between weekly rock volume extracted and the weekly number of recorded events (a) and the weekly cumulative seismic moment released (b). In (b), the dashed lines show the expected relationship for given values of SEFF. Figure 8. View largeDownload slide Cross-plots examining potential correlation between weekly rock volume extracted and the weekly number of recorded events (a) and the weekly cumulative seismic moment released (b). In (b), the dashed lines show the expected relationship for given values of SEFF. McGarr (1976) posited a linear relationship between ΔV and ∑MO:   \begin{equation}{\rm{\Sigma }}{M_{\rm O}} \approx \mu {\rm{\Delta }}V,\end{equation} (1)where μ is the rock shear modulus. This relationship corresponds to the situation whereby all of the deformation produced by the volume change is released seismically. In reality, much of the deformation may occur aseismically. As such, Hallo et al. (2014) proposed a modification to this relationship via a ‘seismic efficiency’ term, SEFF, which describes the portion of the overall deformation that is released as seismic energy:   \begin{equation}{\rm{\Sigma }}{M_O} \approx {S_{\rm EFF}}\mu {\rm{\Delta }}V.\end{equation} (2)In some of the most well-known cases of induced seismicity, values of SEFF have been close to 1 (e.g. McGarr 2014). However, these cases represent outliers: during most industrial operations SEFF is much less that 1 (e.g. Hallo et al. 2014). The dashed lines in Fig. 8(b) show the relationship between ΔV, ∑MO and SEFF, assuming a generic value of μ = 20 GPa. We note that the observed moment release rates correspond to values of SEFF between 0.01 and 0.00001, implying that most of the deformation induced by the mining is released aseismically. This is typical for many cases of seismicity induced by a variety of industrial activities (e.g. Maxwell et al. 2008; Hallo et al. 2014). 4 EVENT MAGNITUDES AND FREQUENCY–MAGNITUDE DISTRIBUTIONS 4.1 Moment magnitude calculation Local magnitudes for the Thoresby Colliery seismicity have been computed by Butcher et al. (2017), who found that the United Kingdom’s existing local magnitude scale (Ottemöller & Sargeant 2013) is not appropriate for use when sources and receivers are within a few kilometres of each other. This is because for nearby receivers, the ray path will be predominantly through the softer, more attenuative sedimentary cover, rather than the underlying crystalline crustal rocks, as will be the case for receivers that are more distant to the event. Butcher et al. (2017) have developed an alternative local magnitude scale based on the Thoresby events, which has been recalibrated to ensure consistency between magnitude measurements made on nearby stations and those made using the UK’s permanent national monitoring network, the nearest stations of which were some distance from the Thoresby site. However, our aim here is to investigate event magnitude distributions in order to understand the length scales of structures being affected by the mining process. This therefore requires the use of moment magnitudes, since seismic moment can be directly related to rupture dimensions. We compute moment magnitudes by fitting a Brune (1970) source model to the observed S-wave displacement amplitude spectra (Fig. 9), following the method described by Stork et al. (2014). The seismic moment is determined from the amplitude of the low-frequency plateau, ΩO. Figure 9. View largeDownload slide Example displacement spectrum used to estimate moment magnitudes. The solid line shows the observed spectrum, while the dashed line shows the best-fitting Brune (1970) source model. The dot–dashed lines show the fC and ΩO values for this model. Figure 9. View largeDownload slide Example displacement spectrum used to estimate moment magnitudes. The solid line shows the observed spectrum, while the dashed line shows the best-fitting Brune (1970) source model. The dot–dashed lines show the fC and ΩO values for this model. Ideally, the measured corner frequency, fC, of the displacement spectra could be used to determine the rupture length. However, to robustly image the corner frequency, it must be significantly lower than the Nyquist frequency, fN of the recording system—Stork et al. (2014) recommend that fN > 4fC to obtain robust estimates of fC. The recording systems at Thoresby had sampling rates of 100 Hz, so fN = 50 Hz. We can use generic values for stress drop and rupture lengths to establish the expected corner frequencies for events with MW < 1. Using the relationships between rupture dimensions, seismic moment and stress drop given by Kanamori & Brodsky (2004), assuming a stress drop of 5 MPa and a rupture velocity of 2000 m s–1, the resulting corner frequency fC ≈ 30 Hz. Evidently, the fN > 4fC criteria is not expected to hold for this particular dataset. However, our observations of event magnitudes, because they are derived from the amplitude spectra at low frequencies, are robust: we therefore use these to make inferences about the length scales of the structures that have generated the observed seismic events. 4.2 Frequency–magnitude distributions The observed event magnitude distribution (EMD) is shown in Fig. 10. We show the EMDs for the overall dataset, as well as individually for the clusters associated with the DS-4 and DS-5 panels. The overall dataset is not well described by the Gutenberg & Richter (1944) distribution log 10N (M) = a − bM, where N(M) is the cumulative number of events larger than a given magnitude M, and a and b are constants to be determined. Such a distribution would be represented by a straight line in M versus log10(N) space. We note that the apparent limit on the largest event size is not an artefact of a short measuring period: while the local array was removed in October 2014, the area continues to be monitored by the BGS National Seismometer Array, which has an estimated detection capability across the United Kingdom of magnitude >2. Larger events occurring after the study period would therefore be detectable, but no such events have occurred. Figure 10. View largeDownload slide Observed frequency-magnitude distributions for the full event population (black), as well as for the DS-4 (light grey) and DS-5 (dark grey) clusters individually. Figure 10. View largeDownload slide Observed frequency-magnitude distributions for the full event population (black), as well as for the DS-4 (light grey) and DS-5 (dark grey) clusters individually. However, fault length and/or earthquake magnitude distributions that are constrained at some upper limit, leading to a fall-off from the power-law (PL) relationship at large values, have been suggested by a number of authors. At the largest scale, Richter (1958) argues that ‘a physical upper limit to the largest possible magnitude must be set by the strength of crustal rocks, in terms of the maximum strain which they are competent to support without yielding’. Similarly, Pacheco et al. (1992) argue that the rupture dimensions of very large earthquakes are limited by the thickness of the Earth's seismogenic zone (the portion of the crust that is capable of undergoing brittle failure). For continental rift zones, Scholz & Contreras (1998) suggested that the maximum length of normal faults would be limited by the flexural restoring stress and friction, and found a good match between their model and faults in the East African Rift and in Nevada. At a much smaller scale, Shapiro et al. (2013) have suggested these effects will also apply to induced seismicity, with the maximum fault size, and therefore earthquake magnitude, determined by the dimensions of the volume stimulated by human activities. To understand the observed EMDs at Thoresby, we consider the statistical distributions of fault rupture areas that might produce them. Typically, rupture areas are assumed to follow a self-similar, PL distribution (e.g. Wesnousky et al. 1983; Bonnet et al. 2001). If stress drops are assumed to be roughly constant (e.g. Abercrombie 1995) then this PL rupture area distribution will result in a PL distribution of magnitudes, that is the Gutenberg–Richter distribution. A cumulative PL distribution for rupture area will take the form:   \begin{equation}N ( L ) = \ C{A^{ - \alpha }},\end{equation} (3)where N(A) is the number of ruptures with area greater than length A, α is the PL exponent, and C is a constant. For a PL distribution, there is no upper limit to the maximum rupture area. Instead, if an upper limit to the rupture area is imposed, for example by the geometry of the mining panels, then a truncated power-law (TPL) distribution results (Burroughs & Tebbens 2001, 2002):  \begin{equation}N ( A ) = \ C\left( {{A^{ - \alpha }} - {A_{\rm MAX}}^{ - \alpha }} \right),\end{equation} (4)where AMAX is the maximum rupture area. To simulate event magnitudes based on rupture area, we use Kanamori & Brodsky (2004):   \begin{equation}{M_{\rm O}} = \ {\rm{\Delta }}\sigma {A^{3/2}},\end{equation} (5)where Δσ is the stress drop. As discussed above, the limitation of a relatively low Nyquist frequency means that we cannot measure the stress drop directly. Therefore, to estimate the PL and TPL parameters that best-fitting our observations, we initially assume a generic and arbitrary stress drop of Δσ = 5 MPa. For each of the DS-4 and DS-5 event clusters, we perform a search over the PL and TPL parameters, finding those that minimise the least-squares misfit between observed and modelled EMDs. The resulting EMDs are shown as the solid lines in Fig. 11, with the PL and TPL parameters, and the misfit for each of the models, listed in Table 2. The resulting rupture area distributions are shown in Fig. 12. Figure 11. View largeDownload slide Fitting PL (black) and TPL (grey) rupture area distributions to the DS-4 (a) and DS-5 (b) EMDs. Observed EMDs are shown by black circles. The solid lines show the best-fitting models for a fixed Δσ value, while the dashed lines show ±2 standard deviations when Δσ is varied stochastically. Figure 11. View largeDownload slide Fitting PL (black) and TPL (grey) rupture area distributions to the DS-4 (a) and DS-5 (b) EMDs. Observed EMDs are shown by black circles. The solid lines show the best-fitting models for a fixed Δσ value, while the dashed lines show ±2 standard deviations when Δσ is varied stochastically. Table 2. Best-fitting power law and truncated power-law distributions for each of the DS-4 and DS-5 clusters, and the resulting normalised misfits.   Dist. Type  α  C  AMAX  Misfit  DS-4  PL  0.47  1707  NA  5.46    TPL  0.1  743  10075  1.23  DS-5  PL  0.74  6861  NA  3.05    TPL  0.38  1536  3870  0.86    Dist. Type  α  C  AMAX  Misfit  DS-4  PL  0.47  1707  NA  5.46    TPL  0.1  743  10075  1.23  DS-5  PL  0.74  6861  NA  3.05    TPL  0.38  1536  3870  0.86  View Large Figure 12. View largeDownload slide Best-fitting rupture area PL (black) and TPL (grey) distributions for the DS-4 (a) and DS-5 (b) clusters. Figure 12. View largeDownload slide Best-fitting rupture area PL (black) and TPL (grey) distributions for the DS-4 (a) and DS-5 (b) clusters. Having established the best-fitting PL and TPL distributions with a fixed stress drop value, we then investigate the impact of a variable range of Δσ. We do this in a stochastic manner, simulating rupture area distributions based on the PL and TPL parameters, assigning stress drops randomly from a uniform distribution of 0.1 < Δσ < 20 MPa. We repeat this process over 100 iterations, and in Fig. 11 the dashed lines show the range encompassing ±2 standard deviations around the resulting mean EMD. From Fig. 11 we observe that both event populations are clearly better modelled by a TPL rupture area distribution, even when stochastic variation in Δσ is considered. Based on these results, it is worth examining whether the best-fitting values for AMAX correspond to any length-scales associated with the mining activities. There are two length scales in play that might affect rupture dimensions: the width of the mining face (approximately 300 m); and the separations between (1) the underlying Parkgate Seam, which is 35 m below the Deep Soft (Fig. 13), and (2) the overlying Top Hard Seam, which is approximately 110 m above the Deep Soft. Both seams have already been mined throughout our study area. The voids left by the longwall mining of these seams will be filled with goaf, the rubble and detritus created as the roof collapses behind the mining face. It is difficult to envisage a mechanism by which ruptures could propagate through such a rubble-filled void. Figure 13. View largeDownload slide Diagrammatic section showing the spacing between the Deep Soft Seam, and the underlying Parkgate Seam, which has already been mined out across the study area. Image taken from UK Coal Authority Mine Abandonment Plans (2017). Figure 13. View largeDownload slide Diagrammatic section showing the spacing between the Deep Soft Seam, and the underlying Parkgate Seam, which has already been mined out across the study area. Image taken from UK Coal Authority Mine Abandonment Plans (2017). Assuming circular ruptures, areas of 10 075 and 3870 m2 correspond to rupture radii of 57 and 35 m. The larger dimension radius is therefore roughly equivalent to a circular rupture extending from the Deep Soft to the Top Hard. Alternatively, assuming a rectangular rupture, the DS-4 AMAX value could correspond to a rupture with dimensions of approximately 35 × 300 m, equivalent to a rupture extending from the Deep Soft to the Parkgate, across the length of the mined face. In reality, ruptures will not be rectangular nor circular. Nevertheless, the general agreement between the dimensions of the maximum rupture area and these distances leads us to suggest that the presence of the overlying and underlying Top Hard and Parkgate seams is indeed limiting the rupture dimensions. Given the similarities between these dimensions, it is not possible to determine whether one of these features in particular is controlling the maximum rupture area. Indeed, it is likely that all three features: the width of the mining face; the distance to the underlying Parkgate Seam; and the distance to the overlying Top Hard Seam, are all playing a role in limiting the maximum rupture dimensions. 5 SEISMIC ANISOTROPY AND SHEAR-WAVE SPLITTING Shallow crustal anisotropy can be generated by several mechanisms, including: alignment of macroscopic fracture networks; the preferential alignment of microcracks due to anisotropic stress field (in practice, the microscopic and macroscopic effects often combine, as both larger-scale fracture networks and microcracks are preferentially opened or closed by the same stress field); and by the alignment of sedimentary bedding planes. Shear-wave splitting (SWS), where the velocity of a shear-wave is dependent upon the direction of travel and the polarity of the wave, is an unambiguous indicator of seismic anisotropy, and has been used previously to image stress changes induced by mining activities (Wuestefeld et al. 2011). Shear waves that propagate near-vertically will not be sensitive to horizontally layered sedimentary fabrics, which produce Vertically Transverse-Isotropy (VTI) symmetry systems. Instead, in the absence of other major structural fabrics, the fast shear wave polarisation orientation can be treated as a proxy for the direction of maximum horizontal stress (e.g. Boness & Zoback 2006). We perform SWS measurements on the Thoresby data. Accurate SWS measurements can only be obtained within the ‘S-wave window’ (Crampin & Peacock 2008), because arrivals at an incidence angle greater than ∼35° from vertical may be disturbed by S-to-P conversions at the free surface. This constraint limits the available data considerably, such that events within the S-wave window are found only on station NOLA, and for only 28 of the recorded events. We perform the SWS measurement using the automated cluster-based approach described by Teanby et al. (2004). Where larger datasets are studied, automated quality assessments such as that described by Wuestefeld et al. (2010) can be used, but in this case, given the small sample size, the quality of measurements were assessed manually. Of the 28 arrivals within the S-wave window at NOLA, 9 provided good-quality, robust results according to the diagnostic criteria specified by Teanby et al. (2004). This is a typical rate-of-return for such studies given the relatively low magnitude (and therefore signal-to-noise) of the events. An example of a robust SWS measurement is provided in Fig. 14. Figure 14. View largeDownload slide Example shear-wave splitting measurement using the method described by Teanby et al. (2004). In (a) we plot the N, E and Z components of the recorded waveforms, where P- and S-wave windows are highlighted by the shaded areas. In (b) we plot the radial and transverse components prior to and after the splitting correction, where the aim of the correction is to minimise energy on the transverse component. In (c) we plot the waveform particle motions before (solid lines) and after (dashed lines) correction. In (d) we plot the error surfaces of the correction method as a function of delay time and fast direction normalised such that the 95 per cent confidence interval (highlighted in bold) is 1. In (e) we plot the best-fitting delay times and fast directions that result from choosing different S-wave window start and end times [as indicated by the light-grey shaded zone of (a)]. Figure 14. View largeDownload slide Example shear-wave splitting measurement using the method described by Teanby et al. (2004). In (a) we plot the N, E and Z components of the recorded waveforms, where P- and S-wave windows are highlighted by the shaded areas. In (b) we plot the radial and transverse components prior to and after the splitting correction, where the aim of the correction is to minimise energy on the transverse component. In (c) we plot the waveform particle motions before (solid lines) and after (dashed lines) correction. In (d) we plot the error surfaces of the correction method as a function of delay time and fast direction normalised such that the 95 per cent confidence interval (highlighted in bold) is 1. In (e) we plot the best-fitting delay times and fast directions that result from choosing different S-wave window start and end times [as indicated by the light-grey shaded zone of (a)]. In Fig. 15(a) we show the measured fast directions in the form of an angle histogram. A dominant fast direction striking NW–SE is clearly observed. The mean fast direction azimuth is 130°. No temporal variations in SWS fast directions or percentage anisotropy were observed. The mean delay time was 43 ms, and the mean percentage S-wave anisotropy was 6.8 per cent. Figure 15. View largeDownload slide SWS and stress anisotropy. In (a) we plot an angle histogram of the measured SWS fast directions. In (b) we show regional measurements of SHmax from the World Stress Map database (Heidbach et al. 2008): ‘+’ symbols represent borehole breakouts, ‘o’ symbols represent focal mechanisms, and ‘☆’ symbols represent hydraulic fracturing data. The Thoresby site is marked by the red square. Measurements are coloured by whether they represent a thrust, normal or strike-slip stress regime (if known). Figure 15. View largeDownload slide SWS and stress anisotropy. In (a) we plot an angle histogram of the measured SWS fast directions. In (b) we show regional measurements of SHmax from the World Stress Map database (Heidbach et al. 2008): ‘+’ symbols represent borehole breakouts, ‘o’ symbols represent focal mechanisms, and ‘☆’ symbols represent hydraulic fracturing data. The Thoresby site is marked by the red square. Measurements are coloured by whether they represent a thrust, normal or strike-slip stress regime (if known). In Fig. 15(b) we compare the measured fast S-wave orientations with independent measurements for SHmax taken from the World Stress Map database (Heidbach et al. 2008). These measurements, mainly from borehole breakouts and hydraulic fracturing tests, also indicate an approximate regional SHmax strike that is to the NW–SE. We conclude that the mean measured S-wave fast polarity of 130° can be used as a proxy for SHmax at this site. 6 SOURCE MECHANISMS We compute event focal mechanisms by inverting the observed P-wave polarities and relative P wave, SH and SV wave amplitudes for the best-fitting double-couple source mechanism. In doing so, we preclude the possibility of non-double-couple sources in our inversion, as might be anticipated during mining-induced seismicity. We do this because the monitoring array consists of only 4 3-C and 3 1-C stations, which limits our ability to robustly constrain non-double-couple events. However, we note that the recovered mechanisms do a reasonable job of fitting the observed polarities, non-double-couple sources do not appear to be necessary to match the majority of our observations. Of the 305 events, a total of 65 had sufficient signal-to-noise ratios such that P-wave polarities could be robustly assigned, and produced reliable and consistent source mechanisms. These strikes, dips and rakes for these events are plotted in Fig. 16. We note three main clusters of event types, representative source mechanisms for which are also plotted. Figure 16. View largeDownload slide Source mechanisms (strike, dip and rake) for each event for which a reliable mechanism could be obtained. Three main clusters of mechanisms can be identified, representative focal spheres for which are shown. These spheres are upper-hemisphere projections where the compressive quadrants are shaded black. Figure 16. View largeDownload slide Source mechanisms (strike, dip and rake) for each event for which a reliable mechanism could be obtained. Three main clusters of mechanisms can be identified, representative focal spheres for which are shown. These spheres are upper-hemisphere projections where the compressive quadrants are shaded black. The most common source mechanism type (numbered 1 in Fig. 16) consists of events with strikes of approximately 50°, high angles of dip and rakes of between 60° and 90°. This source mechanism orientation corresponds to near-vertical planes whose strikes match the strike of the mining face, on which dip-slip movement occurs, with the side of the fault that is towards the mine moving downwards. A second, less populous source mechanism type (numbered 2 in Fig. 16) shows similar strikes and dips, but with the opposite sense of movement such that the side of the fault towards the mining face moves upwards. Similar event mechanisms—near-vertical failure planes striking parallel to the mining face with upward and downward dip-slip motion—were observed by Bischoff et al. (2010) for longwall mines in the Ruhr Area, Germany, and we share their geomechanical interpretation for these events (Fig. 17). As the coal is mined, the surrounding rock mass will collapse to fill the void. This will result in downward motion of the overlying rock (as per source mechanism type 1), and upward motion of the underlying rock (as per source mechanism type 2) along vertical planes that run parallel to the mining face. Figure 17. View largeDownload slide Geomechanical interpretation of the observed source mechanisms. As the surrounding rocks move to fill the void created by mining, dip-slip motion occurs on near-vertical slip planes oriented parallel to the mining face. Adapted from Bischoff et al. (2010). Figure 17. View largeDownload slide Geomechanical interpretation of the observed source mechanisms. As the surrounding rocks move to fill the void created by mining, dip-slip motion occurs on near-vertical slip planes oriented parallel to the mining face. Adapted from Bischoff et al. (2010). A third type of source mechanism is also observed (numbered 3 in Fig. 16), with thrust-type mechanisms occurring on steeply dipping planes that strike approximately north–south. It is possible that they result from the interaction between mining activities and pre-existing structures in the area, since the N–S orientation of these planes does not match the orientation of any feature in the mine. Using the source mechanisms for all events, we use the STRESSINVERSE iterative joint inversion algorithm described by Vavrycuk (2014) to estimate the orientations of principal stresses and the shape ratio, R (Gephart & Forsyth 1984):   \begin{equation}R\ = \frac{{{\sigma _1} - {\sigma _2}}}{{{\sigma _1} - {\sigma _3}}},\end{equation} (4)where σ1, σ2 and σ3 represent the maximum, intermediate and minimum principal stresses. The results of this inversion are listed in Table 3, and shown in Fig. 18. We note that the resulting maximum horizontal stress is sub-horizontal, with an azimuth of 144°. This is consistent, within error, with the maximum horizontal stress orientation estimated from SWS analysis. This implies that, while the orientations of the slip planes are consistent with the geometry of the mining activities, the resulting deformation is also consistent with the regional in situ stress conditions. Table 3. Principal stress orientations and Shape ratio (R) as inverted from event source mechanisms. Stress  Azimuth  Plunge (down from horizontal)  Shape ratio (R)  σ1  144°  31°  0.17  σ2  52°  2°    σ3  319°  59°    Stress  Azimuth  Plunge (down from horizontal)  Shape ratio (R)  σ1  144°  31°  0.17  σ2  52°  2°    σ3  319°  59°    View Large Figure 18. View largeDownload slide Stress tensor inversion results using the STRESSINVERSE algorithm (Vavrycuk 2014). In (a) we show a lower hemisphere projection of the P (dark grey ○) and T (light grey ◊) axes for every event, with the overall estimate for the σ1, σ2 and σ3 axes marked by a large ○, ■ and ◊, respectively. In (b) we show confidence limits for the principle stress axes, assuming ±15° error in source mechanism orientations. Figure 18. View largeDownload slide Stress tensor inversion results using the STRESSINVERSE algorithm (Vavrycuk 2014). In (a) we show a lower hemisphere projection of the P (dark grey ○) and T (light grey ◊) axes for every event, with the overall estimate for the σ1, σ2 and σ3 axes marked by a large ○, ■ and ◊, respectively. In (b) we show confidence limits for the principle stress axes, assuming ±15° error in source mechanism orientations. 7 CONCLUSIONS In this paper, we characterise the seismicity recorded during longwall mining of the Deep Soft Seam at the Thoresby Colliery, Nottinghamshire, UK. A local monitoring network was installed for 8 months, recording 305 events, with the largest event having a local magnitude of ML = 1.7. Event locations are found to track the advance of the mining faces, with most events being located up to 300 m ahead of the face. We conclude that these events are ‘mining-induced’, that is they are directly induced by the mining activity, as opposed to ‘mining-tectonic’ events, which are caused by static stress transfer producing activation of pre-existing tectonic faults. However, comparison between weekly mining rates and the rates of seismic activity do not show strong correlation. Moreover, the amount of deformation released in the form of seismic events is a small percentage of the overall deformation produced by the mining activities (in other words, most of the deformation is released aseismically). Event magnitudes do not follow the expected Gutenberg–Richter distribution. Instead, we find that the observed magnitude distribution can be reproduced by assuming that rupture areas follow a TPL distribution, whereby there is a limit to the maximum size of the rupture area. The observed maximum rupture area could correspond to several controlling features around the seam, including the width of the mining face, and the distances to the underlying Parkgate and overlying Top Hard seams, which have already been excavated. Our inference is that the presence of these rubble-filled voids where the excavated seams have been mined out creates a limit to the maximum rupture dimensions. Event source mechanism analysis shows that most events comprise dip-slip motion along near-vertical planes that strike parallel to the orientation of the mining face. This type of deformation is the expected response to the longwall mining process, and has been observed at other longwall mining sites. The observed source mechanisms are also consistent with the orientation of in situ regional stresses as inferred from SWS analysis. Acknowledgements JPV and JMK are funded by the BGS/University of Bristol Strategic Partnership in Applied Geophysics. This work was performed as part of the Bristol University Microseismicity Project (BUMPS). REFERENCES Abercrombie R.E., 1995. Earthquake source scaling relationships from −1 to 5 ML using seismograms recorded at 2.5-km depth, J. geophys. Res. , 100, 24 015–24 036. https://doi.org/10.1029/95JB02397 Allen R., 1982. Automatic phase pickers: their present use and future prospects, Bull. seism. Soc. Am. , 72, S225– S242. Bischoff M., Cete A., Fritschen R., Meier T., 2010. Coal mining induced seismicity in the Ruhr area, Germany, Pure appl. Geophys. , 167, 63– 75. https://doi.org/10.1007/s00024-009-0001-8 Google Scholar CrossRef Search ADS   Bishop I., Styles P., Allen M., 1993. Mining-induced seismicity in the Nottinghamshire coalfield, Quart. J. Eng. Geol. Hydrogeol. , 26, 253– 279. https://doi.org/10.1144/GSL.QJEGH.1993.026.004.03 Google Scholar CrossRef Search ADS   Boness N.L., Zoback M.D., 2006. Mapping stress and structurally controlled crustal shear velocity anisotropy in California, Geology , 34, 825– 828. https://doi.org/10.1130/G22309.1 Google Scholar CrossRef Search ADS   Bonnet E., Bour O., Odling N.E., Davy P., Main I., Cowie P., Berkowitz B., 2001. Scaling of fracture systems in geological media, Rev. Geophys. , 39, 347– 383. https://doi.org/10.1029/1999RG000074 Google Scholar CrossRef Search ADS   Brune J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75, 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   Burroughs S.M., 2002. The upper-truncated power law applied to earthquake cumulative frequency-magnitude distributions: evidence for a time-independent scaling parameter, Bull. seism. Soc. Am. , 92, 2983– 2993. https://doi.org/10.1785/0120010191 Google Scholar CrossRef Search ADS   Burroughs S.M., Tebbens S.F., 2001. Upper-truncated power laws in natural systems, Pure appl. geophys. , 158, 741– 757. https://doi.org/10.1007/PL00001202 Google Scholar CrossRef Search ADS   Butcher A., Luckett R., Verdon J.P., Kendall J.-M., Baptie B., Wookey J., 2017. Local magnitude discrepancies for near-event receivers: implications for the U.K. traffic-light scheme, Bull. seism. Soc. Am. , 107( 2), 532– 541. https://doi.org/10.1785/0120160225 Google Scholar CrossRef Search ADS   Cook N.G.W., 1976. Seismicity associated with mining, Eng. Geol. , 10, 99– 122. https://doi.org/10.1016/0013-7952(76)90015-6 Google Scholar CrossRef Search ADS   Crampin S., Peacock S., 2008. A review of the current understanding of seismic shear-wave splitting in the Earth's crust and common fallacies in interpretation, Wave Motion , 45, 675– 722. https://doi.org/10.1016/j.wavemoti.2008.01.003 Google Scholar CrossRef Search ADS   Edwards W.N., 1967. Geology of the country around Ollerton. Memoirs of the Geological Survey of Great Britain, Her Majesty's Stationery Office. Available at: http://pubs.bgs.ac.uk/publications.html?pubID=B01568. Gephart J.W., Forsyth D.W., 1984. An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence, J. geophys. Res. , 89, 9305– 9320. https://doi.org/10.1029/JB089iB11p09305 Google Scholar CrossRef Search ADS   Gibowicz S.J., Harjes H.-J., Schäfer M., 1990. Source parameters of seismic events at Heinrich Robert Mine, Ruhr Basin, Federal Republic of Germany: evidence for nondouble-couple events, Bull. seism. Soc. Am. , 80, 88– 109. Gutenberg B., Richter C.F., 1944. Frequency of earthquakes in California, Bull. seism. Soc. Am. , 34, 185– 188. Hallo M., Oprsal I., Eisner L., Ali M.Y., 2014. Prediction of magnitude of the largest potentially induced seismic event, J. Seismol. , 18, 421– 431. https://doi.org/10.1007/s10950-014-9417-4 Google Scholar CrossRef Search ADS   Heidbach O., Tingay M., Barth A., Reinecker J., Kurfeß D., Müller B., 2008. The World Stress Map Database Release 2008. Hudyma M., Potvin Y., Allison D., 2008. Seismic monitoring of the Northparkes Lift 2 block cave—Part 2. Production caving, J. S. Afr. Inst. Min. Metall. , 108, 421– 430. Kanamori H., Brodsky E.E., 2004. The physics of earthquakes, Rep. Prog. Phys. , 67, 1429– 1496. https://doi.org/10.1088/0034-4885/67/8/R03 Google Scholar CrossRef Search ADS   Kwiatek G., Plenkers K., Dresen G., JAGUARS Research Group, 2011. Source parameters of picoseismicity recorded at Mponeng Deep Gold Mine, South Africa: implications for scaling relations, Bull. seism. Soc. Am. , 101, 2592– 2608. https://doi.org/10.1785/0120110094 Google Scholar CrossRef Search ADS   Maxwell S.C., Shemeta J., Campbell E., Quirk D., 2008. Microseismic deformation rate monitoring, in Proceedings of the SPE Annual Technical Conference , Denver, SPE 116596. McGarr A., 1976. Seismic moments and volume changes, J. geophys. Res. , 81, 1487– 1494. https://doi.org/10.1029/JB081i008p01487 Google Scholar CrossRef Search ADS   McGarr A., 2014. Maximum magnitude earthquakes induced by fluid injection, J. geophys. Res. , 119, 1008– 1019. Google Scholar CrossRef Search ADS   Ottemoller L., Sargeant S., 2013. A local magnitude scale ML for the United Kingdom, Bull. seism. Soc. Am. , 103, 2884– 2893. https://doi.org/10.1785/0120130085 Google Scholar CrossRef Search ADS   Pacheco J.F., Scholz C.H., Sykes L.R., 1992. Changes in frequency-size relationship from small to large earthquakes, Nature , 355, 71– 73. https://doi.org/10.1038/355071a0 Google Scholar CrossRef Search ADS   Podvin P., Lecomte I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, Geophys. J. Int. , 105, 271– 284. https://doi.org/10.1111/j.1365-246X.1991.tb03461.x Google Scholar CrossRef Search ADS   Redmayne D.W., 1988. Mining induced seismicity in UK coalfields identified on the BGS national seismograph network, in Engineering Geology of Underground Movements , eds Bell F.G., Culshaw M.G., Cripps J.C., Lovell M.A., pp. 405– 413, Geological Society Engineering Geology Special Publication no. 5. Google Scholar CrossRef Search ADS   Richter C.F., 1958. Elementary Seismology , Freeman and Co. Sambridge M., 1999. Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space, Geophys. J. Int. , 138, 479– 494. https://doi.org/10.1046/j.1365-246X.1999.00876.x Google Scholar CrossRef Search ADS   Scholz C.H., Contreras J.C., 1998. Mechanics of continental rift architecture, Geology , 26, 967– 970. https://doi.org/10.1130/0091-7613(1998) Google Scholar CrossRef Search ADS   Sen A.T., Cesca S., Bischoff M., Meier T., Dahm T., 2013. Automated full moment tensor inversion of coal mining-induced seismicity, Geophys. J. Int. , 195, 1267– 1281. https://doi.org/10.1093/gji/ggt300 Google Scholar CrossRef Search ADS   Shapiro S.A., Krüger O.S., Dinske C., 2013. Probability of inducing given-magnitude earthquakes by perturbing finite volumes of rocks, J. geophys. Res. , 118, 3557– 3575. Google Scholar CrossRef Search ADS   Stec K., 2007. Characteristics of seismic activity of the Upper Silesian Coal Basin in Poland, Geophys. J. Int. , 168, 757– 768. https://doi.org/10.1111/j.1365-246X.2006.03227.x Google Scholar CrossRef Search ADS   Stork A.L., Verdon J.P., Kendall J.-M., 2014. The robustness of seismic moment and magnitudes estimated using spectral analysis, Geophys. Prospect. , 62, 862– 878. https://doi.org/10.1111/1365-2478.12134 Google Scholar CrossRef Search ADS   Teanby N.A., 2004. Automation of shear-wave splitting measurements using cluster analysis, Bull. seism. Soc. Am. , 94, 453– 463. https://doi.org/10.1785/0120030123 Google Scholar CrossRef Search ADS   Turvill W., 2014. Welcome to Britain's earthquake capital: Sleepy Nottinghamshire town has been hit by 36 tremoers in just 50 days – and geologists say mining is to blame, The Daily Mail, Available at: http://www.dailymail.co.uk/news/article-2548146/Welcome-Britains-EARTHQUAKE-capital-Sleepy-Nottinghamshire-town-hit-36-tremors-just-50-days-geologists-say-mining-blame.html, last accessed 21 January 2017. UK Coal Authority Mine Abandonment Plans, 2017. Mine Abandonment Plans are available upon application to the UK Coal Authority. Available at: https://www.gov.uk/guidance/coal-mining-records-data-deeds-and-documents. Vavrycuk V., 2014. Iterative joint inversion for stress and fault orientations from focal mechanisms, Geophys. J. Int. , 199, 69– 77. https://doi.org/10.1093/gji/ggu224 Google Scholar CrossRef Search ADS   Wesnousky S.G., Scholz C.H., Shimazaki K., Matsuda T., 1983. Earthquake frequency distribution and the mechanics of faulting, J. geophys. Res. , 88, 9331– 9340. https://doi.org/10.1029/JB088iB11p09331 Google Scholar CrossRef Search ADS   Wilson M.P., Davies R.J., Foulger G.R., Julian B.R., Styles P., Gluyas J.G., Almond S., 2015. Anthropogenic earthquakes in the UK: a national baseline prior to shale exploitation, Mar. Petrol. Geol. , 68, 1– 17. https://doi.org/10.1016/j.marpetgeo.2015.08.023 Google Scholar CrossRef Search ADS   Wuestefeld A., Al-Harrasi O., Verdon J.P., Wookey J., Kendall J.M., 2010. A strategy for automated analysis of passive microseismic data to image seismic anisotropy and fracture characteristics, Geophys. Prospect. , 58, 755– 773. https://doi.org/10.1111/j.1365-2478.2010.00891.x Google Scholar CrossRef Search ADS   Wuestefeld A., Kendall J.M., Verdon J.P., Van As A., 2011. In situ monitoring of rock fracturing using shear wave splitting analysis: an example from a mining setting, Geophys. J. Int. , 187, 848– 860. https://doi.org/10.1111/j.1365-246X.2011.05171.x Google Scholar CrossRef Search ADS   Younger P.L., 2016. How can we be sure fracking will not pollute aquifers? Lessons from a major longwall coal mining analogue (Selby, Yorkshire, UK), Earth Environ. Sci. Trans. R. Soc. Edinburgh , 106, 89– 113. https://doi.org/10.1017/S1755691016000013 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Seismicity induced by longwall coal mining at the Thoresby Colliery, Nottinghamshire, U.K.

Loading next page...
 
/lp/ou_press/seismicity-induced-by-longwall-coal-mining-at-the-thoresby-colliery-9wuzp01L7F
Publisher
Oxford University Press
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx465
Publisher site
See Article on Publisher Site

Abstract

Abstract The United Kingdom has a long history of deep coal mining, and numerous cases of mining-induced seismicity have been recorded over the past 50 yr. In this study, we examine seismicity induced by longwall mining at one of the United Kingdom’s last deep coal mines, the Thoresby Colliery, Nottinghamshire. After public reports of felt seismicity in late 2013 a local seismic monitoring network was installed at this site, which provided monitoring from February to October 2014. This array recorded 305 seismic events, which form the basis of our analysis. Event locations were found to closely track the position of the mining face within the Deep Soft Seam, with most events occurring up to 300 m ahead of the face position. This indicates that the seismicity is being directly induced by the mining, as opposed to being caused by activation of pre-existing tectonic features by stress transfer. However, we do not observe correlation between the rate of excavation and the rate of seismicity, and only a small portion of the overall deformation is being released as seismic energy. Event magnitudes do not follow the expected Gutenberg–Richter distribution. Instead, the observed magnitude distributions can be reproduced if a truncated power-law distribution is used to simulate the rupture areas. The best-fitting maximum rupture areas correspond to the distances between the Deep Soft Seam and the seams that over- and underlie it, which have both previously been excavated. Our inference is that the presence of a rubble-filled void (or goaf) where these seams have been removed is preventing the growth of larger rupture areas. Source mechanism analysis reveals that most events consist of dip-slip motion along near-vertical planes that strike parallel to the orientation of the mining face. These mechanisms are consistent with the expected deformation that would occur as a longwall panel advances, with the under- and overburdens moving upwards and downwards respectively to fill the void created by mining. This further reinforces our conclusion that the events are directly induced by the mining process. Similar mechanisms have been observed during longwall mining at other sites. Geomechanics, Earthquake source observations, Induced seismicity, Seismic anisotropy 1 INTRODUCTION Seismicity induced by coal mining has been a common occurrence in the United Kingdom (e.g. Redmayne 1988). Indeed, Wilson et al. (2015) estimated that between 20 and 30 per cent of all earthquakes recorded in the United Kingdom between 1970 and 2012 were induced by coal mining. From the late 1980s onwards the rate of coal production has declined significantly, as has the rate of associated earthquakes (Fig. 1). Figure 1. View largeDownload slide Deep mined coal production in the United Kingdom by year (bars) and the number of induced earthquakes per year associated with coal mining (grey line), as categorised by Wilson et al. (2015). The drop in both production and induced seismicity in 1984 is associated with the U.K. miner's strike. Figure 1. View largeDownload slide Deep mined coal production in the United Kingdom by year (bars) and the number of induced earthquakes per year associated with coal mining (grey line), as categorised by Wilson et al. (2015). The drop in both production and induced seismicity in 1984 is associated with the U.K. miner's strike. Nevertheless, seismicity associated with deep coal mining still occurs in the United Kingdom. Between December 2013 and January 2014, the United Kingdom’s national seismometer network detected a series of over 40 earthquakes near to the village of New Ollerton, Nottinghamshire. The largest of these events had a magnitude of ML = 1.7. Given the generally low levels of seismicity in the United Kingdom, the village was dubbed the ‘United Kingdom’s Earthquake Capital’ (Turvill 2014). The area has a history of seismic activity relating to coal mining (e.g. Bishop et al. 1993), and it was soon identified that the events were likely to be associated with longwall coal mining at the nearby Thoresby Colliery, which at the time was one of the few remaining deep coal mining sites in the United Kingdom. In response to the felt earthquakes, a temporary local monitoring network of surface seismometers was deployed between the 5th February and the October 30th, 2014 by the British Geological Survey (BGS). This network recorded a further 300 seismic events. The high quality of the data recorded by the local network permits a detailed study into the nature of seismicity and deformation induced by the longwall mining process. 1.1 Longwall coal mining at Thoresby The Thoresby Colliery opened in 1925. Over the history of the site, at least four different seams have been mined, including the High Hazels, Top Hard, Deep Soft and Parkgate Seams, in order from shallowest to deepest: see Edwards (1967) for a stratigraphic section showing the positions of these and other seams in the region. The Deep Soft Seam was the last to be developed, with work beginning in 2010: this was the only seam being actively mined during the study period. The colliery closed entirely in mid-2015. This was for economic reasons related to the low price of coal, not because of the induced seismicity. The Deep Soft Seam was mined using standard longwall methods: hydraulic jacks are used to support the roof while a shearing device cuts coal from the face. As the face advances, the jacks are moved forward, allowing the roof to collapse into the cavity that is left behind. The collapsed, brecciated roof material filling this void is known as goaf (e.g. Younger 2016). At Thoresby, each longwall panel has dimensions of approximately 300 m width, between 1000 and 3000 m length, and approximately 2.5 m height. 1.2 Seismicity associated with longwall coal mining Seismicity has often been associated with the longwall mining process (e.g. Cook 1976; Gibowicz et al. 1990; Bishop et al. 1993; Stec 2007; Bischoff et al. 2010; Sen et al. 2013). Seismic events associated with coal mining have often been divided into two categories: ‘mining-tectonic’ activity, produced by activation of pre-existing tectonic faults, and ‘mining-induced’ activity, directly associated with the mining excavations (e.g. Stec 2007). Observed magnitudes have typically ranged from 0.5 < ML < 3.5. At some sites event magnitudes have followed the Gutenberg & Richter (1944) distribution (e.g. Bishop et al. 1993; Kwiatek et al. 2011), while in other cases bimodal or other frequency–magnitude distributions have been observed (e.g. Stec 2007; Hudyma et al. 2008; Bischoff et al. 2010). These non-Gutenberg–Richter distributions have been attributed to the presence of characteristic length scales (the dimensions of the mined panels, for example) that provide a control on rupture dimensions and thereby event magnitudes. Analysis of event focal spheres has revealed a variety of source mechanisms in different settings (e.g. Stec 2007; Bischoff et al. 2010; Sen et al. 2013) including: non-double-couple events, indicating a volumetric component of deformation usually associated with the roof collapse process; double-couple events showing a direct relationship to mined panels, with vertical fault planes running parallel to the mining face, on which dip-slip motion occurs; and double-couple events that correspond to regional fault orientations and in situ tectonic stress conditions. In this paper we follow the processes developed in the aforementioned studies to characterize the seismicity induced by mining at the Thorseby Colliery. We begin by locating events, comparing the event locations to the propagation of the mining faces with time, and seismicity rates with the volume of coal extracted from the mine. We investigate the source characteristics of the events, using spectral analysis combined with event frequency–magnitude distributions to assess the length-scales of structures that have generated the observed events. We use shear-wave splitting analysis to image in situ stress orientations at the site, and we calculate focal mechanisms for the events to establish the orientations of fault planes and slip directions generated by the mining process. 2 EVENT DETECTION AND LOCATION 2.1 Monitoring array and event detection The local surface network deployed to monitor seismicity at the Thoresby Colliery comprised of 4 3-component Guralp 3ESP broad-band seismometers (stations NOLA, NOLD, NOLE and NOLF) and 3 vertical-component Geotech Instruments S13J short-period seismometers (NOLB, NOLC and NOLG)). The station positions are shown in Fig. 2. Events were detected using the BGS’s in-house event detection algorithm, which is based on identification of peaks in running short-time/long-time averages (STA/LTA), as described by Allen (1982). A total of 305 events were identified during the deployment of the local monitoring network. Figure 2. View largeDownload slide Map of event hypocentres, with events coloured by occurrence date. Also shown are the positions of the monitoring network (triangles) and the mining panels (brown rectangles). Panels DS-4 and DS-5 were active during the monitoring period, and the coloured bars running across these panels show the forward movement of the mining faces with time. The position of the cross-section A–B (Fig. 5) is marked by the dashed line. Figure 2. View largeDownload slide Map of event hypocentres, with events coloured by occurrence date. Also shown are the positions of the monitoring network (triangles) and the mining panels (brown rectangles). Panels DS-4 and DS-5 were active during the monitoring period, and the coloured bars running across these panels show the forward movement of the mining faces with time. The position of the cross-section A–B (Fig. 5) is marked by the dashed line. P- and S-wave arrival times were repicked manually for every event (e.g. Fig. 3). For most event-station pairs the P-wave arrival was clear and unambiguous, and so could be accurately picked (83 per cent of station-event pairs where a pick could be manually assigned). Stations NOLB, NOLC and NOLG were single, vertical component stations, so S-wave picks were not made for these stations. For smaller events with lower signal-to-noise ratios, clear S-wave arrivals were sometimes difficult to identify, resulting in a lower number of picks (74 per cent of station-event pairs where a pick could be manually assigned). Figure 3. View largeDownload slide Recorded waveforms for a larger event (ML = 1.3). The N (red), E (blue) and Z (green) components for each station are overlain. Stations NOLB, NOLC and NOLG are single (Z) component stations. The P- and S-wave picks are marked by the solid and dashed tick marks. Figure 3. View largeDownload slide Recorded waveforms for a larger event (ML = 1.3). The N (red), E (blue) and Z (green) components for each station are overlain. Stations NOLB, NOLC and NOLG are single (Z) component stations. The P- and S-wave picks are marked by the solid and dashed tick marks. The velocity model used to locate the events is taken from Bishop et al. (1993), and is listed in Table 1. The arrival time picks were inverted for the best-fitting location that minimizes the least-squares residual between modelled and picked arrival times. The search for the best-fitting location was performed using the Neighbourhood Algorithm (Sambridge 1999), and the modelled travel times were calculated using an eikonal solver (Podvin & Lecomte 1991). A map of event hypocentres is shown in Fig. 2, in which the mining panels and the position of the mining face with time are also shown. Table 1. 1D, layered, isotropic velocity model used to locate events. Model is based on that used by Bishop et al. (1993). Layer no.  Depth to layer top (m)  VP (ms−1)  VS (ms−1)  1  0  1900  1280  2  60  2750  1540  3  135  3100  1740  4  275  3500  1970  5  1019  4200  2360  6  1351  5250  2920  7  2751  6000  3370  Layer no.  Depth to layer top (m)  VP (ms−1)  VS (ms−1)  1  0  1900  1280  2  60  2750  1540  3  135  3100  1740  4  275  3500  1970  5  1019  4200  2360  6  1351  5250  2920  7  2751  6000  3370  View Large In Fig. 4 we show histograms of the event location uncertainties laterally and in depth. Note that these uncertainties pertain solely to the residuals between picked and modelled arrival times, and do not account for velocity model uncertainties. The velocity model used is based on limited site-specific data, relying mainly on regional seismic refraction surveys (Bishop et al. 1993). Figure 4. View largeDownload slide Histograms showing the lateral and depth uncertainties for the located events. Figure 4. View largeDownload slide Histograms showing the lateral and depth uncertainties for the located events. A brief sensitivity analysis suggested that velocity model uncertainties of up to 10 per cent may affect depth locations by as much as 150 m, while lateral locations are relatively unaffected. This reflects the geometry of the array, which provides reasonable azimuthal coverage but with surface stations only, such that an uncertain velocity model will primarily affect the event depths. Fig. 5 shows a cross-section of event depths relative to the coal seams. We note that, while it appears that the events are located below the seam depths, given the likely velocity model uncertainties; it is not possible to rule out that these events are actually located at the same depths as the Deep Soft Seam being mined. Figure 5. View largeDownload slide Events depths shown along cross-section A–B (see Fig. 2). The positions of the Top Hard, Deep Soft and Parkgate Seams are also marked. Note that velocity uncertainties mean that the event depths may not be particularly well constrained. Figure 5. View largeDownload slide Events depths shown along cross-section A–B (see Fig. 2). The positions of the Top Hard, Deep Soft and Parkgate Seams are also marked. Note that velocity uncertainties mean that the event depths may not be particularly well constrained. 2.2 Event locations with respect to mining activities The positions of the mining panels, and the progress of the mining face with time, have been provided by the UK Coal Authority in their Mine Abandonment Plans (2017). The position of the mining face with respect to the events can be seen in Fig. 2. It is immediately apparent that the event locations are tracking the position of the face as it moves SE along panel DS-4, before switching to DS-5 and again tracking the mining front to the SE. The monitoring period ceases when the events have propagated approximately half-way along the length of panel DS-5. We investigate the position of events in relation to the mining face in greater detail in Fig. 6, which shows a histogram of event positions relative to the mining face, along an axis parallel to the mining panels. Most events are found to occur ahead of the face, with most events occurring within 300 m of the face. This close correlation between events and the mining face implies that the events are being directly induced by mining activities, as opposed to the activation of pre-existing tectonic features, in which case we would expect the events to align along an activated fault. As per the categorisation described by Stec (2007), we characterise these as mining-induced events. Figure 6. View largeDownload slide Histograms showing the lateral position of each event relative to the mining face at the time of event occurrence, where a positive distance represents events occurring in advance of the face. Figure 6. View largeDownload slide Histograms showing the lateral position of each event relative to the mining face at the time of event occurrence, where a positive distance represents events occurring in advance of the face. However, we also note small cluster of five events that is found at greater depths (>2000 m), to the SW of the DS-4 panel. Four of these five events occurred within a single 7-hr period. Establishing the causality of these events is more difficult. It is possible that these events have been have been triggered by the static transfer of stress changes to greater depths, leading to fault activation. As per the Stec (2007) categorisation, these may be mining-tectonic events. However, it is not possible to rule out that these deeper events may in fact have a natural origin. 3 CORRELATION BETWEEN SEISMICITY AND MINING RATES? In Fig. 7 we show the volume of rock removed from the mine on a weekly basis (ΔV), the number of events per week (NE), and the cumulative seismic moment (∑MO) released per week. The volume of rock removed per week is estimated from the forward progress of the mining face, multiplied by its dimensions (width and height). To further investigate any correlation between the extracted volume and seismicity, in Fig. 8 we cross-plot these parameters. From Fig. 8 it is apparent that there is little immediate correlation between ΔV and NE and ∑MO on a weekly basis. Figure 7. View largeDownload slide Weekly rock volume extracted (black lines) compared with (a) the weekly number of recorded events and (b) the weekly cumulative seismic moment released (grey lines). Figure 7. View largeDownload slide Weekly rock volume extracted (black lines) compared with (a) the weekly number of recorded events and (b) the weekly cumulative seismic moment released (grey lines). Figure 8. View largeDownload slide Cross-plots examining potential correlation between weekly rock volume extracted and the weekly number of recorded events (a) and the weekly cumulative seismic moment released (b). In (b), the dashed lines show the expected relationship for given values of SEFF. Figure 8. View largeDownload slide Cross-plots examining potential correlation between weekly rock volume extracted and the weekly number of recorded events (a) and the weekly cumulative seismic moment released (b). In (b), the dashed lines show the expected relationship for given values of SEFF. McGarr (1976) posited a linear relationship between ΔV and ∑MO:   \begin{equation}{\rm{\Sigma }}{M_{\rm O}} \approx \mu {\rm{\Delta }}V,\end{equation} (1)where μ is the rock shear modulus. This relationship corresponds to the situation whereby all of the deformation produced by the volume change is released seismically. In reality, much of the deformation may occur aseismically. As such, Hallo et al. (2014) proposed a modification to this relationship via a ‘seismic efficiency’ term, SEFF, which describes the portion of the overall deformation that is released as seismic energy:   \begin{equation}{\rm{\Sigma }}{M_O} \approx {S_{\rm EFF}}\mu {\rm{\Delta }}V.\end{equation} (2)In some of the most well-known cases of induced seismicity, values of SEFF have been close to 1 (e.g. McGarr 2014). However, these cases represent outliers: during most industrial operations SEFF is much less that 1 (e.g. Hallo et al. 2014). The dashed lines in Fig. 8(b) show the relationship between ΔV, ∑MO and SEFF, assuming a generic value of μ = 20 GPa. We note that the observed moment release rates correspond to values of SEFF between 0.01 and 0.00001, implying that most of the deformation induced by the mining is released aseismically. This is typical for many cases of seismicity induced by a variety of industrial activities (e.g. Maxwell et al. 2008; Hallo et al. 2014). 4 EVENT MAGNITUDES AND FREQUENCY–MAGNITUDE DISTRIBUTIONS 4.1 Moment magnitude calculation Local magnitudes for the Thoresby Colliery seismicity have been computed by Butcher et al. (2017), who found that the United Kingdom’s existing local magnitude scale (Ottemöller & Sargeant 2013) is not appropriate for use when sources and receivers are within a few kilometres of each other. This is because for nearby receivers, the ray path will be predominantly through the softer, more attenuative sedimentary cover, rather than the underlying crystalline crustal rocks, as will be the case for receivers that are more distant to the event. Butcher et al. (2017) have developed an alternative local magnitude scale based on the Thoresby events, which has been recalibrated to ensure consistency between magnitude measurements made on nearby stations and those made using the UK’s permanent national monitoring network, the nearest stations of which were some distance from the Thoresby site. However, our aim here is to investigate event magnitude distributions in order to understand the length scales of structures being affected by the mining process. This therefore requires the use of moment magnitudes, since seismic moment can be directly related to rupture dimensions. We compute moment magnitudes by fitting a Brune (1970) source model to the observed S-wave displacement amplitude spectra (Fig. 9), following the method described by Stork et al. (2014). The seismic moment is determined from the amplitude of the low-frequency plateau, ΩO. Figure 9. View largeDownload slide Example displacement spectrum used to estimate moment magnitudes. The solid line shows the observed spectrum, while the dashed line shows the best-fitting Brune (1970) source model. The dot–dashed lines show the fC and ΩO values for this model. Figure 9. View largeDownload slide Example displacement spectrum used to estimate moment magnitudes. The solid line shows the observed spectrum, while the dashed line shows the best-fitting Brune (1970) source model. The dot–dashed lines show the fC and ΩO values for this model. Ideally, the measured corner frequency, fC, of the displacement spectra could be used to determine the rupture length. However, to robustly image the corner frequency, it must be significantly lower than the Nyquist frequency, fN of the recording system—Stork et al. (2014) recommend that fN > 4fC to obtain robust estimates of fC. The recording systems at Thoresby had sampling rates of 100 Hz, so fN = 50 Hz. We can use generic values for stress drop and rupture lengths to establish the expected corner frequencies for events with MW < 1. Using the relationships between rupture dimensions, seismic moment and stress drop given by Kanamori & Brodsky (2004), assuming a stress drop of 5 MPa and a rupture velocity of 2000 m s–1, the resulting corner frequency fC ≈ 30 Hz. Evidently, the fN > 4fC criteria is not expected to hold for this particular dataset. However, our observations of event magnitudes, because they are derived from the amplitude spectra at low frequencies, are robust: we therefore use these to make inferences about the length scales of the structures that have generated the observed seismic events. 4.2 Frequency–magnitude distributions The observed event magnitude distribution (EMD) is shown in Fig. 10. We show the EMDs for the overall dataset, as well as individually for the clusters associated with the DS-4 and DS-5 panels. The overall dataset is not well described by the Gutenberg & Richter (1944) distribution log 10N (M) = a − bM, where N(M) is the cumulative number of events larger than a given magnitude M, and a and b are constants to be determined. Such a distribution would be represented by a straight line in M versus log10(N) space. We note that the apparent limit on the largest event size is not an artefact of a short measuring period: while the local array was removed in October 2014, the area continues to be monitored by the BGS National Seismometer Array, which has an estimated detection capability across the United Kingdom of magnitude >2. Larger events occurring after the study period would therefore be detectable, but no such events have occurred. Figure 10. View largeDownload slide Observed frequency-magnitude distributions for the full event population (black), as well as for the DS-4 (light grey) and DS-5 (dark grey) clusters individually. Figure 10. View largeDownload slide Observed frequency-magnitude distributions for the full event population (black), as well as for the DS-4 (light grey) and DS-5 (dark grey) clusters individually. However, fault length and/or earthquake magnitude distributions that are constrained at some upper limit, leading to a fall-off from the power-law (PL) relationship at large values, have been suggested by a number of authors. At the largest scale, Richter (1958) argues that ‘a physical upper limit to the largest possible magnitude must be set by the strength of crustal rocks, in terms of the maximum strain which they are competent to support without yielding’. Similarly, Pacheco et al. (1992) argue that the rupture dimensions of very large earthquakes are limited by the thickness of the Earth's seismogenic zone (the portion of the crust that is capable of undergoing brittle failure). For continental rift zones, Scholz & Contreras (1998) suggested that the maximum length of normal faults would be limited by the flexural restoring stress and friction, and found a good match between their model and faults in the East African Rift and in Nevada. At a much smaller scale, Shapiro et al. (2013) have suggested these effects will also apply to induced seismicity, with the maximum fault size, and therefore earthquake magnitude, determined by the dimensions of the volume stimulated by human activities. To understand the observed EMDs at Thoresby, we consider the statistical distributions of fault rupture areas that might produce them. Typically, rupture areas are assumed to follow a self-similar, PL distribution (e.g. Wesnousky et al. 1983; Bonnet et al. 2001). If stress drops are assumed to be roughly constant (e.g. Abercrombie 1995) then this PL rupture area distribution will result in a PL distribution of magnitudes, that is the Gutenberg–Richter distribution. A cumulative PL distribution for rupture area will take the form:   \begin{equation}N ( L ) = \ C{A^{ - \alpha }},\end{equation} (3)where N(A) is the number of ruptures with area greater than length A, α is the PL exponent, and C is a constant. For a PL distribution, there is no upper limit to the maximum rupture area. Instead, if an upper limit to the rupture area is imposed, for example by the geometry of the mining panels, then a truncated power-law (TPL) distribution results (Burroughs & Tebbens 2001, 2002):  \begin{equation}N ( A ) = \ C\left( {{A^{ - \alpha }} - {A_{\rm MAX}}^{ - \alpha }} \right),\end{equation} (4)where AMAX is the maximum rupture area. To simulate event magnitudes based on rupture area, we use Kanamori & Brodsky (2004):   \begin{equation}{M_{\rm O}} = \ {\rm{\Delta }}\sigma {A^{3/2}},\end{equation} (5)where Δσ is the stress drop. As discussed above, the limitation of a relatively low Nyquist frequency means that we cannot measure the stress drop directly. Therefore, to estimate the PL and TPL parameters that best-fitting our observations, we initially assume a generic and arbitrary stress drop of Δσ = 5 MPa. For each of the DS-4 and DS-5 event clusters, we perform a search over the PL and TPL parameters, finding those that minimise the least-squares misfit between observed and modelled EMDs. The resulting EMDs are shown as the solid lines in Fig. 11, with the PL and TPL parameters, and the misfit for each of the models, listed in Table 2. The resulting rupture area distributions are shown in Fig. 12. Figure 11. View largeDownload slide Fitting PL (black) and TPL (grey) rupture area distributions to the DS-4 (a) and DS-5 (b) EMDs. Observed EMDs are shown by black circles. The solid lines show the best-fitting models for a fixed Δσ value, while the dashed lines show ±2 standard deviations when Δσ is varied stochastically. Figure 11. View largeDownload slide Fitting PL (black) and TPL (grey) rupture area distributions to the DS-4 (a) and DS-5 (b) EMDs. Observed EMDs are shown by black circles. The solid lines show the best-fitting models for a fixed Δσ value, while the dashed lines show ±2 standard deviations when Δσ is varied stochastically. Table 2. Best-fitting power law and truncated power-law distributions for each of the DS-4 and DS-5 clusters, and the resulting normalised misfits.   Dist. Type  α  C  AMAX  Misfit  DS-4  PL  0.47  1707  NA  5.46    TPL  0.1  743  10075  1.23  DS-5  PL  0.74  6861  NA  3.05    TPL  0.38  1536  3870  0.86    Dist. Type  α  C  AMAX  Misfit  DS-4  PL  0.47  1707  NA  5.46    TPL  0.1  743  10075  1.23  DS-5  PL  0.74  6861  NA  3.05    TPL  0.38  1536  3870  0.86  View Large Figure 12. View largeDownload slide Best-fitting rupture area PL (black) and TPL (grey) distributions for the DS-4 (a) and DS-5 (b) clusters. Figure 12. View largeDownload slide Best-fitting rupture area PL (black) and TPL (grey) distributions for the DS-4 (a) and DS-5 (b) clusters. Having established the best-fitting PL and TPL distributions with a fixed stress drop value, we then investigate the impact of a variable range of Δσ. We do this in a stochastic manner, simulating rupture area distributions based on the PL and TPL parameters, assigning stress drops randomly from a uniform distribution of 0.1 < Δσ < 20 MPa. We repeat this process over 100 iterations, and in Fig. 11 the dashed lines show the range encompassing ±2 standard deviations around the resulting mean EMD. From Fig. 11 we observe that both event populations are clearly better modelled by a TPL rupture area distribution, even when stochastic variation in Δσ is considered. Based on these results, it is worth examining whether the best-fitting values for AMAX correspond to any length-scales associated with the mining activities. There are two length scales in play that might affect rupture dimensions: the width of the mining face (approximately 300 m); and the separations between (1) the underlying Parkgate Seam, which is 35 m below the Deep Soft (Fig. 13), and (2) the overlying Top Hard Seam, which is approximately 110 m above the Deep Soft. Both seams have already been mined throughout our study area. The voids left by the longwall mining of these seams will be filled with goaf, the rubble and detritus created as the roof collapses behind the mining face. It is difficult to envisage a mechanism by which ruptures could propagate through such a rubble-filled void. Figure 13. View largeDownload slide Diagrammatic section showing the spacing between the Deep Soft Seam, and the underlying Parkgate Seam, which has already been mined out across the study area. Image taken from UK Coal Authority Mine Abandonment Plans (2017). Figure 13. View largeDownload slide Diagrammatic section showing the spacing between the Deep Soft Seam, and the underlying Parkgate Seam, which has already been mined out across the study area. Image taken from UK Coal Authority Mine Abandonment Plans (2017). Assuming circular ruptures, areas of 10 075 and 3870 m2 correspond to rupture radii of 57 and 35 m. The larger dimension radius is therefore roughly equivalent to a circular rupture extending from the Deep Soft to the Top Hard. Alternatively, assuming a rectangular rupture, the DS-4 AMAX value could correspond to a rupture with dimensions of approximately 35 × 300 m, equivalent to a rupture extending from the Deep Soft to the Parkgate, across the length of the mined face. In reality, ruptures will not be rectangular nor circular. Nevertheless, the general agreement between the dimensions of the maximum rupture area and these distances leads us to suggest that the presence of the overlying and underlying Top Hard and Parkgate seams is indeed limiting the rupture dimensions. Given the similarities between these dimensions, it is not possible to determine whether one of these features in particular is controlling the maximum rupture area. Indeed, it is likely that all three features: the width of the mining face; the distance to the underlying Parkgate Seam; and the distance to the overlying Top Hard Seam, are all playing a role in limiting the maximum rupture dimensions. 5 SEISMIC ANISOTROPY AND SHEAR-WAVE SPLITTING Shallow crustal anisotropy can be generated by several mechanisms, including: alignment of macroscopic fracture networks; the preferential alignment of microcracks due to anisotropic stress field (in practice, the microscopic and macroscopic effects often combine, as both larger-scale fracture networks and microcracks are preferentially opened or closed by the same stress field); and by the alignment of sedimentary bedding planes. Shear-wave splitting (SWS), where the velocity of a shear-wave is dependent upon the direction of travel and the polarity of the wave, is an unambiguous indicator of seismic anisotropy, and has been used previously to image stress changes induced by mining activities (Wuestefeld et al. 2011). Shear waves that propagate near-vertically will not be sensitive to horizontally layered sedimentary fabrics, which produce Vertically Transverse-Isotropy (VTI) symmetry systems. Instead, in the absence of other major structural fabrics, the fast shear wave polarisation orientation can be treated as a proxy for the direction of maximum horizontal stress (e.g. Boness & Zoback 2006). We perform SWS measurements on the Thoresby data. Accurate SWS measurements can only be obtained within the ‘S-wave window’ (Crampin & Peacock 2008), because arrivals at an incidence angle greater than ∼35° from vertical may be disturbed by S-to-P conversions at the free surface. This constraint limits the available data considerably, such that events within the S-wave window are found only on station NOLA, and for only 28 of the recorded events. We perform the SWS measurement using the automated cluster-based approach described by Teanby et al. (2004). Where larger datasets are studied, automated quality assessments such as that described by Wuestefeld et al. (2010) can be used, but in this case, given the small sample size, the quality of measurements were assessed manually. Of the 28 arrivals within the S-wave window at NOLA, 9 provided good-quality, robust results according to the diagnostic criteria specified by Teanby et al. (2004). This is a typical rate-of-return for such studies given the relatively low magnitude (and therefore signal-to-noise) of the events. An example of a robust SWS measurement is provided in Fig. 14. Figure 14. View largeDownload slide Example shear-wave splitting measurement using the method described by Teanby et al. (2004). In (a) we plot the N, E and Z components of the recorded waveforms, where P- and S-wave windows are highlighted by the shaded areas. In (b) we plot the radial and transverse components prior to and after the splitting correction, where the aim of the correction is to minimise energy on the transverse component. In (c) we plot the waveform particle motions before (solid lines) and after (dashed lines) correction. In (d) we plot the error surfaces of the correction method as a function of delay time and fast direction normalised such that the 95 per cent confidence interval (highlighted in bold) is 1. In (e) we plot the best-fitting delay times and fast directions that result from choosing different S-wave window start and end times [as indicated by the light-grey shaded zone of (a)]. Figure 14. View largeDownload slide Example shear-wave splitting measurement using the method described by Teanby et al. (2004). In (a) we plot the N, E and Z components of the recorded waveforms, where P- and S-wave windows are highlighted by the shaded areas. In (b) we plot the radial and transverse components prior to and after the splitting correction, where the aim of the correction is to minimise energy on the transverse component. In (c) we plot the waveform particle motions before (solid lines) and after (dashed lines) correction. In (d) we plot the error surfaces of the correction method as a function of delay time and fast direction normalised such that the 95 per cent confidence interval (highlighted in bold) is 1. In (e) we plot the best-fitting delay times and fast directions that result from choosing different S-wave window start and end times [as indicated by the light-grey shaded zone of (a)]. In Fig. 15(a) we show the measured fast directions in the form of an angle histogram. A dominant fast direction striking NW–SE is clearly observed. The mean fast direction azimuth is 130°. No temporal variations in SWS fast directions or percentage anisotropy were observed. The mean delay time was 43 ms, and the mean percentage S-wave anisotropy was 6.8 per cent. Figure 15. View largeDownload slide SWS and stress anisotropy. In (a) we plot an angle histogram of the measured SWS fast directions. In (b) we show regional measurements of SHmax from the World Stress Map database (Heidbach et al. 2008): ‘+’ symbols represent borehole breakouts, ‘o’ symbols represent focal mechanisms, and ‘☆’ symbols represent hydraulic fracturing data. The Thoresby site is marked by the red square. Measurements are coloured by whether they represent a thrust, normal or strike-slip stress regime (if known). Figure 15. View largeDownload slide SWS and stress anisotropy. In (a) we plot an angle histogram of the measured SWS fast directions. In (b) we show regional measurements of SHmax from the World Stress Map database (Heidbach et al. 2008): ‘+’ symbols represent borehole breakouts, ‘o’ symbols represent focal mechanisms, and ‘☆’ symbols represent hydraulic fracturing data. The Thoresby site is marked by the red square. Measurements are coloured by whether they represent a thrust, normal or strike-slip stress regime (if known). In Fig. 15(b) we compare the measured fast S-wave orientations with independent measurements for SHmax taken from the World Stress Map database (Heidbach et al. 2008). These measurements, mainly from borehole breakouts and hydraulic fracturing tests, also indicate an approximate regional SHmax strike that is to the NW–SE. We conclude that the mean measured S-wave fast polarity of 130° can be used as a proxy for SHmax at this site. 6 SOURCE MECHANISMS We compute event focal mechanisms by inverting the observed P-wave polarities and relative P wave, SH and SV wave amplitudes for the best-fitting double-couple source mechanism. In doing so, we preclude the possibility of non-double-couple sources in our inversion, as might be anticipated during mining-induced seismicity. We do this because the monitoring array consists of only 4 3-C and 3 1-C stations, which limits our ability to robustly constrain non-double-couple events. However, we note that the recovered mechanisms do a reasonable job of fitting the observed polarities, non-double-couple sources do not appear to be necessary to match the majority of our observations. Of the 305 events, a total of 65 had sufficient signal-to-noise ratios such that P-wave polarities could be robustly assigned, and produced reliable and consistent source mechanisms. These strikes, dips and rakes for these events are plotted in Fig. 16. We note three main clusters of event types, representative source mechanisms for which are also plotted. Figure 16. View largeDownload slide Source mechanisms (strike, dip and rake) for each event for which a reliable mechanism could be obtained. Three main clusters of mechanisms can be identified, representative focal spheres for which are shown. These spheres are upper-hemisphere projections where the compressive quadrants are shaded black. Figure 16. View largeDownload slide Source mechanisms (strike, dip and rake) for each event for which a reliable mechanism could be obtained. Three main clusters of mechanisms can be identified, representative focal spheres for which are shown. These spheres are upper-hemisphere projections where the compressive quadrants are shaded black. The most common source mechanism type (numbered 1 in Fig. 16) consists of events with strikes of approximately 50°, high angles of dip and rakes of between 60° and 90°. This source mechanism orientation corresponds to near-vertical planes whose strikes match the strike of the mining face, on which dip-slip movement occurs, with the side of the fault that is towards the mine moving downwards. A second, less populous source mechanism type (numbered 2 in Fig. 16) shows similar strikes and dips, but with the opposite sense of movement such that the side of the fault towards the mining face moves upwards. Similar event mechanisms—near-vertical failure planes striking parallel to the mining face with upward and downward dip-slip motion—were observed by Bischoff et al. (2010) for longwall mines in the Ruhr Area, Germany, and we share their geomechanical interpretation for these events (Fig. 17). As the coal is mined, the surrounding rock mass will collapse to fill the void. This will result in downward motion of the overlying rock (as per source mechanism type 1), and upward motion of the underlying rock (as per source mechanism type 2) along vertical planes that run parallel to the mining face. Figure 17. View largeDownload slide Geomechanical interpretation of the observed source mechanisms. As the surrounding rocks move to fill the void created by mining, dip-slip motion occurs on near-vertical slip planes oriented parallel to the mining face. Adapted from Bischoff et al. (2010). Figure 17. View largeDownload slide Geomechanical interpretation of the observed source mechanisms. As the surrounding rocks move to fill the void created by mining, dip-slip motion occurs on near-vertical slip planes oriented parallel to the mining face. Adapted from Bischoff et al. (2010). A third type of source mechanism is also observed (numbered 3 in Fig. 16), with thrust-type mechanisms occurring on steeply dipping planes that strike approximately north–south. It is possible that they result from the interaction between mining activities and pre-existing structures in the area, since the N–S orientation of these planes does not match the orientation of any feature in the mine. Using the source mechanisms for all events, we use the STRESSINVERSE iterative joint inversion algorithm described by Vavrycuk (2014) to estimate the orientations of principal stresses and the shape ratio, R (Gephart & Forsyth 1984):   \begin{equation}R\ = \frac{{{\sigma _1} - {\sigma _2}}}{{{\sigma _1} - {\sigma _3}}},\end{equation} (4)where σ1, σ2 and σ3 represent the maximum, intermediate and minimum principal stresses. The results of this inversion are listed in Table 3, and shown in Fig. 18. We note that the resulting maximum horizontal stress is sub-horizontal, with an azimuth of 144°. This is consistent, within error, with the maximum horizontal stress orientation estimated from SWS analysis. This implies that, while the orientations of the slip planes are consistent with the geometry of the mining activities, the resulting deformation is also consistent with the regional in situ stress conditions. Table 3. Principal stress orientations and Shape ratio (R) as inverted from event source mechanisms. Stress  Azimuth  Plunge (down from horizontal)  Shape ratio (R)  σ1  144°  31°  0.17  σ2  52°  2°    σ3  319°  59°    Stress  Azimuth  Plunge (down from horizontal)  Shape ratio (R)  σ1  144°  31°  0.17  σ2  52°  2°    σ3  319°  59°    View Large Figure 18. View largeDownload slide Stress tensor inversion results using the STRESSINVERSE algorithm (Vavrycuk 2014). In (a) we show a lower hemisphere projection of the P (dark grey ○) and T (light grey ◊) axes for every event, with the overall estimate for the σ1, σ2 and σ3 axes marked by a large ○, ■ and ◊, respectively. In (b) we show confidence limits for the principle stress axes, assuming ±15° error in source mechanism orientations. Figure 18. View largeDownload slide Stress tensor inversion results using the STRESSINVERSE algorithm (Vavrycuk 2014). In (a) we show a lower hemisphere projection of the P (dark grey ○) and T (light grey ◊) axes for every event, with the overall estimate for the σ1, σ2 and σ3 axes marked by a large ○, ■ and ◊, respectively. In (b) we show confidence limits for the principle stress axes, assuming ±15° error in source mechanism orientations. 7 CONCLUSIONS In this paper, we characterise the seismicity recorded during longwall mining of the Deep Soft Seam at the Thoresby Colliery, Nottinghamshire, UK. A local monitoring network was installed for 8 months, recording 305 events, with the largest event having a local magnitude of ML = 1.7. Event locations are found to track the advance of the mining faces, with most events being located up to 300 m ahead of the face. We conclude that these events are ‘mining-induced’, that is they are directly induced by the mining activity, as opposed to ‘mining-tectonic’ events, which are caused by static stress transfer producing activation of pre-existing tectonic faults. However, comparison between weekly mining rates and the rates of seismic activity do not show strong correlation. Moreover, the amount of deformation released in the form of seismic events is a small percentage of the overall deformation produced by the mining activities (in other words, most of the deformation is released aseismically). Event magnitudes do not follow the expected Gutenberg–Richter distribution. Instead, we find that the observed magnitude distribution can be reproduced by assuming that rupture areas follow a TPL distribution, whereby there is a limit to the maximum size of the rupture area. The observed maximum rupture area could correspond to several controlling features around the seam, including the width of the mining face, and the distances to the underlying Parkgate and overlying Top Hard seams, which have already been excavated. Our inference is that the presence of these rubble-filled voids where the excavated seams have been mined out creates a limit to the maximum rupture dimensions. Event source mechanism analysis shows that most events comprise dip-slip motion along near-vertical planes that strike parallel to the orientation of the mining face. This type of deformation is the expected response to the longwall mining process, and has been observed at other longwall mining sites. The observed source mechanisms are also consistent with the orientation of in situ regional stresses as inferred from SWS analysis. Acknowledgements JPV and JMK are funded by the BGS/University of Bristol Strategic Partnership in Applied Geophysics. This work was performed as part of the Bristol University Microseismicity Project (BUMPS). REFERENCES Abercrombie R.E., 1995. Earthquake source scaling relationships from −1 to 5 ML using seismograms recorded at 2.5-km depth, J. geophys. Res. , 100, 24 015–24 036. https://doi.org/10.1029/95JB02397 Allen R., 1982. Automatic phase pickers: their present use and future prospects, Bull. seism. Soc. Am. , 72, S225– S242. Bischoff M., Cete A., Fritschen R., Meier T., 2010. Coal mining induced seismicity in the Ruhr area, Germany, Pure appl. Geophys. , 167, 63– 75. https://doi.org/10.1007/s00024-009-0001-8 Google Scholar CrossRef Search ADS   Bishop I., Styles P., Allen M., 1993. Mining-induced seismicity in the Nottinghamshire coalfield, Quart. J. Eng. Geol. Hydrogeol. , 26, 253– 279. https://doi.org/10.1144/GSL.QJEGH.1993.026.004.03 Google Scholar CrossRef Search ADS   Boness N.L., Zoback M.D., 2006. Mapping stress and structurally controlled crustal shear velocity anisotropy in California, Geology , 34, 825– 828. https://doi.org/10.1130/G22309.1 Google Scholar CrossRef Search ADS   Bonnet E., Bour O., Odling N.E., Davy P., Main I., Cowie P., Berkowitz B., 2001. Scaling of fracture systems in geological media, Rev. Geophys. , 39, 347– 383. https://doi.org/10.1029/1999RG000074 Google Scholar CrossRef Search ADS   Brune J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75, 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   Burroughs S.M., 2002. The upper-truncated power law applied to earthquake cumulative frequency-magnitude distributions: evidence for a time-independent scaling parameter, Bull. seism. Soc. Am. , 92, 2983– 2993. https://doi.org/10.1785/0120010191 Google Scholar CrossRef Search ADS   Burroughs S.M., Tebbens S.F., 2001. Upper-truncated power laws in natural systems, Pure appl. geophys. , 158, 741– 757. https://doi.org/10.1007/PL00001202 Google Scholar CrossRef Search ADS   Butcher A., Luckett R., Verdon J.P., Kendall J.-M., Baptie B., Wookey J., 2017. Local magnitude discrepancies for near-event receivers: implications for the U.K. traffic-light scheme, Bull. seism. Soc. Am. , 107( 2), 532– 541. https://doi.org/10.1785/0120160225 Google Scholar CrossRef Search ADS   Cook N.G.W., 1976. Seismicity associated with mining, Eng. Geol. , 10, 99– 122. https://doi.org/10.1016/0013-7952(76)90015-6 Google Scholar CrossRef Search ADS   Crampin S., Peacock S., 2008. A review of the current understanding of seismic shear-wave splitting in the Earth's crust and common fallacies in interpretation, Wave Motion , 45, 675– 722. https://doi.org/10.1016/j.wavemoti.2008.01.003 Google Scholar CrossRef Search ADS   Edwards W.N., 1967. Geology of the country around Ollerton. Memoirs of the Geological Survey of Great Britain, Her Majesty's Stationery Office. Available at: http://pubs.bgs.ac.uk/publications.html?pubID=B01568. Gephart J.W., Forsyth D.W., 1984. An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence, J. geophys. Res. , 89, 9305– 9320. https://doi.org/10.1029/JB089iB11p09305 Google Scholar CrossRef Search ADS   Gibowicz S.J., Harjes H.-J., Schäfer M., 1990. Source parameters of seismic events at Heinrich Robert Mine, Ruhr Basin, Federal Republic of Germany: evidence for nondouble-couple events, Bull. seism. Soc. Am. , 80, 88– 109. Gutenberg B., Richter C.F., 1944. Frequency of earthquakes in California, Bull. seism. Soc. Am. , 34, 185– 188. Hallo M., Oprsal I., Eisner L., Ali M.Y., 2014. Prediction of magnitude of the largest potentially induced seismic event, J. Seismol. , 18, 421– 431. https://doi.org/10.1007/s10950-014-9417-4 Google Scholar CrossRef Search ADS   Heidbach O., Tingay M., Barth A., Reinecker J., Kurfeß D., Müller B., 2008. The World Stress Map Database Release 2008. Hudyma M., Potvin Y., Allison D., 2008. Seismic monitoring of the Northparkes Lift 2 block cave—Part 2. Production caving, J. S. Afr. Inst. Min. Metall. , 108, 421– 430. Kanamori H., Brodsky E.E., 2004. The physics of earthquakes, Rep. Prog. Phys. , 67, 1429– 1496. https://doi.org/10.1088/0034-4885/67/8/R03 Google Scholar CrossRef Search ADS   Kwiatek G., Plenkers K., Dresen G., JAGUARS Research Group, 2011. Source parameters of picoseismicity recorded at Mponeng Deep Gold Mine, South Africa: implications for scaling relations, Bull. seism. Soc. Am. , 101, 2592– 2608. https://doi.org/10.1785/0120110094 Google Scholar CrossRef Search ADS   Maxwell S.C., Shemeta J., Campbell E., Quirk D., 2008. Microseismic deformation rate monitoring, in Proceedings of the SPE Annual Technical Conference , Denver, SPE 116596. McGarr A., 1976. Seismic moments and volume changes, J. geophys. Res. , 81, 1487– 1494. https://doi.org/10.1029/JB081i008p01487 Google Scholar CrossRef Search ADS   McGarr A., 2014. Maximum magnitude earthquakes induced by fluid injection, J. geophys. Res. , 119, 1008– 1019. Google Scholar CrossRef Search ADS   Ottemoller L., Sargeant S., 2013. A local magnitude scale ML for the United Kingdom, Bull. seism. Soc. Am. , 103, 2884– 2893. https://doi.org/10.1785/0120130085 Google Scholar CrossRef Search ADS   Pacheco J.F., Scholz C.H., Sykes L.R., 1992. Changes in frequency-size relationship from small to large earthquakes, Nature , 355, 71– 73. https://doi.org/10.1038/355071a0 Google Scholar CrossRef Search ADS   Podvin P., Lecomte I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, Geophys. J. Int. , 105, 271– 284. https://doi.org/10.1111/j.1365-246X.1991.tb03461.x Google Scholar CrossRef Search ADS   Redmayne D.W., 1988. Mining induced seismicity in UK coalfields identified on the BGS national seismograph network, in Engineering Geology of Underground Movements , eds Bell F.G., Culshaw M.G., Cripps J.C., Lovell M.A., pp. 405– 413, Geological Society Engineering Geology Special Publication no. 5. Google Scholar CrossRef Search ADS   Richter C.F., 1958. Elementary Seismology , Freeman and Co. Sambridge M., 1999. Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space, Geophys. J. Int. , 138, 479– 494. https://doi.org/10.1046/j.1365-246X.1999.00876.x Google Scholar CrossRef Search ADS   Scholz C.H., Contreras J.C., 1998. Mechanics of continental rift architecture, Geology , 26, 967– 970. https://doi.org/10.1130/0091-7613(1998) Google Scholar CrossRef Search ADS   Sen A.T., Cesca S., Bischoff M., Meier T., Dahm T., 2013. Automated full moment tensor inversion of coal mining-induced seismicity, Geophys. J. Int. , 195, 1267– 1281. https://doi.org/10.1093/gji/ggt300 Google Scholar CrossRef Search ADS   Shapiro S.A., Krüger O.S., Dinske C., 2013. Probability of inducing given-magnitude earthquakes by perturbing finite volumes of rocks, J. geophys. Res. , 118, 3557– 3575. Google Scholar CrossRef Search ADS   Stec K., 2007. Characteristics of seismic activity of the Upper Silesian Coal Basin in Poland, Geophys. J. Int. , 168, 757– 768. https://doi.org/10.1111/j.1365-246X.2006.03227.x Google Scholar CrossRef Search ADS   Stork A.L., Verdon J.P., Kendall J.-M., 2014. The robustness of seismic moment and magnitudes estimated using spectral analysis, Geophys. Prospect. , 62, 862– 878. https://doi.org/10.1111/1365-2478.12134 Google Scholar CrossRef Search ADS   Teanby N.A., 2004. Automation of shear-wave splitting measurements using cluster analysis, Bull. seism. Soc. Am. , 94, 453– 463. https://doi.org/10.1785/0120030123 Google Scholar CrossRef Search ADS   Turvill W., 2014. Welcome to Britain's earthquake capital: Sleepy Nottinghamshire town has been hit by 36 tremoers in just 50 days – and geologists say mining is to blame, The Daily Mail, Available at: http://www.dailymail.co.uk/news/article-2548146/Welcome-Britains-EARTHQUAKE-capital-Sleepy-Nottinghamshire-town-hit-36-tremors-just-50-days-geologists-say-mining-blame.html, last accessed 21 January 2017. UK Coal Authority Mine Abandonment Plans, 2017. Mine Abandonment Plans are available upon application to the UK Coal Authority. Available at: https://www.gov.uk/guidance/coal-mining-records-data-deeds-and-documents. Vavrycuk V., 2014. Iterative joint inversion for stress and fault orientations from focal mechanisms, Geophys. J. Int. , 199, 69– 77. https://doi.org/10.1093/gji/ggu224 Google Scholar CrossRef Search ADS   Wesnousky S.G., Scholz C.H., Shimazaki K., Matsuda T., 1983. Earthquake frequency distribution and the mechanics of faulting, J. geophys. Res. , 88, 9331– 9340. https://doi.org/10.1029/JB088iB11p09331 Google Scholar CrossRef Search ADS   Wilson M.P., Davies R.J., Foulger G.R., Julian B.R., Styles P., Gluyas J.G., Almond S., 2015. Anthropogenic earthquakes in the UK: a national baseline prior to shale exploitation, Mar. Petrol. Geol. , 68, 1– 17. https://doi.org/10.1016/j.marpetgeo.2015.08.023 Google Scholar CrossRef Search ADS   Wuestefeld A., Al-Harrasi O., Verdon J.P., Wookey J., Kendall J.M., 2010. A strategy for automated analysis of passive microseismic data to image seismic anisotropy and fracture characteristics, Geophys. Prospect. , 58, 755– 773. https://doi.org/10.1111/j.1365-2478.2010.00891.x Google Scholar CrossRef Search ADS   Wuestefeld A., Kendall J.M., Verdon J.P., Van As A., 2011. In situ monitoring of rock fracturing using shear wave splitting analysis: an example from a mining setting, Geophys. J. Int. , 187, 848– 860. https://doi.org/10.1111/j.1365-246X.2011.05171.x Google Scholar CrossRef Search ADS   Younger P.L., 2016. How can we be sure fracking will not pollute aquifers? Lessons from a major longwall coal mining analogue (Selby, Yorkshire, UK), Earth Environ. Sci. Trans. R. Soc. Edinburgh , 106, 89– 113. https://doi.org/10.1017/S1755691016000013 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off