Numerical modelling of microseismicity associated with longwall coal mining

Microseismicity has long been a precursor for underground mining hazards such as rockbursts and coal and gas outbursts. In this research, a methodology combining deterministic stress and failure analysis and stochastic fracture slip evaluation, based upon the widely accepted fracture slip seismicity-generation mechanism, has been developed to simulate microseismic events induced by longwall mining. Using the built-in DFN facility in FLAC, discrete fractures following a power law size distribution are distributed throughout a 3D continuum model in a probabilistic way to account for the stochastic nature of microseismicity. The DFN-based modelling approach developed was adopted to simulate the evolution of microseismicity induced by the progressive face advance in a longwall top coal caving (LTCC) panel at Coal Mine Velenje, Slovenia. At each excavation step, global stress and failure analysis with reference to the strain-softening post-failure behaviour characteristic of coal, and fracture slip evaluation for microseismicity are conducted sequentially. The model findings are compared to the microseismic event data recorded during a long-term field monitoring campaign conducted at the same LTCC panel. It was found that the released energy and frequency-magnitude distribution of microseismicity are associated with the slipped fracture sizes and fracture size distribution. These features for recorded microseismic events were fairly constant until a xylite rich heterogeneous zone ahead of the working face was approached, which indicates that fractures within the extracted coal seam follow the same size distribution. The features obtained from modelled microseismic events were consistent over the production period, and matched well the field observations. Furthermore, the model results indicate that the power law fracture size distribution can be used to model longwall-mining-induced microseismicity. This model provides a unique prospective to understand longwall coal mining-induced microseismicity and lays a foundation to predict microseismicity, or even rockburst potential in specific geological realisations.


Introduction
Mining-induced seismicity is an ever present problem in most mines, especially at great depths. Underground rocks and coal strata are stressed by large overburden pressure and tectonic stresses. The sequence of mining excavation perturbs the initial stress field, which may trigger the formation or re-opening and movement of microfractures around excavation openings and beyond, manifesting as microseismic activities. In this process, the stored strain energy in the immediate vicinity of the event focal region is radiated (Cook, 1963). Most of the energy is dissipated in plastic deformation of rock material or transformed to heat or sound, whilst the remaining part is released in the form of seismic waves and propagates in the subsurface. The vibration induced can be recorded by geophones and marked as one microseismic event. In adverse geological environments, the nucleation and connection of large amounts of neighbouring microseismic emissions may give rise to catastrophic failure of the rock mass in a violent manner, such as rockbursts.
Extensive research has been conducted to investigate the potential source mechanisms of mining induced microseismic events and the findings (mainly in hard rock mines) have suggested that they are predominantly caused by slippage of pre-existing fractures or flaws under the combined effects of tectonic and mining-induced stresses (Stiller et al., 1983;Rorke and Roering, 1984;Spottiswoode, 1984). Seismicity as a result of slip along weak planes gives normal seismic signatures as observed in earthquakes, whereas those resulting from crushing of rock material under high stress state generates "anomalous" seismic signatures (Ryder, 1988;Wong and McGarr, 1990).
Mining experience indicates that microseismic events display close correlation with mining activities in both space and time (Sato and Fujii, 1988;Styles et al., 1988). In view of this correlation, numerical simulation has become a promising means for the assessment and even the prediction of mining-induced microseismicity. It has been shown that, using conventional stress and rock failure modelling, the predicted local overstressed zones (Senfaute et al., 2001;Abdul-Wahed et al., 2006) and the damaged and failed zones (Read, 2004;Vatcher et al., 2014) have a strong spatial correlation with the intensive microseismic zones monitored in the mines. The shortcoming of the conventional modelling approach is that microseismic events generated at a distance away from the excavation openings could not be represented.
A brief review of the advances made in modelling methods developed to simulate both mining induced microseismic events and acoustic emission (AE) generated in laboratory rock failure experiments is presented in the next section. Resulting from the generation of slippage or the propagation of fractures, both microseismicity and AE are low-energy seismic events (Cai et al., 2007), even though acoustic emissions in rocks have a much higher frequency (10 4 -10 5 Hz) than that of the seismic motion wave frequencies (10 2 -10 3 Hz, audible by human ears). Many AE simulation methods have been developed to reproduce and interpret laboratory experimental results, including continuum damage models, particle-based discrete element method (DEM) models, lattice models, and the combined finite-discrete element method (FEM/DEM) models (Lisjak et al., 2013).

Advances in modelling of seismic activities in rock formations
Broadly speaking, four approaches may be identified with respect to the source mechanisms for the generation of microseismicity (microseismicity potential or fracture slip) and its model representation (elemental or explicit fracture) in rock formations. The first approach can be traced back to the 1990s, where an evaluation index was introduced to estimate the microseismicity potential at all the grid blocks in the model domain representing an isotropic medium. Fujii et al. (1997) used the total maximum shear seismic moment of fractured elements during one mining step as the indicator in three coal mining case studies, and successfully predicted burst-susceptible areas, and even rockburst occurrences. Nazimko and Babenko (2011) proposed an empirical criterion comprising of equivalent stress and stress relaxation rate, in order to account for dynamic failure of rock mass effectively. The spatial distribution of experimental and simulated events demonstrated good agreement. Shnorhokian et al. (2014) constructed four spatial clusters in a FLAC 3D model based on 24 events within three geological formations, and used the maximum shear stress, brittle shear ratio (BSR), and the continuous change in differential stress (CC-DS) as potential criterion of mining-induced seismicity. They found that CC-DS achieved the best result by conducting a quantitative evaluation for all the locations in the model. A major limitation of this approach, which assumes isotropic rock properties, is that the predicted microseismic events are usually located in close proximity of mine openings, and thus the spatial dispersion nature of seismic events cannot be realised.
The second approach aimed at taking into account the presence of pervasive heterogeneity in rock and its effect on the occurrence of microseismicity during mining, or laboratory rock failure tests. Tang (1997) used a Weibull distribution to describe the heterogeneity of rock strength in a finite element code RFPA, and simulated seismicity generated during the progressive failure of rock in the laboratory. Fang and Harrison (2002) developed a local degradation model in FLAC to simulate brittle fracture of heterogeneous rock in the laboratory, whereby material strength of zone elements are represented by the Weibull distribution, and the bulk and shear moduli of zone elements degrade once the Mohr-Coulomb failure criterion is reached. Boltz (2014) adopted a strain-hardening/softening ubiquitous-joint model in the same simulator, while the heterogeneity in the rocks is represented by scaling-down rock properties of the elements representing the joints in the model. One limitation of this approach is that the size of the fractures generated is constrained by the size of the elements, which results in the simulated energy release falling into a much narrower range compared with that typically recorded in laboratory experiments.
The third approach, which simulates large displacements and represents rock fracture phenomena explicitly using discrete element method (DEM) or its coupling with finite element method (FEM), allows micro-fracture nucleation and propagation during rock failure process to be modelled in a more realistic manner compared with the continuum mechanics-based methods described above. Young (2000, 2004) used a DEM code, Particle Flow Code, to simulate rock failure, which is represented by the breakage of bonds between neighbouring particles. The clustering of the microcracks created, both in space and time, is deemed to form an AE event. Lisjak et al. (2013) developed a quasi-dynamic technique, whereby the fracture is represented by means of cohesive discrete elements in a FEM/DEM code to model AE activities resulting from brittle failure of rock. Although capable of capturing the intrinsic nature of AE activities, this approach is computationally intensive, which limits its application to field-scale models.
The fourth approach was developed to address the knowledge gaps in the simulation of microseismicity generation in field-scale models in a more efficient and realistic manner. In the pioneering work of Salamon (1993), a reasonable approximation of microfractures, along which slip can occur to form a microseismic event according to the Mohr-Coulomb failure criterion, was implemented in the model prior to mining. Board (1994) applied the excess-shear-stress (ESS) criterion for fracture slip, and found that the produced seismic responses, such as slip radii and stress drop, were correlated well with the observed data. Linkov (1997) developed an elasticity-softening-creep (ESC) model to simulate both seismic (unstable) and aseismic (stable) events. Wiles et al. (2001) used monitored field seismic events data, believed to be caused by the slip of an existing fault, to estimate its location and orientation. This inferred fault was then implemented in the numerical model. Malovichko and Basson (2014) introduced a set of sub-flaws to account for the partial slip of a mother flaw, so that rather complex flaw shapes can be modelled.
This last approach has so far been applied mostly to hard rocks, whose response to excavation remains predominantly elastic, and where the induced microseismicity usually concentrates around underground openings such as tunnels and caverns. In contrast, microseismic events induced from longwall coal mining can occur up to 100 m ahead of a working face. Jackson (1984) reported that the majority of microseismic activities are located at a distance 10-80 m ahead of longwall faces, and yet little activity was detected in close proximity of the faces in Polish collieries. Sato and Fujii (1988) found that up to two thirds of the microseismic events were located in a zone 60 m in front of and 20 m behind the longwall face, and that the hypocentre distribution along the direction of advance was about symmetrical with regard to a position 10-30 m ahead of the longwall face.
To the authors' best knowledge, simulation studies of microseismicity resulting from longwall coal mining have so far not been reported in the literature. In this research, a microseismicity modelling approach, which is applicable to longwall coal mining, was developed based on the shear slip mechanism of mining-induced microseismicity. This approach makes use of the built-in discrete fracture network (DFN) capability in the commercial software FLAC 3D to form stochastically distributed hypocentres for generating synthetic seismic events in a 3D continuum model. After each longwall face advance, deterministic global elastoplastic stress analysis is first performed, followed by evaluation for fracture slip at all discrete fractures. The released energy from a microseismic event is estimated by considering both the fracture size and stress drop along the fracture surface. A distinctive feature of the developed methodology is the use of a strain-softening constitutive model to better represent the effect of plastic zones around underground openings on microseismicity distributions. This model provides a tool for numerical interpretation of longwall coal mining-induced microseismicity, and its implications for the assessment of seismic response in different geological realisations.
The model developed during the research reported here was applied to simulate mining induced microseismicity at Coal Mine Velenje in Slovenia, where long-term microseismic monitoring was performed over a twelve-week period at a longwall top coal caving panel. Section 2 below presents the development of the modelling methodology. The recorded data of microseismic events at the study mine is described and analysed in Section 3. Section 4 presents a numerical model of the study area representing the specific mining layout and geological conditions, which was developed to simulate microseismic events induced by the progressive advance of the coal face. The results obtained are compared to the field monitored data in terms of spatial distribution, released energy magnitude and frequency-magnitude distribution in Section 5. In Section 6, the spatial distribution of microseismic events, and the fracture size range and distribution are discussed.

A DFN-based model to simulate microseismicity in longwall coal mining
The microseismicity modelling procedure developed involves three main steps.
The model set-up and initialisation starts with the building of a 3D model representing the working longwall panel and the surrounding strata in FLAC 3D . A group of discrete fractures (DFN) statistically characterised as described in Section 2.1 are implemented into the 3D model using the built-in logic in FLAC 3D . Mechanical and strength properties for stress and failure analysis are assigned to the grid blocks representing different rock materials. The boundary and initial conditions are then applied to the model so as to reach an initial equilibrium state.
Next, the simulation of mining induced stress-redistribution and associated failure analysis step is performed. Progressive longwall mining is simulated by a series of coal excavation steps, each of which typically reflects face advance in one or a few working days. In each excavation step, coal extraction is mimicked by "removing" grid elements (or cells) representing the coal to be extracted by a shearer. This is achieved by assigning null mechanical and strength properties to these cells. In this way, the loads previously withstood by the "removed" cells are redistributed to the surrounding cells and new stress equilibrium is established. Failure analysis is then performed throughout the model domain to determine whether local shear or tensile failure has occurred, which is described in more detail in Section 2.2 below. At the next excavation step, concurrent to the simulation of coal extraction, these cells are "reinstated" to simulate the recompaction of the goaf created by the previous excavation step, and this is achieved by assigning the same cells with the mechanical and strength properties of roof goaf materials.
In the final modelling step, microseismic event generation and associated energy release estimation are achieved. At the end of each coal extraction step, a specially-written subroutine is invoked to check whether slip, and thus a microseismic event, has occurred for each individual fracture based upon the prevailing stress states. This process is described in Section 2.3. The evaluation of energy release that accompanies a microseismic event is presented in Section 2.4.
Steps 2 and 3 described above are repeated in each coal excavation step until the completion of longwall mining simulation in one study area. A flow chart of the modelling procedure is presented in Fig. 1. Note that the model progress from Step 2 to Step 3 is a one-way process, i.e. the initiation of microseismic events are determined by stress states, while it is assumed that the stress states of the continuum rock mass are not affected by fracture slippage. In this respect, the implemented DFN is self-contained, which means that the fracture slip can be evaluated at the end of each excavation step, rather than interacting with model elements during the cycling process. Nevertheless, in view of the mesh size used in the current model and computational efficiency achieved, it was found that the one-way coupling scheme implemented was adequate to represent the key characteristics of microseismicity associated with longwall mining without a significant loss of accuracy.

DFN in FLAC 3D
Using a set of pre-defined attributes, a randomly-distributed network of disc-shaped discrete planar fractures can be created in a 3D continuum model in FLAC 3D by using various random seeds. The set of attributes includes fracture intensity, fracture size range, fracture size distribution (power law, Gaussian, uniform and bootstrapped distributions), and fracture orientation (random or pre-specified).

Fracture intensity
The facture intensity (total number of fractures) is related to the average fracture spacing d: where N is the total number of fractures, and V is the volume of the 3D model domain. Thus, the total number of fractures may be determined if the average fracture spacing data are available from field measurements, such as the use of borehole/scanline surveys, face/bench/outcrop mapping, or 3D seismic data. In the absence of these data, N is first assigned a value, which is then varied through trial and error until a satisfactory match to the recorded field data is achieved in numerical simulations of microseismicity.

Fracture size range
As discussed in Section 2.4, the energy released by fracture slip forming a microseismic event is proportional to the cube of the fracture size. The range of the fracture sizes may thus be constrained by the range of recorded energy release magnitude. Fig. 2a presents an example of the spatial distribution and the energy level, represented by the size of individual circles, of the microseismic events recorded at Coal Mine Velenje over one production week (Si et al., 2015a). The corresponding histogram of the magnitude of the recorded energy is presented in Fig. 2b, where the recorded energy values range from 10 2 to 10 6 J. In addition, a small number of recorded events fell into low energy ranges, which is due to the limited energy detection range of the geophones. This suggest that only microseismic events above the magnitude level with the maximum event counts (>~10 3.8 J) can be completely recorded, from which the corresponding fracture size range may be estimated.

Fracture size distribution
In previous studies, fracture sizes were randomly sampled from a Weibull probability density distribution (Salamon, 1993) or a negativeexponential distribution (Board, 1994). In recent years, the power law relationship, which more closely represents the natural distribution of fractures, has been widely accepted to describe the frequency distribution of fracture sizes (Bonnet et al., 2001), which was also applied in this study: where n(l) is the number of fractures with size l per unit volume, α is the density term which gives the number of fractures with unit fracture size, and a (> 1) is the power law scaling exponent which defines the relative abundance of fractures with different sizes within the population. a usually ranges between 2 and 4. The total number of fractures within the range [l min , l max ] may be obtained by integrating Eq. (2), The cumulative distribution function is then given by As illustrated in Fig. 3, as a increases, the distribution progressively skews to small fracture sizes. In this particular case, the fracture length varies from 1.8 m to 13 m.

Fracture spatial distribution and orientation
The fracture position is generated from a Poisson process, leading to a uniform spatial distribution throughout the volume. The fracture orientation is determined by the dip and dip direction. Either preferential joint orientation or random fracture distribution can be assigned to the rock material, depending on the availability of geological measurements. The latter is used in this study.

Strain-softening post-failure constitutive model for coal
Laboratory triaxial experiments have shown that coal samples generally exhibit strain-softening behaviour after failure under confinement, whereby the strength of coal reduces steadily with further straining until a plateau, i.e. the residual strength, is reached (Fig. 4a). As an example shown in Fig. 4, the residual strength is reached when the plastic strain approaches to 0.001. The Mohr-Coulomb failure criterion is widely used for evaluating rock shear failure. It states that shear failure would occur if the prevailing shear stress on a plane exceeds its shear strength, where μ s is the static friction coefficient, σ is the normal stress on the plane, and c p is the peak cohesion (for intact coal). The corresponding post failure residual strength may be expressed as where μ d is the dynamic friction coefficient, and c r is the residual cohesion (for failed coal). The dynamic friction coefficient is typically assumed to be 0.9 times its static friction coefficient. The Mohr-Coulomb failure criteria for intact and failed coal are illustrated in Fig. 4b. Coal has a relatively low strength and strain-softening post-failure behaviour, therefore, a sizeable shear failure zone is usually formed ahead of an advancing longwall face as a result of the stress redistribution around working coal faces. Si et al. (2015b,c) reported a failure zone extending up to 40 m ahead of the LTCC faces in Coal Mine Velenje, where the thick coal seam is mined. In addition, a front abutment pressure zone is formed ahead of the failure (de-stressed) zone, and the stress disturbance can be discerned as far as 100 m above and 50 m below the working panel (Durucan, 1981).

Fracture slip evaluation
Previous modelling work on mining or rock excavation induced microseismicity has been mainly concerned with hard rocks where rock  mass around underground openings behaved predominantly as an elastic medium (Salamon, 1993;Board, 1994). In longwall coal mining, it has been reported (Jackson, 1984;Yu et al., 2016) that microseismic events tended to occur at a distance ahead of the longwall face, where the coal is likely to remain elastic (intact). In this study, it is assumed that microseismicity is predominantly generated in regions where coal remains intact, and behaves elastically.
The Mohr-Coulomb slip condition is used to evaluate if a fracture has slipped or not. The fracture is considered to have slipped if the shear stress along the fracture plane exceeds its resistance τ sf , which can be described as The resistance to slip, which is provided by the bonding and frictional properties of the fracture surfaces, is given by where σ is the normal stress on the fracture surface, and μ sf and c f are the static friction coefficient and the cohesion along the fracture, respectively. Alternatively, the fracture shear strength in Eq. (8) can be estimated by applying a strength reduction factor β to the intact rock strength given by Eq. (9) = τ βτ sf s The shear and normal stress used in Eqs. (7) and (8) may be calculated using the prevailing stresses on the element closest to the centroid of a fracture, hereafter referred to as the zone element. Salamon (1993) used this method for fracture slip evaluation. This approach is valid if the fracture size is small compared to that of the zone element. Otherwise more test points for slip within a fracture are warranted. Board (1994) employed square-shaped fractures, each containing 9 test points for slip, with 8 of them located on the periphery of the square (4 corner points plus 4 in the middle), and one in the centre. In this study, circular discrete fractures provided in FLAC 3D were used and, following Board (1994), a total of nine test points, eight of them uniformly distributed around the periphery of the circle and one at the centre were adopted (Fig. 5).
The test points within a fracture are examined for slip sequentially, using the prevailing stresses of its zone element (closest to the test point). Slip is deemed to have occurred over the entire fracture face and a microseismic event is generated if Eq. (7) is satisfied at any test point, and the fracture is not considered for slip again in subsequent simulations. The workflow implemented in this modelling study for testing for fracture slip is presented in Fig. 6. The detailed stress computation at a test point is described in Section 2.3.1 below.

Computation of normal and shear stresses for a slip test point
For convenience the test points within a discrete fracture are specified in local coordinates with its origin at (x 0 , y 0 , z 0 ) in the global coordinate system (Fig. 5). The direction cosines of axis OX' (or the strike direction vector) are l 1 , m 1 and n 1 , and those of axis OY' and OZ' (or the fracture normal vector) are l 2 , m 2 and n 2 , and l 3 , m 3 and n 3 respectively. The transformation for the ith test point from its local coordinates to the global ones is given by The normal stresses along the local-coordinate axis directions at the test point can be written as where σ x , σ y , σ z and τ xy , τ xz , τ yx , τ yz , τ zx , τ zy are the normal and shear stresses of its zone element, respectively. The corresponding normal and shear stresses on the fracture plane associated with the test point are then given by The computed normal and shear stresses are used in Eqs. (7) and (8) to evaluate for fracture slippage.

Microseismic energy release
The concept of excess shear stress (ESS) proposed by Ryder (1988), which has been extensively used for numerical evaluation of rupture along weak planes (Salamon, 1993), is used to evaluate energy released during fracture slip. The excess shear stress, i.e., the stress drop during slip, is defined as the difference between the prevailing shear stress prior to slip and the dynamic strength of the fracture surfaces: where τ df = μ df σ is the dynamic strength of the fracture surfaces, and μ df is the dynamic friction coefficient along the fracture. Considering that slip occurs when Eq. (7) is met for the first time for any test point within the fracture, it may be argued that the resistance τ sf is likely to be exceeded by the prevailing shear stress τ only marginally at an excavation step. Thus τ may be approximated by τ sf in Eq. (10).
The stress drop is graphically illustrated in Fig. 7. Due to the small difference between μ sf and μ df , the stress drop is expected to be strongly influenced by the bonding properties of the fracture surface.
The resultant seismic moment M 0 is given by (Salamon, 1993): where υ is the Poisson's ratio of the host rock, and R is the fracture radius.
The local Richter magnitude M L is expressed in terms of M 0 (Salamon, 1993;Board, 1994): The released kinetic energy E k from the fracture slip is related to M 0 by (Salamon, 1993): where G is the shear modulus of the host rock. Substituting Eq. (13) into Eq. (15), one has Eq. (16) shows that the released kinetic energy is proportional to the cube of fracture size (radius) and square of the stress drop, suggesting that it is strongly affected by the fracture size and to a less extent by the rock strength (through stress drop) and mechanical properties.

Analysis of mining-induced microseismic events at Coal Mine Velenje
Coal Mine Velenje is located in the north of Slovenia. The main coal seam is lens shaped, about 8.3 km long, 1.5-2.5 km wide and over 100 m thick in the central part (Markič and Sachsenhofer, 2010). The Velenje Coal is a humic type lignite, the matrix of which is detritedominated, dark brownish to dusky brown in colour. Other lithotypes are composed by varying amounts of xylite, fusinite and mineral matter such as alumosilicates and carbonate minerals.
For efficient coal extraction from the ultra-thick coal seam, multilevel longwall top coal caving (LTCC) is practiced at the mine. Specifically, the entire coal deposit is divided into a number of levels ranging from 10 to 20 m in thickness, and extracted sequentially from the top to the bottom. At each level, the lower section of the seam (3-4 m in thickness) is cut by a shearer under the hydraulic supports, and the top coal (7-17 m in thickness) is allowed to cave by gravity and be recovered in front of the supports (Si et al., 2015b). During coal extraction, the extending canopy is collapsed and the top coal is caved after several cuts at the face, allowing for a steady face advance.
The LTCC panel K. -50/C was selected as the field site for monitoring mining-induced seismic activities. This longwall panel lies at a depth of −350 m, underlain by thick floor coal, and overlain half by an intact clay layer and half by a layer of caved roof goaf, formed as a result of previous coal extraction (Si et al., 2015b). Coal at LTCC panel K. -50/C was mined from November 2010 to October 2011. The average face advance rate was approximately 11-12 m per week, yielding a daily coal production of   7930 t. During the period between 27 April and 30 August 2011, a host of monitoring techniques, including microseismic monitoring, inseam borehole stress and pressure measurements, and ventilation environment monitoring were implemented at the K. -50/C LTCC panel to investigate the dynamic response of the coal seam to face advance.
A 32-channel flameproof automated seismic observation system designed by the Department of Geology and Geophysics of Central Mining Institute (GIG) in Poland was used for the underground microseismic monitoring campaign. A total of eight geophones (four uniaxial and four triaxial) were installed at the panel, placed parallel to the coal face in pairs, with a spacing of around 100 m. Recorded seismic data were collected automatically in these geophones, transferred to surface data loggers via fibreoptic cables, and then processed to obtain source parameters and event locations. The geophones were all located in the same level of the LTCC panel, and only microseismic activities that occurred close to this horizontal plane were accurately identified. In addition, time-lapse active seismic tomography over a 100 m × 141 m area ahead of the advancing face was also performed on the same panel (Si et al., 2015a).
Over 2000 microseismic events were recorded around the panel K. -50/C during the monitoring period. In a previous study by the authors (Si et al., 2015a), it was observed that the range of seismic energy released spanned by up to four orders of magnitude (from 10 2 to~10 6 J). Furthermore, the spatial distribution of weekly microseismic events generally scattered over an area approximately 70-80 m ahead of the advancing face-line (Fig. 8). From the two time-lapse seismic tomography measurement campaigns conducted concurrently, Si et al. (2015a) also noted that there was an anomaly zone with relatively higher velocities ahead of the advancing face. The authors suggested that this anomaly was caused by local heterogeneity in the coal deposit, specifically related to the relative richness of the stronger component xylite in the seam.
In the current work, the microseismic events were further analysed in terms of spatial distribution, frequency-magnitude distribution and seismic energy distribution, focusing on the events recorded during twelve coal production weeks from 23 May to 30 August 2011 (excluding a two-week summer holiday from 18 to 31 July).
In earthquake seismology, the seismicity frequency-magnitude distribution, i.e., relative abundance of small seismic events with respect to large ones, has been found to obey the Gutenberg-Richter relationship (Richter, 1958) where N is the cumulative number of seismic events with an amplitude greater than A, and a and b are constants. The parameter b gives a measure of the scaling of seismicity magnitude distribution, with a higher value of b indicating a higher proportion of small seismic events in the population.
Eq. (17) may be rewritten in log-log coordinates as where M = log 10 A. Thus, b may be estimated by fitting a straight line in a logN − logA plot. Monitoring changes in the value of b has been proposed as a potential precursor to tectonic earthquakes (Seggern, 1980). A reduction in b associated with macroscopic rupture of rock samples in laboratory experiments has been reported (Lockner et al., 1991;Lisjak et al., 2013). It has been reported that the magnitude of b for mining-induced microseismicity typically falls in the range between around 0.6 and over 1.0 (Board, 1994).
A frequency-magnitude analysis was carried out on the recorded microseismicity around LTCC panel K. -50/C at Coal Mine Velenje. It was decided to group the microseismic data in two-week intervals to ensure that there is sufficient number of events in the population. Fig. 9 presents curves fitted on the frequency-magnitude distribution. In order to minimise the influence of incomplete detection of small seismic events and statistical variations due to small samples of large seismic events, the magnitude range on which a linear regression fitting can be applied was taken as the range of regression fitting, which mostly lies between the magnitude −0.6 and 0. The fitted bi-weekly b values were found to be fairly consistent with each other (~1.0) throughout the first eight weeks, until the xylite rich heterogeneous zone is approached during the week 11 to 17 July 2011 (Fig. 10). A drop in b value to~0.8 was noted.
A range of factors influencing the fracturing process of rock masses contribute to the complexity of the seismic energy distribution (Lasocki and Orlecka-Sikora, 2008). Zhang et al. (2013) conducted a statistical energy analysis on different seismic datasets and found that no common distribution patterns (normal, Weibull, exponential and lognormal distributions) could be universally applied to describe seismic energy distribution. Nonetheless, a number of statistical distribution models, such as truncated Pareto distribution (Lasocki, 1992) and Weibull distribution (Lasocki, 1993), have been used to fit the logarithm of energy released from microseismic events taking place in specific geological and mining conditions.  Si et al. (2015a) reported that the weekly microseismic energy released at longwall panel K. -50/C recorded at Coal Mine Velenje followed a Gaussian distribution. A Gaussian distribution is characterised by a mean value μ and standard deviation σ as illustrated in Fig. 11. In the current study, Gaussian distribution curves were fitted to the weekly histograms over the monitoring period between 27 April and 30 August 2011. With the exception of two weeks' data, for which there were too few events (8 to 21 August 2011), a good fit was achieved for 10 weekly histograms. The goodness of fitting using the Chi-square test has shown that a p-value > 0.05 (95% confidence) was achieved for 8 out of the 10 weeks evaluated (except for weeks May 30 to 5 June, and 6 to 19 June 2011). The fitted parameter values and their variation with face advance are plotted in Fig. 12. It was observed that the mean value μ and the standard deviation σ remained fairly constant during the first 7 weeks.
It has been shown that the recorded microseismic event energy is, to a large extent, reflecting the size of slipped fractures (Eq. (16). Therefore, the fitted Gaussian distribution mean value μ may be related to the mean sizes of slipped fractures, and both the standard deviation σ of the logarithmic event energy and the b value for the Gutenberg-Richter relationship reflect the fracture size distribution within the coal seam. Consequently, the consistency in both the σ and b values before the week 11 to 17 July 2011 implies that fracture size distribution within the extracted coal panel was fairly constant before the longwall faceline reached the inferred xylite rich heterogeneous zone.
In view of the above mentioned consistency in σ and b values, a numerical model built for the LTCC panel K. -50/C can be assumed as homogeneous at the element scale, and a group of discrete fractures characterised by one set of intensity term and size distribution can be applied throughout the model domain. The following sections present the findings of a numerical modelling study carried out at panel K. -50/ C. The model results are compared to the recorded data in terms of weekly histograms of released event energy and the microseismic frequency-magnitude distribution.

Model set-up and generation of a discrete fracture network
A three-dimensional model measuring 300 m × 350 m × 200 m in x, y and z directions respectively was constructed in FLAC 3D to represent the LTCC panel K. -50/C, the floor coal and the overlying strata (Fig. 13a). In the model the coal panel was at a depth of −350 m and had a thickness of 20 m, the top 15 m of which is caved when the bottom 5 m of coal is extracted by the shearer. The longwall panel modelled was 150 m wide, including both the intake and return gateroads. The coal panel was underlain by a 100 m thick floor coal, and overlain by a 20 m thick mixed roof layer, half of which was clay on the intake gateroad side, and the other half was caved roof goaf on the return gateroad side. This mixed layer was further overlain by a clay layer of 50 m in thickness. The model comprised of 280,000 elements, each with a dimension of 3 m × 5 m × 5 m.
Three geological structures, i.e., coal, clay, and failed clay (for roof goaf) were used in the model. Elastic (bulk modulus, shear modulus) and strength (cohesion) properties for coal, clay and failed clay were obtained from laboratory experiments on samples from Coal Mine Velenje, which are presented in Table 1. Strain-softening constitutive model was used to describe post-failure rock behaviour. The residual cohesions of coal and clay represent a reduction of 83% from their intact values. In the post-failure phase, it was assumed that the residual  cohesion of coal and clay is reached when the plastic strain exceeds 0.001. Considering the mesh size dependency of the strain-softening constitutive model, a preliminary sensitivity study was performed to ensure that the mesh size used in the model did not have a notable effect on the model performance.
The in situ vertical stress in the model was induced by the overburden weight (using a density of 2360 kg/m 3 ), and the horizontal stresses were estimated by the Poisson's effect under gravity loading of the overburden. The vertical and horizontal stress gradients were computed based on the Poisson's ratio of the rock types used in the model. The boundary conditions were set in such a way that normal stresses calculated from the horizontal stress gradient were applied to the lateral boundaries, and the bottom boundary was fixed.
The methodology followed in generating the DFN for modelling microseismicity in FLAC 3D has been described in Section 2.1. The attributes of fractures required for model implementation include total number of fractures, their size range, the scaling exponent in power law distribution, and spatial distribution and orientation. Out of these, random fracture positions and orientations were assigned to the rock types used in the study. The total number of fractures and the scaling exponent needed to be determined during the simulation to match the field monitored data. A strength reduction factor of 0.4 was used in Eq. (9) to determine the shear strength of each fracture.
It has been shown that the recorded event energy range may be used to estimate the range of fracture size in the model implementation using Eq. (16). Based on the complete detection range of released energy (~10 3.8 to~10 6 J) from Fig. 14, the fracture radius was estimated to be between 0.9 and 6.5 m. It can be seen from the figure that the minimum complete detection magnitude is around 10 3.6 J for the modelled microseismicity.
It was found that a total fracture number of 45,000 and scaling exponent of 2.7 yielded a satisfactory match of the recorded event energy distribution (Fig. 14). Fig. 15 presents the implementation of DFN in the model domain.

Modelling of longwall top coal caving mining
The model run covered a period of eight weeks, starting on 23 May 2011, when the face was 81 m from the right-hand side model boundary (Fig. 13b), and finished on 17 July 2011. In order to establish the initial stress conditions to account for the effects of previous mining activities, the two development headings driven around the panel were first extracted (Fig. 13a). Next, a pre-modelling step, whereby a 12 m length (average weekly advance distance) of coal at the longwall behind the 23 May 2011 face-line was simulated. A total of eight longwall extraction steps were then modelled (from right to left in Fig. 13b), each extraction step representing one production week during which 12 m of coal (4 elements in the x direction) was "removed" and then "backfilled" in the following week.
A subroutine, specially written to determine the failure state and the stress state of zone elements hosting the fractures, is invoked after an equilibrium state is reached at each excavation step, and whether the shear slip condition is met by each fracture to generate a microseismic event is ascertained. In order to account for the stochastic nature of the DFN, a total of five series of DFN realisations with different random seeds were generated and used in the model.   profiles in front of the longwall face during longwall face advance are plotted for four typical production weeks (Fig. 17). It is observed that, for each production week, the major principle stress σ 1 first increases from 0 to the peak value of over 15 MPa at the interface between elastic and plastic zones, and then decreases slowly to the in situ vertical stress level with the distance from the longwall face. The minor principle stress σ 3 also reaches a peak at the edge of the failed zone, and gradually converges to the in situ horizontal stress level. -50/C during one typical mining week, each microseismic event being represented by one point, with its colour representing the logarithmic magnitude of the released kinetic energy. It can be seen that, due to the relatively higher stresses in the floor compared to the roof, more microseismic events are generated in the floor coal.

Spatial distribution of simulated microseismicity
In order to be comparable with the field monitored data available from Coal Mine Velenje (Fig. 8)

Table 1
Rock mechanical properties used in the K.-50/C LTCC model (after Zavšek, 1993 andSi et al., 2015c   more than100 m ahead of the face-line were filtered out. It was observed that the microseismic events concentrated in the area 40 m to 70 m ahead of the face-line (Fig. 19), whereas the field recorded events were more evenly distributed over a much wider region (Fig. 8). This relative discrepancy is believed to be caused by the omission of the presence of a small scale xylite rich heterogeneous zone in the coal panel, which will be further discussed in Section 6.

Frequency-magnitude distribution and released energy
Fig. 20 presents the frequency-magnitude distribution of the recorded and modelled microseismic events. The event magnitude obtained from numerical modelling ranges between −2.0 and 1.0, which is in agreement with the field observations. In this study, both the monitored and modelled datasets were fitted to Eq. (18), with the corresponding b values calculated and marked in the figure. The figure shows that a relatively good agreement has been achieved between the model results and field observations. Frequency-magnitude distribution and associated b values for other production weeks are presented in Appendix A. The b values for recorded seismicity fall within one standard deviation of those for the modelled seismicity, suggesting the robustness of applying DFN models in simulating microseismicity (Fig. 21). In addition, b values for both field observations and model results only experience slight fluctuations during the modelling period, which means that the size scaling exponent for fractures in the extracted coal seam is quite consistent.
It has been shown earlier that the recorded weekly microseismic energy released around the K.-50/C LTCC panel at Coal Mine Velenje follows a Gaussian distribution. Curve fitting of a Gaussian distribution to the weekly histograms of the magnitude of microseismic energy predicted was also performed. Two examples of this are shown in Fig. 22, with the remaining results being presented in Appendix B. It can be seen that the goodness of fit for the modelled results to the Gaussian distribution is not as good as that obtained for the field monitored data. This is probably due to microseismic events generated from fractures smaller than the minimum fracture size (0.9 m) being ignored in the simulation, which were incompletely detected in field observations, leading to a discrepancy in the number of events falling into low energy ranges.
As shown in Fig. 23, the fitted mean value μ and standard deviation σ values during each week are fairly consistent over the eight production weeks. Moreover, the simulated values are compared favourably with the field observations over the same period.

Stress drop frequency
As shown in Eq. (16), the stress drop during rock failure is an important parameter influencing the released seismic energy. The histograms of modelled stress drops during the first two production weeks are plotted in Fig. 24. It can be seen that two clusters are formed, with the right hand side one ranging from 0.78 to 1.10 MPa, and the left hand side one from 0.28 to 0.45 MPa. Most stress drops fall within the right hand side cluster.
Stress drops are associated with the dyanmic strength of the fracture surfaces, and thus with lithology as presented in Eq. (12). Given that cohesion and friction coefficient values assigned to coal and clay roof are high and are not significantly different, stress drops occurring in these layers form the right hand side cluster in the histograms. In contrast, the left hand side cluster of stress drops represent the microseismic events occuring in the caved roof. In addition, it can be deduced that once failure has occurred in the xylite-rich heterogeneous zone, the associated stress drop, and thus seismic energy released would be much larger than that in the surrounding coal.

Discussion
Since the coal panel modelled was assumed to be homogeneous at element scale, the spatial distribution of small-scale heterogeneity was not considered in the model. Therefore, the plastic zones induced by coal extraction extend into the whole region approximately 40 m ahead of the working face, as shown in Fig. 16. In the present methodology,   only fractures located in elastic areas are evaluated for the occurrence of microseismicity, thus a few events occur in the immediate vicinity of the longwall face modelled. This is consistent with the monitoring results of Jackson (1984) and Yu et al. (2016).
However, small-scale heterogeneity in rock mechanical properties do exist in geological media. Regions with high rock strength in the vicinity of a working longwall face may remain intact, where microseismic events can still be generated. Therefore, the spatial distribution of recorded weekly microseismic events scattered over an extended area ahead of the working face at LTCC panel K.-50/C.
In addition, the xylite rich heterogeneous zone within the coal panel has a large influence on mining-induced microseismic characteristics, the influence of which will be investigated further in the future.
The cumulative number of implemented fractures can be obtained by integrating Eq. (2): Considering that there is no upper threshold for fracture size in natural geological media, Eq. (19) can be reduced to This equation has a similar form as the Gutenberg-Richter relationship (Eq. (17)), which is a power law relationship in linear scale. The scaling exponent for cumulative fracture number is 1 − a, compared to -b for the Gutenberg-Richter relationship. This further confirms that microseismic events and related energy release are closely associated with the fracture systems in the coal seam. It is noteworthy that, due to the limited sensitivity of monitoring equipment, a considerable number of low energy events emitted by the initiation of small-scale fractures are overlooked. Therefore, the histogram of recorded microseismic energy follows a Gaussian distribution, rather than the power law distribution.
Given the complete detection range of released energy, the fracture size range in the geological media can be estimated. In previous studies, Board (1994) used fracture sizes ranging from < 1 m to approximately 50 m. The mean fracture size used by Salamon (1993) reached up to 100 m. In the model described here, the average fracture spacing was around 7.76 m, and the fracture size ranged from 0.9 m to 6.5 m. Following the power law distribution, most fracture sizes tended to be close to the lower fracture size limit. For example, the scaling exponent a used in the current model was 2.7, indicating that 90% of embedded fracture sizes are < 5.94 m. In addition, a smaller a value indicates a larger proportion of large fractures over small ones, and thus a larger fitted mean value μ for the weekly histogram of released energy.

Concluding remarks
Continuous microseismic monitoring of the effects of a LTCC face at Coal Mine Velenje was carried out during a period of twelve weeks. It was found that the bi-weekly frequency-magnitude distribution of microseismic events and the weekly logarithmic energy release can be described by the Gutenberg-Richter relationship and a Gaussian distribution, respectively. It was found that the fitted mean value μ of the Gaussian distribution is related to the sizes of slipped fractures during face advance, and both the fitted standard deviation σ of the Gaussian distribution and b value for the Gutenberg-Richter relationship are associated with the fracture size distribution. These features were fairly constant during the first eight weeks monitoring period (before the advancing face approached a heterogeneous zone), which implied that fractures within the extracted coal panel could be assumed to follow the same size distribution.  In the current research, a methodology developed to simulate longwall mining-induced microseismicity based upon the widely accepted fracture slip seismicity-generation mechanism was used. Using the built-in DFN facility in FLAC 3D , discrete fractures following a power law size distribution were distributed throughout the 3D continuum model of a LTCC panel in a probabilistic way to account for the stochastic nature of microseismicity occurrence. A geomechanical model was constructed using FLAC 3D to model the progressive advance of a longwall face during the first eight weeks, and the developed DFNbased modelling approach was applied to simulate LTCC mining-induced microseismicity. The model was homogeneous at the element scale, and only one set of fracture size distribution was implemented throughout the model domain. At each excavation step, global stress and failure analysis with reference to the strain-softening post-failure behaviour characteristic of coal, and fracture slip evaluation for microseismicity are conducted sequentially.
Model results have shown that a good match was achieved between recorded and modelled microseismic events in terms of spatial distribution, released energy, and frequency-magnitude distribution. The energy released from a microseismic event is shown to be proportional to the cube of fracture size. This suggests that the resulting microseismic energy distribution is expected to be strongly correlated with the fracture size distribution implemented in the model. The energy released from a microseismic event is also affected by the rock strength (through stress drop) and its mechanical properties. The model findings indicate that the power law fracture size distribution can be used to model microseismic generation around the longwall coal face. This modelling approach provides a new alternative to interpret and predict longwall coal mining-induced microseismicity in field applications.