Caving 2018 – Y Potvin and J Jakubec (eds)

© 2018 Australian Centre for Geomechanics, Perth, ISBN 978-0-9924810-9-4

https://papers.acg.uwa.edu.au/p/1815_50_Jonsson/

Evaluating the effect on seismicity of a hydraulic fracturing

trial using Bayesian data analysis

K Jonsson Luossavaara-Kiirunavaara Aktiebolag (LKAB), Sweden

J Martinsson Luleå University of Technology and Luossavaara-Kiirunavaara Aktiebolag (LKAB), Sweden

Abstract
During the summer of 2017, a hydraulic fracturing trial was performed in the crown pillar of one of the
orebodies in LKAB’s mine in Malmberget. Mucking in the orebody beneath the crown pillar recommenced in
September 2017, after being shut down since April 2015. The mucking was shut down in 2015 because a
number of large seismic events had caused significant disturbances in the nearby town, with vibration up to
35 mm/s. A research project was started with the final goal to start up production again with less
disturbance to the town. Hydraulic fracturing was determined to be an interesting method for improving
caveability of the crown pillar, and a trial was set up. The trial consisted of partially fracturing three
boreholes drilled from the ground surface into the crown pillar, which is between 550 and 650 m thick at the
trial location. Hydraulic fractures were placed with 2.5 m spacing at the bottom 50 m of each hole.
As part of the evaluation of the effectiveness of the trial, the seismicity before and after the hydraulic
fracturing trial has been compared using standard Bayesian techniques to infer differences. The analysis show
a difference in b-value before and after the trial, which is demonstrated by the parameter 𝛥b having a positive
value. The 95% high density interval (HDI) of 𝛥b excludes the value of zero, which gives confidence in the
conclusion. The difference in b-value depends on the choice of magnitude of completeness, so it is
recommended that the analysis is repeated with a method using all recorded events, i.e. using the full dataset.
Keywords: hydraulic fracturing, seismicity, Bayesian data analysis

1 Introduction
Deep mining in hard, brittle rock can be challenging. With high stress and seismicity impacting on mining at
LKAB’s Kiruna and Malmberget mines, the management of mining-induced seismicity is crucial in order to
continue deep mining.
The Malmberget underground mine consists of about 20 orebodies, of which 14 are mined at present.
Some orebodies have been mined by surface mining, while mining in others started some distance below
the ground surface, leaving a crown pillar. In the Printzsköld orebody, significant seismic events occurred
in the crown pillar causing large vibrations in the town of Malmberget, located adjacent to this orebody.
Hydraulic fracturing has been successful elsewhere in enhancing caveability (Lowther et al. 2016), so a
trial was performed during the summer of 2017, inducing 93 fractures in three holes drilled from the
ground surface.
While the hydraulic fracturing was ongoing in dayshift, mucking recommenced in the afternoon shifts in
drifts located below the fractured volume. After the trial was completed in late August, mucking was
gradually allowed to go up to full production. The working areas both on the ground surface and
underground were continuously monitored.
This paper comprises a comparison of the seismicity before and after the hydraulic fracturing trial to
determine if the procedure has changed the properties of the rock mass in the crown pillar.

Caving 2018, Vancouver, Canada 651

Evaluating the effect on seismicity of a hydraulic fracturing trial K Jonsson and J Martinsson
using Bayesian data analysis

2 Background and description of Bayesian model

Analysing seismic data involves handling large amounts of data, and statistical methods are well-suited for
comparing two large datasets with each other. In this study, Bayesian data analysis was chosen over more
traditional methods since it allows modelling of the posterior distribution of several parameters and gives
possibilities for richer inference. It is also possible to account for the uncertainty or noise in the analysis.
Kruschke (2015) has defined the five basic steps of Bayesian analysis as identify the data, define a
descriptive model, specify a prior, compute the posterior distribution, interpret the posterior distribution,
and check that the model is a reasonable description of the data.

2.1 Data selection

To reduce the number of uncertainties between the studied time periods, a volume was defined (Figure 1)
from which the events were exported for the chosen time periods. The time periods were selected so that
the production rate in the orebody was the same, although the production sequencing was different. It has
been noted before that it is mucking and not blasting that controls the seismic activity in the crown pillar.
The before-dataset contains 23,786 events and the after-dataset 20,569 events. These datasets are used as
input to the modelling. There were a few sensors added to the system between the two time periods,
which has been accounted for by using only events above the magnitude of completeness. This precaution
eliminates the effects of change, such as the sensitivity of the system. When comparing the two time
periods, there is no obvious migration of events, so the assumption of similar conditions holds true.

Figure 1 Studied volume (blue shape). The grey shape is an outline of the cave estimated from the
seismicity, the red lines are the fracturing holes drilled from the ground surface

2.2 Define a descriptive model

Seismic data is often plotted in a Gutenberg–Richter (GR) diagram where the slope of the curve (the
b-value) indicates the relative frequency of small and large events. A higher b-value means less extreme
events can occur. In this study, we want to investigate if there is a difference in the b-value for the two time
periods, before and after the hydraulic fracturing trial. The geographical spread of the geophones used for
seismic monitoring means that the seismic system has varying sensitivity in different areas of the mine, and
that the sensitivity changes over time as the system is expanded or geophones are lost. Small events have
high frequencies and they are therefore often under-recorded by the seismic system, which means that the

652 Caving 2018, Vancouver, Canada

 Seismicity – seismic hazard case studies

linear GR law does not hold for small magnitudes. This problem is usually handled by discarding all events
of a magnitude lower than a specific threshold, called the magnitude of completeness (mc). The value of mc
can be determined by visual inspection of the GR curve, or be provided by software such as mXrap (Harris &
Wesseloo 2015) (Figure 2). For this dataset (before fracturing), mXrap determines mc to be -1.6, which
leaves 12,679 out of 23,786 events in the dataset. This means that almost 50% of the total number of
events are not used for analysis. Choosing mc to be -1.7 gives 15,330 events in the dataset and mc = -1.8
gives 18,254 events. Determination of the b-value is very sensitive to the choice of mc (see e.g. Martinsson
& Jonsson 2018). For the dataset after fracturing, mXrap determines mc to be -1.8 (Figure 3).

mc

b = 0.96

Figure 2 Frequency–magnitude plot showing the magnitude of completeness, mc, and b-value for the
before-fracturing events

mc

b = 1.02

Figure 3 Frequency–magnitude plot showing the magnitude of completeness, mc, and b-value for the
after-fracturing events

Caving 2018, Vancouver, Canada 653

Evaluating the effect on seismicity of a hydraulic fracturing trial K Jonsson and J Martinsson
using Bayesian data analysis

The GR curve is plotted on a log scale, which gives a linear behaviour above the magnitude of completeness
(mc). Plotting the observed magnitudes without the log scale gives an exponential curve, which can be
described by:
𝑁 = 10(𝑎−𝑏𝑚) (1)
where N is the number of events having a magnitude ≥ m, and a and b are parameters describing the
activity and b-value, respectively.
Above the magnitude of completeness threshold, the magnitudes follow an exponential distribution.
𝑓(𝑦, 𝜆) = 𝜆𝑒 −𝜆𝑦 (2)
where λ is the parameter, y = m-m and the relationship to the b-value is
c

b = λ/ln(10) (3)
In this study, the comparison was done by adding a parameter 𝛥𝜆 to the 𝜆 of the after-fracturing group of
data. If 𝛥𝜆 is zero, there is no difference between the before and after groups of data, but if 𝛥𝜆 is different
from zero, there is a difference. A positive 𝛥𝜆 indicates that λ, and hence the b-value, has increased for the
after group.
The analysis was done using PyStan (PyStan 2018), where Stan (STAN 2018) is a statistical package, and
Py is short for the programming language Python. Stan uses Hamiltonian Monte Carlo simulation which
makes it more effective for large or complex samples compared to other samplers. In Stan, the model is
built up according to a fixed structure. Figure 4 illustrates the Bayesian statistical model using the
graphical model representation proposed by Kruschke (2015). The data, yi = mi+mc, consist of magnitude
values of events of events exceeding mc from before and after the fracturing trial, to which the value of
magnitude of completeness has been added to shift the distribution along the x-axis. The data is
modelled by an exponential distribution, whose shape is determined by 𝜆 and the possible difference 𝛥𝜆.
The possible difference 𝛥𝜆 between the groups is multiplied by a time period index (s[i]) depending on
whether the events belong to the before-fracturing group (s[i] = 0) or the after-fracturing group (s[i] = 1),
for the ith magnitude.

Figure 4 Bayesian statistical model of the dependencies between the data and the parameters. The data
is truncated at mc and then translated as yi = mi+mc. so that a standard exponential distribution
can be applied in Stan

The magnitude of completeness (mc) was chosen to be -1.7 to ensure enough recorded events were used to
inform the model. The model was set to run a chain of 100,000 samples.

654 Caving 2018, Vancouver, Canada

 Seismicity – seismic hazard case studies

2.3 Specify prior

In modelling, the choice of prior can have a great influence on the results. A vague prior distribution means
the data has a greater influence on the posterior distribution than the prior, while an informed
(well-defined and specific) prior gives the opposite. Prior distributions have to be assigned for the two
parameters 𝜆 and 𝛥𝜆. In this model, the prior distribution for 𝜆 has been assumed to be a uniform
distribution between 0 and 1,000, with the restriction that the 𝜆 has to be greater than zero. The prior for
𝛥𝜆 is also assumed to be a uniform distribution that can have values between -1,000 and 1,000.

2.4 Posterior distribution

Running the model produces posterior distributions for the two groups, before and after, as well as credible
values of the parameters 𝜆 and 𝛥𝜆. A way to check the uncertainty of the parameter values is to determine
the span of values that are most credible and cover 95% of the distribution (Kruschke 2015). This is called
the highest density interval (HDI). Another important check of the model is the effective sample size (ESS),
which indicates how stable and accurate the representation of posterior distribution is based on the
samples generated. Out of 100,000 posterior samples using Stan, the ESS for both 𝜆 and 𝛥𝜆 are above
45,000. Kruschke (2015) recommend that the ESS should be at least 10,000 to get an accurate
determination of the 95% HDI.

3 Results and discussion

The output from the model is posterior distributions for all parameters. All results discussed here are in
relation to the chosen magnitude of completeness, mc = -1.7. The results are shown in terms of b-value,
which has been determined from Equation 3. The sensitivity of the results to the choice of magnitude of
completeness is also discussed.

3.1 Evaluation of model

The difference between the groups is shown by the parameter 𝛥b. The posterior distribution for 𝛥b is
shown in Figure 5. The mean value of 𝛥b, for mc = -1.7 is 0.096, and the 95% HDI interval is 0.072 to 0.121,
which clearly excludes zero.

Figure 5 Posterior distribution of 𝛥b = 𝛥λ/ln(10) for mc = -1.7

The effect on the b-value is shown in Figure 6. The mean of the b-value for the before fracturing events is
0.928 with the 95% HDI ranging from 0.913 to 0.942. The mean of the b-value for the after-fracturing
events is 1.024 with the 95% HDI ranging from 1.004 to 1.044.

Caving 2018, Vancouver, Canada 655

Evaluating the effect on seismicity of a hydraulic fracturing trial K Jonsson and J Martinsson
using Bayesian data analysis

Figure 6 Posterior distribution of b-values before (black) and after (red) including lines indicating mean
and 95% HDI limits

Figure 7 shows the exponential distribution for the before (black) and after (red) groups, plotted together
with the corresponding histograms of the input data (dots). The two panels are identical except for the
difference in vertical scale. The log-scale representation clearly shows a difference in slope of the two
groups. The curves in corresponding colour represent the estimated exponential distribution of the two
groups. The thin lines are credible distribution curves generated using a random subset of 1,000 posterior
samples from the histogram in Figure 6, and the bold lines mark the distribution using the sample mean.

Figure 7 Exponential distribution for the frequency of events–magnitude relationship before and after
the trial (lines), plotted together with measured data (dots). The bold lines represent the
exponential distribution using the corresponding mean b-value of the two groups in Figure 6,
and the several thin lines surrounding the bold line are credible distributions using a subset of
1,000 random samples from the posterior distributions in Figure 6

656 Caving 2018, Vancouver, Canada

 Seismicity – seismic hazard case studies

The exponential distributions based on a random subset of credible parameter values show minor deviation
from the distributions obtained using the sample means in both groups. The effect of the estimation
uncertainty of the b-value is best spotted in log scale at larger magnitudes. To conclude, we can say that
the posterior uncertainty of the b-value shown in Figure 6 has minor impact on the deviation of the
GR curve.
Note that Figures 2 and 3 are cumulative histograms and Figure 7 is plain histograms of the magnitudes with
bin size 0.1. The cumulative operator in Figures 2 and 3 will accumulate all errors starting from the left, and
hence more errors are summed up towards the right of the data in the figures. The effect of this accumulation
of errors is shown as a ‘random walk’ behaviour of the data in these figures, and hence we should be careful
when drawing conclusions regarding the validity of the GR model based on cumulative histograms.
Conclusions regarding the validity of the GR model to the data should instead be based on the
non-cumulative histograms shown in Figure 7, and specifically in linear vertical scale (i.e. the top panel in
Figure 7). When drawing conclusions based on the bottom panel in logarithmic scale, we must have the
logarithmic scale in mind.

3.2 Sensitivity of model

The same model was run for different values of mc to determine the effect on b-value and 𝛥b. Figure 8
shows b-values before and after the fracturing trial for a range of magnitude of completeness between -1.9
and -1.4. For this range of mc, there is a difference in b-value for all mc tested, and the 95% HDI for the
before and after value does not overlap for any mc, indicating that hydraulic fracturing has the desired
effect on the b-value at reasonable values of mc.

Figure 8 Difference in b-value before (black dots) and after (red dots) fracturing trial for -1.9 < mc < -1.4,
also showing 95% HDI as vertical lines

Figure 9 shows 𝛥b-values for different values of magnitude of completeness. The smallest difference occurs
for mc = -1.55, but for no choice of mc does the HDI of 𝛥b include zero. It can also be noted that for low mc,
the HDI is smaller than for large mc. This increased uncertainty in the estimate is caused by the fact that a
smaller dataset is used to inform the model than for smaller mc.

Caving 2018, Vancouver, Canada 657

Evaluating the effect on seismicity of a hydraulic fracturing trial K Jonsson and J Martinsson
using Bayesian data analysis

Figure 9 Difference in 𝛥b = 𝛥𝜆/ln(10) for -1.9 < mc < -1.4, also showing 95% HDI as vertical lines

In summary, the analysis shows that there is a difference in seismic behaviour when studying the time
periods before and after the hydraulic fracturing trial. The sensitivity analysis with regard to choice of
magnitude of completeness suggests that this value should be chosen with care. Regardless of choice of m c,
the 95% HDI 𝛥b (or equivalently 𝛥λ) did not in any case include zero, so the conclusion that there is a
difference is valid.

4 Conclusion
The 95% HDI for 𝛥b did not include zero for any value of mc. This gives confidence in stating that there is a
difference in the seismic behaviour of the crown pillar after the hydraulic fracturing trial. Some large seismic
events are still expected to happen, but not as frequently as before. The magnitude of completeness has an
effect on the results. The choice of mc = -1.7 is a compromise that gives a reasonably large dataset to inform
the model, reasonable accuracy in the result, and gives reasonable computational time.

5 Recommendations
Since the estimation of the b-value is sensitive to the choice of magnitude of completeness, analysis of the
same dataset should be done using the method outlined in Martinsson and Jonsson (2018) where all events
are included. That would give a more robust estimation of the b-value. The analysis should be repeated
when a larger tonnage has been mucked from the drifts below the hydraulically fractured volume in the
crown pillar. Only then can the question of whether a large enough volume was fractured in the trial be
answered. Another study that could be conducted is if the relationship between production and the seismic
activity in the crown pillar has changed after the hydraulic fracturing trial.

References
Harris, PC & Wesseloo, J 2015, mXrap, computer software, Australian Centre for Geomechanics, The University of Western
Australia, Perth, https://www.mXrap.com
Kruschke, JK 2015, Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan, 2nd edn, Academic Press, Cambridge.
Lowther, RJ, Olivier, L, Lett, JL & Brunton, I 2016, ‘Implementation of a surface-based hydraulic fracturing program to successfully
propagate a large cave through hard, competent near-surface rock masses to achieve breakthrough’, Proceedings of
MassMin 2016, The Australasian Institute of Mining and Metallurgy, Melbourne, pp. 83–95.
Martinsson, J & Jonsson, A 2018, ‘A new model for the distribution of observable earthquake magnitudes and application to b-value
estimation’, IEEE Geoscience and Remote Sensing Letters, vol. 15, no. 6, https://dx.doi.org/10.1109/LGRS.2018.2812770
PyStan 2018, PyStan, computer software, https://pystan.readthedocs.io/en/latest/index.html
Stan 2018, Stan, computer software, http://mc-stan.org

658 Caving 2018, Vancouver, Canada