The robustness of seismic moment and magnitudes estimated using spectral analysis

We present an assessment of how microseismic moment magnitude, MW , estimates vary with the method and parameters used to calculate seismic moment. This is an important topic for operators and regulators who require good magnitude estimates when monitoring induced seismicity. It is therefore imperative that these parties know and understand what errors exist in given magnitude values, something that is poorly reported. This study concentrates on spectral analysis techniques and compares MW computed in the time and frequency domains. Using recordings of MW>−1.5 events at Cotton Valley, east Texas, the maximum discrepancy between MW estimated using the different methods is 0.6 units, a significant variation. By adjusting parameters in the MW calculation we find that the radiation pattern correction term can have the most significant effect on MW , generally up to 0.8 units. Following this investigation we make a series of recommendations for estimating microseismic MW using spectral methods. Noise should be estimated and removed from recordings and an attenuation correction should be applied. The spectral level can be measured by spectral fitting or taken from the low frequency level. Significant factors in obtaining reliable microseismic MW estimates include using at least four receivers recording at ⩾1000 Hz and making radiation pattern corrections based on focal mechanism solutions, not average values.

small events (M W < 0), indicating increased fracture density around the injection point, and no larger events (M W > 0) because this could indicate that pre-existing faults are being activated, which may lead to a dissipation of injection fluid and increased seismic hazard. To assess seismic hazard an estimate of the expected magnitude of induced earthquakes can be determined by assuming the Gutenberg-Richter law for frequency-magnitude distribution (e.g., Shapiro et al. 2011;Dinske and Shapiro, 2013). This model is given by where N is the number of earthquakes with magnitude greater than M and a and b are constants to describe the equation of a straight line when plotted semi-logarithmically. The 'b-value' is the slope of the line and globally, for naturally occurring earthquakes, b ∼ 1. Variations in b-values are being used to characterize the stress regime (Schorlemmer, Wiemer and Wyss 2005) and hydraulic stimulation complexity (Wessels et al. 2011). To verify the Gutenberg-Richter law and obtain accurate b-values, it is obviously important that consistent magnitude estimates are available.
For regulatory regimes where injection activities depend on the magnitude of induced events, estimates of the size of an induced seismic event are also important. Such a 'traffic-light' system has been adopted for gas injection in the Netherlands (Dost et al. 2012) and has been proposed to regulate hydraulic fracturing activities in the UK (Green, Styles and Baptie 2012;Davey 2012). If magnitude-dependent traffic-light systems are to be imposed, then it is important that regulators understand why uncertainties in reported magnitudes exist and how large these uncertainties might be.
Various methods exist to estimate magnitude. Empirical methods have long been employed, for example in local magnitude scales (e.g., Richter 1935). However, M W estimates, derived from M 0 , are required to study the physical properties of earthquakes, such as source parameters and scaling relations, because these properties cannot be deduced from empirical scales. Unfortunately earthquake magnitudes are often reported without any indication of the type of magnitude or how the magnitude is calculated. This lack of information makes it difficult to compare datasets or verify results. Seismic event magnitudes are routinely estimated for earthquakes induced by industrial activities such as hydraulic fracturing, geothermal exploration, waste water injection or CO 2 sequestration, but these values are usually given without any processing information or any indication of errors in M 0 estimates.
The aim of this paper is to assess the importance and influence of different parameters on estimated M W . Understanding the effect of using different methods to calculate magnitudes will aid operators to assess the usefulness of stimulated volume estimates and allow regulators to understand what errors there might be in reported magnitude estimates when they come to enforce traffic-light systems for operations. There has been little previous discussion in the hydraulic fracturing community of the errors in estimated microseismic magnitudes and how this might affect operations. Shemeta and Anderson (2010) report large variations in magnitudes calculated by different processors for one dataset and discuss why this might be. Studies quantifying the effect of different approaches to model earthquake spectra and to estimate M 0 also report significant differences for individual earthquakes (Ide et al. 2003;Oye, Hilmar and Roth 2005;Sonley and Aber-crombie 2006), although these studies concentrate mostly on attenuation, Q, effects on spectra and the resulting impact on source parameters such as stress drop and apparent stress. In the present study we assess variations in M W estimated using spectral analysis because this is the method currently routinely used in microseismic studies of induced seismicity. Richter (1935) developed the original local magnitude scale for earthquakes in California and recorded by a Wood-Anderson seismograph. Further magnitude scales have been developed and proposed for different regions, instruments and source-receiver distances. These scales use maximum recorded amplitudes, A, of seismic waves as a measure of earthquake size, and their general form is given by

M A G N I T U D E C A L C U L A T I O N S
where T is the dominant period; F is a correction for the variation in amplitude with earthquake depth, h, and distance, , from the recording instrument; and C is a scale factor dependent on the region of interest. Since these magnitude scales are empirical they do not allow the derivation of any physical earthquake properties. To resolve this issue, the moment magnitude scale was proposed by Hanks and Kanamori (1979): with M 0 in Nm. Estimating M 0 requires more complicated analysis than measuring A, which is why local magnitude scales are still often used and empirical relationships have been developed to relate the two scales (e.g., Pearson 1982;Edwards et al. 2010). However, to investigate the physical properties (e.g., source parameters) of induced and natural seismicity, M 0 must be estimated and therefore M W used.

S E I S M I C M O M E N T C A L C U L A T I O N S
M 0 is the best measure of earthquake size and is defined as where μ is the shear modulus andD is the average slip on the fault with area S. M 0 can be calculated through spectral analysis or waveform inversion. For large earthquakes, M 0 is routinely calculated through waveform inversion by the Global Centroid-Moment-Tensor (CMT) Project (www.globalcmt.org; Dziewonski, Chou and Woodhouse 1981;Ekström, Nettles and Dziewonski 2012). For microseismic events, spectral analysis methods are more often used to calculate M 0 (e.g., Urbancic et al. 1996;Kwiatek et al. 2011), although full waveform inversions for moment tensor solutions have recently been, and are currently, under development (Song and Toksöz 2011;O'Toole et al. 2013). M 0 can be estimated from seismograms using recordings in the time domain or frequency domain and we investigate both these methods to provide general recommendations on M W calculation and to illustrate the relative effect that different parameters have on magnitude estimates. This study does not consider the conversion of empirical magnitude scales to M W because these scales are region dependent and the conversions for small and microseismic magnitudes are discussed by several authors (e.g., Pearson 1982;Jost et al. 1998;Clinton, Hauksson and Solanki 2006;Edwards et al. 2010;Bethmann, Deichmann and Mai 2011).

Time domain calculations
As given by Aki and Richards (2002), the far-field displacements from a point source for P-and S-waves in a homogeneous whole space are respectively, where R P,S are the P-or S-wave radiation pattern correction terms, ρ is the rock density, v is the P-or S-wave velocity, r is the source-receiver distance andṀ 0 t − r v is the moment rate function. Figure 1 (b) illustrates the model moment rate function in the time domain. Body-wave radiation patterns are non-uniform, but if no focal mechanism solution is available then the average radiation pattern correction, R, 0.44 for P-waves and 0.60 for S-waves (Boore and Boatwright 1984), can be used. The seismic moment can then be estimated from the average P-and S-wave moments given by This calculation approximates geometrical spreading as 1/r (valid for the short distances in the calculations below) and assumes recordings are corrected for instrument response. If the receiver is at the surface, a free-surface correction must also be made (Aki and Richards 2002). Thus, to estimate M 0 in the time domain, the seismogram must be converted to displacement and integrated over 0 → ∞. In practice the integral is over the length of the displacement pulse, between times t 1 and t 2 (as illustrated in Fig.1 (a,b)). The main factors affecting this calculation are the accuracy of the arrival time pick (t 1 ), the estimated pulse width (t 2 − t 1 ), the signal-to-noise ratio (SNR) and any artefacts resulting from the conversion to a displacement seismogram. Due to scattering, the energy does not arrive as a single pulse on the displacement seismogram, and t 2 is chosen to include this scattered energy, which results in a window usually about 2 cycles long. It should also be noted that the estimate of M 0 in equation (6) is uncorrected for intrinsic attenuation (Q).

Frequency domain calculations
A Fourier transform allows M 0 to be estimated in the frequency domain (see Fig. 1(c)) and the equivalent of equation (6) is where 0 is the low frequency level of the amplitude spectrum and is equal to the area under the displacement pulse. This can be measured directly from the amplitude spectra assuming low frequencies are correctly recorded. Alternatively, the spectral level can be estimated by fitting a source model spectrum to amplitude spectra. The most frequently applied source models are: 1. the Brune (1970) model, where t is the travel-time, f is the frequency and f c is the corner frequency; 2. the Boatwright (1980) model of the form, A fitting algorithm to solve for 0 , f c and Q simultaneously allows further investigations into source parameters and scaling relationships and provides a measure of intrinsic attenuation.

C A S E S T U D Y -C O T T O N V A L L E Y
To conduct this experiment to vary input parameters and observe their effect on magnitude estimates, we use microseismic data recorded at Cotton Valley, east Texas, USA. This dataset was chosen because it includes 30 microseismic events located The seismic moment rate function, M 0 (t), is proportional to the far-field displacement pulse (Aki and Richards 2002). (c) In the frequency ( f ) domain, the moment rate spectrum is the Fourier transform ofṀ 0 (t). At low frequencies the spectral amplitude becomes constant and is proportional to seismic moment. f c is the corner frequency. FFT is the fast Fourier transform.
less than 1 km from the receivers and with reported magnitudes > −1.5 (Rutledge, Phillips and Mayerhofer 2004) and has some recordings with good SNR. In addition, two borehole arrays were in operation and all 30 events used were recorded by at least 11 receivers, thereby allowing a test of how the array configuration might affect magnitude estimates. Since M W for these events were reported by Rutledge et al. (2004) we can compare independent estimates of magnitude. Rutledge et al. (2004) applied Andrews' (1986) spectral analysis technique to estimate M 0 and reported mean M W values obtained from P-and S-wave recordings from several stations. In this method a Brune spectrum is assumed and 0 is given by where S D is an integral of the displacement power spectrum and S V is an integral of the velocity power spectrum. Rutledge et al. (2004) also calculated focal mechanism solutions for the events studied here. This allows us to demonstrate the effect on M W of assuming an average radiation pat-tern for all recordings compared with using the corrections given by Aki and Richards (2002). It should also be noted that Rutledge et al. compute double-couple solutions only and any non-double-couple components would result in errors in the radiation pattern coefficients used here. However, for the purposes of assessing the dependence of M W estimates on radiation pattern corrections, the double-couple solutions provide suitable data.
A series of six hydraulic fracture tests were conducted in the Carthage Cotton Valley gas field in east Texas in 1997. This was done to test fracture imaging methods, including microseismic methods, and microseismic data were collected and processed in two wells (Mayerhofer et al. 2000). Figure 2 shows the location of the 20 geophones in Arrays 1 and 2 used in this study and the related treatment interval. The locations of the 30 events studied here with M W > −1.5 (as given by Rutledge et al. 2004) are also shown in Fig. 2. Array 2 spans the treatment interval and the event depths and therefore, combined with a second array, provides good coverage to record the seismicity. The data are recorded at 1000 Hz by 30 Hz geophones (Walker 1997).

D A T A T R E A T M E N T
In the processing of the seismograms we remove the mean; rotate the seismograms into the P, S H and SV directions; and convert to displacement seismograms. An appropriate time window is then chosen to perform M 0 calculations. This depends on the source-receiver distance and earthquake magnitude because larger magnitude events have longer source durations, and a longer coda develops the further the waves travel. The start of the time window, t 1 , is the arrival time for the P-or S-waves, and the end of the time window, t 2 , is at 2+ cycles of the P-or S-wave seismograms, while the signal is visible above the noise level. This is to maximize the event energy included in the calculation while minimizing the contribution of other sources of energy. Figure  3(a) shows an example displacement seismogram to illustrate the time window selection for the M 0 calculations, given by equation (6).
For the frequency domain calculations, after a Fourier transform zero-padded to a length of 2048, the amplitude spectra are smoothed using a five-point window, and the P-wave pre-signal noise is subtracted from the spectra; an example is shown in Fig. 3(b). We then estimate M 0 using equation (7) with 0 estimated by four different methods. Two of these estimates derive from least-squares best-fit Brune (1970) and Boatwright (1980) models, as given by equations (8) and (9), respectively. The difference in shape between the two models is illustrated in Fig. 3(b). The methodology used is similar to that described in Abercrombie (1995) and Prejean and Ellsworth (2001). First a grid search was conducted to find approximate values for 0 , f c and Q. Subsequently, a Levenberg-Marquardt algorithm was applied to find the best-fitting model.
The final two estimates for 0 were obtained from direct estimation of the low frequency level of the amplitude spectrum. For these estimates an average of the amplitude spectrum between 40 Hz and 50 Hz was taken with spectra corrected and uncorrected for Q. The Q values were taken from the results of the Brune model fit. The spectral level was not taken at lower frequencies and we do not perform least-squares fits below 30 Hz because we do not correct for instrument response below 30 Hz. The instruments have a flat frequency response above 30 Hz and we correct for this in the processing. The amplitude spectra are observed to be flat up to at least 60 Hz, and therefore taking 0 between 40 and 50 Hz gives similar results as would be achieved by using lower frequencies.
These treatments lead to five estimates for M 0 and, hence, M W : In addition we compare these M W values to those estimated by Rutledge et al. (2004), M R W . Some parameters affect the above calculations independent of seismogram processing in the frequency or time domain, and in this study we consider how event location accuracy (affecting r ), the velocity model accuracy, the radiation pattern correction and the rock density affect microseismic M W estimates. We also investigate the effect the following factors have on the calculations: r the length of the time window (t 2 − t 1 ); r the recording frequency; r using P-and/or S-waves; r SNR; r the number of recording stations.
Other factors affecting magnitude calculations include the accuracy of instrument calibration, how well receivers are coupled to the borehole and scattering effects. The calibration of instruments and their set-up are outside the scope of this study but it is an important consideration that should be understood before installing microseismic monitoring arrays (Baig et al. 2013). Scattering results in the development of coda waves, which can be used to estimate source parameters through spectral ratio techniques (e.g., Hough 1997;Ide et al. 2003;Oye et al. 2005). Spectral ratios are not considered in this study because these methods require co-located events with similar mechanisms and therefore are not always applicable.
In the results presented below, M W is estimated for the available recordings that satisfy the stated conditions. Each of the 30 events is recorded by at least 11 receivers, providing at least 330 potential results. However, depending on the conditions imposed (e.g. , SNR > 3.0), only up to 98 recordings meet each set of conditions and the figures below illustrate up to 98 results.

Method dependence
We consider how M W estimates (averaged for P-and S-waves) vary depending on the method chosen to compute M 0 for the Cotton Valley data. The estimates made in this study use recordings with SNR > 3.0 between 30 Hz and 300 Hz and time windows individually picked for each 1000 Hz recording, with a radiation pattern correction applied based on the focal mechanism solutions reported by Rutledge et al. (2004). These conditions form our baseline against which we perform tests in the following results sections. Figure 4(a-c) illustrate how M Q W estimated for individual receivers compares to a) M t W , b) M Br W and c) M W . The greatest discrepancy is between the time and frequency domains (Fig. 4a), up to 0.3 units, mostly due to the fact that the results in the frequency domain are corrected for Q. As expected, M W estimated using Q-corrected spectra are larger than the values estimated from uncorrected spectra, with a maximum difference of 0.25 units and a mean difference of 0.13 units (Fig. 4c). In general, the different methods used in this study produce M W differences <0.3 magnitude units. The effect of spectral model assumptions on M W values is illustrated in Fig. 4(e), and the difference between M Br W and M Bo W is <0.1 units. This is expected because the 0 estimate should be independent of the fall off in amplitude at higher frequencies.
The maximum difference between any of the values given in Fig. 4 Fig. 4 are used to make comparisons, and results are mostly presented as deviations from these baseline M W , M W . It is stated in the results below where estimates differing from this so-called baseline are presented.

Radiation pattern correction
M W values were estimated using R calculated from the focal mechanism solutions given in Rutledge et al. (2004) and using the average values for the radiation pattern (0.44 for P-waves and 0.60 for S-waves). The difference between these magnitude estimates, M W , is illustrated in Fig. 5 for P-and S-phases separately. For the particular event mechanisms and station configuration in this study, R deviates significantly from the average value and this results in large variations in estimated magnitudes, up to 1.1 magnitude units. It appears in Fig. 5 that M W increase with magnitude for S-waves. This effect occurs because the larger events are reported by a greater number of more distributed sensors. Some of these sensors are closer to a node in the S-wave radiation pattern than the P-waves radiation pattern. Therefore, using R calculated from the Rutledge et al. (2004) focal mechanism solutions affects S-phase M W estimates to a greater extent than P-phase estimates. The single result with M W > 0.8 has an SV correction <0.01. This highlights that the radiation pattern correction is particularly important if the receiver is situated near a node in the radiation pattern because R changes very rapidly with a small change in fault plane solution or event location. Accurate and precise event location and fault plane solutions are required to avoid large errors in estimated M W . Note, in reporting the other results presented in this paper, the minimum correction is set to 0.01 because if the correction were truly smaller than this it is unlikely that the arrivals would be observed. This only affects two SV M W estimates.
In practice, reported M W values are usually calculated as an average of the values from P-and S-phases which, in turn, are averaged over all sensors available for an event. If such average values are taken, the maximum difference between M W values estimated using R calculated from the focal mechanism solutions and using the average values for the radiation pattern is reduced to 0.4 units (triangles in Fig. 5), which is a smaller, but still significant, difference than is obtained using values estimated for individual phase arrivals.

Time window selection
We consider how the time window (t 2 − t 1 ) chosen to represent the signal affects M W . First, M W is estimated using time windows picked individually for each recording (as illustrated in Fig. 3a) so the length of these windows varies from recording to recording. These are the so-called baseline values. M W is also estimated using different fixed-length windows. We test P-phase windows 0.02 s, 0.03 s and 0.04 s in length. The maximum fixed P-phase window length was determined by the minimum arrival time difference between P-and S-waves, t sp . We also test S-phase windows 0.03 s, 0.08 s and 0.10 s in length.
Most P-phase estimates are within 0.2 units, irrespective of the time window length and, on average, longer time windows give slightly smaller magnitude estimates than the shorter time windows for S-phases (Fig. 6). The P-and Sphase estimates each have one outlier, up to 0.4 units. These occur because the arrivals are impulsive and the longer time windows do not include significant extra event energy but reduce SNR in the selected time window. These occasional outliers highlight why it is important that time windows are carefully chosen. The maximum differences using different time windows are nonetheless smaller than the differences between methods (Fig. 4) or for different radiation pattern corrections (Fig. 5). Note also that, based on the SNR criterium of SNR > 3.0, 98 estimates of M W are available for individual time windows but this number decreases with increasing window length.

Recording frequency
Intuitively the recording frequency should not make a significant difference to the estimated low frequency spectral level; it should only affect attempts to estimate the corner frequency as illustrated in Fig. 7(b). However, analysis of this Cotton Valley dataset is restricted by the P-and S-arrival time separation, t sp . The minimum t sp in the data is 0.04 s and, at a sampling frequency of 250 Hz, this allows only 10 samples to be included in the P-window. A time window containing a very small number ( 20) of samples can severely affect the estimated frequency spectrum: the lower the sampling frequency and shorter the time window, the lower the estimated frequency level (Fig. 7). To test how recording frequency affects M W estimates, the seismograms recorded at 1000 Hz were resampled at 500 Hz and 250 Hz and the estimates repeated for SNR > 1.5 and individually picked time windows. For most estimates, the P-phase values (Fig. 8a) are more affected than those for the S-phases (Fig. 8b), which is expected because the S-phase time windows are longer. P-phase estimates at 250 Hz are up to 0.3 units smaller than estimates at 1000 Hz. The differences between S-phase estimates are <0.1 units. These results show that the recording frequency can have a significant effect on reported M W because signals are not represented correctly at lower sampling rates.
If 0 , Q and f c are to be estimated through spectral fitting methods, then lower recording frequencies have a serious effect on Q and f c estimates because corner frequencies of microearthquakes (M W < 0) are greater than ∼100 Hz and, as f c approaches the Nyquist frequency, f N , the amplitude spectrum cannot be correctly represented or modelled at these higher frequencies. Authors have made different assertions about what f c : f N ratio is required to correctly characterize the source (e.g., McGarr 1984;Baig and Urbancic 2010). It is clear that, as a minimum, the condition f N > 2 f c should be satisfied, and it is recommended that f N > 4 f c .

P -and/or S-waves
The values for M W estimates reported in Fig. 4 are the mean of P-and S-phase estimates. Comparing the P-and S-phase Figure 9 M W estimated using P-arrivals only versus M W estimated using S-arrivals only for M 0 computed using an average radiation pattern correction, R, 0.60 for S-waves and 0.44 for P-waves, and using a correction computed from the focal mechanism solutions given by Rutledge et al. (2004). estimates, computed using the so-called baseline estimate of 0 , illustrates that, in general, S-wave estimates are larger than P-wave estimates (Fig. 9). The maximum difference between the P-and S-phase estimates is 0.8 units, a significant difference. This discrepancy between P-and S-phase estimates arises largely because there are errors in the radiation pattern coefficients. Around half the SV estimates have radiation pattern correction terms <0.20, close to nodes in the radiation pattern compared with the P-waves, and, near nodes, any errors in the radiation pattern coefficient will have a significant effect on M 0 estimates. It is unlikely that the estimated focal mechanism solutions are exact; for example, the events may have a non-double-couple component, and SV estimates of M 0 may be significantly different to the actual M 0 . Larger variations in estimated M W are observed if sensors are near nodal in the radiation pattern and so there is improved correlation between P-and S-wave estimates when average radiation coefficients are used (Fig. 9).
The estimates from S-waves are expected to be more reliable because usually the SNR for S-waves is higher, although, depending on the proximity of the receiver to the event, the S-wave arrival could be affected by scattered P energy. In this study, the P-phase estimates are less variable than the S-phase estimates because the sensors are closer to an antinode in the P-wave radiation pattern than in the S-wave radiation pattern (as noted above in Section 6.2).

Signal-to-noise ratio (SNR)
Using recordings of the event with the highest SNR, seismograms with SNR > 10 between 30 Hz and 300 Hz were available from 13 receivers. Noise levels were increased on these recordings to obtain data with SNR = 10, 5, 3 and 1.5, and M W were estimated for all four sets of seismograms. Figure 10 illustrates how these estimates vary with SNR, with SNR = 10 used as the reference magnitude. There are no large variations in M W with SNR (<0.1 units) because recordings have the same noise spectra with different amplitudes. If noise can successfully be separated from the signal, the SNR does not affect M W . However, lower SNR makes correct estimates of the noise spectrum more important and the estimation of arrival times more difficult, so results are more reliable when using higher SNR recordings.
To assess the potential impact of picking errors on M W , we vary arrival time estimates according to a Gaussian distribution with a standard deviation, σ , of 0.01 s for P-waves and 0.02 s for S-waves. This results in a standard deviation in M W estimates of <0.2 units at individual stations (Fig. 11), much smaller than the difference in estimates between sensors (illustrated by the minimum and maximum M W values in Fig. 10) and between P-and S-phase estimates (Fig. 9).

Number of recording stations
To test how M W estimates vary between stations, the values calculated for M Q W were randomly sampled for between 2 and 10 receivers. The resulting standard deviation in M W values between stations is plotted in Fig. 12. σ approaches a stable value of ∼ 0.02 magnitude units for 4 or more receivers. However, even with fewer receivers, σ ≤ 0.1 magnitude units. Our observed maximum M W between different stations recording the same event is <0.5 units, illustrated by the minimum and maximum M W values in Fig. 10.
The test indicates that a small number of stations can produce a stable magnitude estimate if receivers are positioned in the right place. However, multiple recordings provide increased confidence in estimates. Larger and multiple arrays with good coverage of the focal sphere are also desirable for other reasons, for example to obtain accurate locations and focal mechanism solutions.

Velocity model, location accuracy and rock density
Before a magnitude estimate is even attempted, pre-existing work may affect the accuracy of any future M W estimates.  Equations (6) and (7) require that the velocity model, the rock density and the source location are known. Any errors in these parameters will propagate into the M W estimates. Using a simple propagation of errors given by and assuming a 10% error in velocity, distance and density, the percentage error in M 0 estimates would be 22%. Similar errors at M 0 = 10 8 Nm (M W ≈ −0.7) would affect the M W es-timate by 0.1-0.2 magnitude units. This effect is smaller than some of the parameters discussed above but could prove significant, particularly for traffic-light systems where operations are restricted at a precise magnitude.

Other factors
Other factors, not considered here, may also affect magnitude estimates. We assume here that the events have pure doublecouple mechanisms. This will introduce an error in M W for . Although it is impossible to make any conclusions from one event, the M W estimated in this study are between the two values estimated independently. This particular event has a non-double-couple component that is unaccounted for when using spectral analysis and this highlights the need to move towards routine implementation of methods that take into account non-double-couple source components.
We report magnitudes calculated using spectra corrected for frequency-independent Q. In fact it is likely that Q is frequency dependent (e.g., Campillo, Plantet and Bouchon 1985;Stork and Ito 2004;Oye et al. 2005), but this is not thought to significantly affect our results. The correction for Q in Fig. 4(c) increases reported magnitudes by up to 0.25 units, and any additional correction for frequency dependence is expected to be less than this absolute Q correction.
As mentioned above, Baig et al. (2013) investigate how frequency and bandwidth limitations, instrument coupling and resonance affect interpretation and magnitude estimates. If larger magnitude events than those considered here (M W > 0) are expected, then frequencies <15 Hz are required to allow 0 to be correctly estimated.
In the present study we only consider borehole recordings with short (<1.0 km) source-receiver distances. Understanding attenuation effects becomes increasingly important as the source-receiver distance increases or if the receivers are lo-cated on or close to the surface. Performing the calculations using surface instruments would also require a free surface correction which could be an extra source of inaccuracy. In addition, scattering effects would gain relative importance with greater source-receiver distance.

Summary
The dataset considered here is small but it is fairly typical of a microseismic dataset that allows some general conclusions to be made. The work presented in this study does not place absolute bounds for errors in magnitudes but it is used to identify the most important factors to consider when estimating magnitudes, thus allowing some recommendations for calculating M W . This study also highlights the fact that reported magnitude values are far from error-free.
Overall, the magnitude estimates made in this study show consistency within 0.3 units between the different methods implemented to estimate these values. If Q is known and unless further information is required regarding source parameters, we recommend that M 0 is calculated using the measured low frequency level of the Q-corrected amplitude spectrum. The difference between 0 estimated this way and estimated using spectral fitting is small compared with other errors in M W estimates. If high sampling rates are available and Q is unknown, then spectral fitting provides estimates of 0 and Q. This also requires instruments capable of recording a wide frequency range and well beyond the corner frequency. There is little difference in this study between choosing the Brune or Boatwright model but, in general, the appropriate model is region and data specific, as previously reported (e.g., Abercrombie 1995; Prejean and Ellsworth 2001). The time domain calculation can be used if no Q values are available. However, a correction for attenuation should be made where possible.
This paper also explores the effect of various parameters required to estimate microseismic M W . This is not an exhaustive study but it considers how the choices made in design, acquisition and processing could produce variation in magnitude estimates. We have highlighted some important factors, such as the correction for radiation pattern, that can significantly affect reported microseismic M W and our observations are summarized in Table 1. Some aspects, for example the number of receivers and recording frequency, must be considered carefully in the design stage to ensure sufficient coverage of the focal sphere, correct signal characterization and a sufficient number of recordings to provide reliable locations. In general a distributed coverage of the focal sphere is more important for focal mechanism and event location estimates than the number of receivers, if the number of receivers is >4 (e.g., Foulger and Julian 2011, and references therein). Multiple sensor arrays should therefore be used to provide this distributed coverage. If estimates of seismic energy release, stimulated rock volume and regulatory traffic-light systems are to be based on magnitude estimates then it is important in the interpretation of results that the interpreter understands the inaccuracies in given magnitude estimates. From a regulatory point of view, the method used to calculate magnitudes must be agreed by all parties to avoid disputes over which calculation is correct. Ideally, for a trafficlight system, independent estimates of magnitude should be made.

Recommendations for seismic moment and magnitude estimates
Finally, from the observations made in this study we make the following recommendations for estimating microseismic M W : 1. Use focal mechanism solutions to compute the radiation pattern correction. However, the results presented here show that any inaccuracies in focal mechanism solutions can result in significant variations in M W estimates between sensors and between phases. Therefore values should be averaged over Pand S-phases and the available sensors.

2.
Use M 0 estimated from Q-corrected amplitude spectra with noise removed. 0 should be estimated at low frequencies.
This gives similar results to 0 estimated through fitting with a model spectrum.
3. Use at least 4 receivers to make estimates. The addition of extra sensors does not reduce the errors in magnitude estimates if the sensors are well-distributed, but to ensure good coverage of the focal sphere and good coverage for location procedures and to allow for the occasional malfunction of equipment, >4 sensors are recommended. 4. Good SNR recordings allow precise arrival time picks to be made. In addition, the availability of high SNR data reduces the importance of any inaccuracies in the estimated noise spectrum. It is therefore generally recommended to use SNR > 3 recordings (e.g., Hardt and Scherbaum 1994).

5.
If source-sensor distances are short (<1 km), use recording frequencies ≥1000 Hz for M W < 0 to allow the signal to be correctly represented in the time and frequency domains. If M W > 0 are expected then sensors with eigenfrequencies ≤15 Hz are required to allow 0 to be correctly estimated. 6. Use individually estimated time windows to avoid occasional large errors in M W . If this is not possible due to data volume, check window length is appropriate for a subset of typical events. 7. Understand the limitations imposed by the velocity model or event location errors.

A C K N O W L E D G E M E N T S
We thank Jim Rutledge for providing his event locations and focal mechanism solutions for the data. ALS is funded by a NERC Partnership Research Grant (No. NE/I010904); JPV is a NERC Early-Career Research Fellow (Grant NE/I021497/1).