Speeded simulation of seismicity accompanying mining and hydrofracture

Purpose – At present numerical simulation of seismicity, used in mining and hydraulic fracturing practice, is quite time expensive what hampers its combined employing with observed seismicity in real time. The purpose of this paper is to suggest amean for drastic speeding up numerical modeling seismic and aseismic events. Design/methodology/approach – The authors propose the means to radically decrease the time expense for the bottleneck stage of simulation: calculations of stresses, induced by a large group of already activated flaws (sources of events), at locations of flaws of another large group, which may be activated by the stresses. This is achieved by building a hierarchical tree and properly accounting for the sizes of activated flaws, excluding check of their influence on flaws, which are beyond strictly defined near-regions of strong interaction. Findings – Comparative simulations of seismicity by conventional and improved methods demonstrate high efficiency of the means developed. When applied to practical mining and hydrofracturing problems, it requires some two orders less time to obtain practically the same output results as those of conventional methods. Originality/value – The proposed improvement provides a means for simulation of seismicity in real time of mining steps and hydrofracture propagation. It can be also used in other applications involving seismic and aseismic events and acoustic emission.


Introduction
Since the pioneering work by Salamon (1993), numerical simulation of seismicity made remarkable progress.In the past decade, modeling of seismic and aseismic events has got increasing applications in industry to better interpret observed seismicity, quantify uncertain parameters, play scenarios, improve design and planning, control safety in mines and increase efficiency of hydraulic fracturing (Du et al., 2008;Dobroskok and Linkov, 2006;Cipolla et al., 2011;Kresse et al., 2011;Malovichko, 2013;Malovichko and Basson, 2014;Linkov et al., 2016).
In all the cases, computations are performed by complementing a conventional (basic) code of the finite element method (FEM), the discrete element method (DEM), the finitedifference method (FDM) or the boundary element method (BEM) with a number of "seismic" subroutines.The basic code provides mechanical quantities, in particular stresses, at points of rock mass."Seismic subroutines," which complement the basic code, serve to: generate statistical input data on sources of possible events; simulate seismic and/or aseismic events by using data on stresses on time steps; and process statistical output data on the simulated seismicity.
A basic code, as well as seismic subroutines, has its specific bottlenecks, which may dramatically increase the time expense for calculations.For a conventional code of FEM, DEM, FDM or BEM, it is the stage of solving an algebraic system with large number of unknowns.It is common for any problem involving the need to consider systems with large number of degrees of freedom (DOF), and there are well-developed approaches to overcome the difficulty.Particularly, in methods involving sparse matrices, like finite elements, discrete elements, finite differences and molecular dynamics, the sparsity is used to speed up calculations (Pissanetzky, 1984).For methods involving densely populated matrices, like BEM and those used to solve potential problems, the needed time reduction is reached by using fast multipole methods (FMM), developed in analytical (Rokhlin,1983;Greengard and Rokhlin, 1987;Liu and Nushimura, 2006) and kernel independent (Ying et al., 2004) forms.These methods have already gained great progress, and their development is out of the scope of this paper.It is assumed that some of them are implemented if needed.Still, below we shall use some concepts and findings of the FMM (Rokhlin, 1983;Greengard and Rokhlin, 1987;Ying et al., 2004;Liu and Nushimura, 2006), which may be used to speed up simulation of seismicity.
In this paper, we focus on the bottleneck, which is specific for numerical modeling of seismic and aseismic events and which had not been tackled yet.It is caused by the need to account for the influence of a large number N a of already activated flaws (sources of events) on a large number N na of flaws, which are non-activated but may be activated by events of the group N a (Salamon, 1993;Malovichko, 2013;Linkov et al., 2016).The total number of arithmetic operations on a single time step, being proportional to N a Â N na , leads to extreme time cost on multiple time steps when the numbers N a and N na are of orders used in practical calculations (hundreds, thousands and tens of thousands).
Our objective is to overcome the difficulty by suggesting a means for drastic speeding up numerical modeling of seismic and aseismic events.The purpose is reached by: building a hierarchical tree of octants and special influence lists; and proper accounting for the sizes of activated flaws with excluding check of their influence on remaining flaws, located beyond strictly defined near-regions of strong influence.
Comparative simulations of seismicity by conventional and improved methods demonstrate high efficiency of the means developed.When applied to practical mining and hydrofracturing problems, it reduces simulation time by two orders of magnitude to obtain practically the same output results as those of conventional methods.

Simulation of seismicity
The need in speeding up simulations of seismicity for assessment seismic hazard in mining is clearly comprehended (Malovichko, 2017).Specifically, it is recognized that probabilistic descriptions are presently not feasible due to excessive time cost of repeated calculations with statistically seeded input parameters.Therefore, it is reasonable to develop a means, which reduces the time expense and makes possible multiple simulations needed for the risk assessment.Speeded simulations become also of growing significance for modeling hydraulic fractures (HFs) because the complexity of the problem increases when accounting for inhomogeneity of rocks.
The remaining of this paper is organized as follows.In Section 2, the flowchart, explaining the general line of numerical modeling of seismic and aseismic events, is briefly discussed.Section 3 presents the components of the proposed procedure, which accelerates accounting for flaw interactions.In Section 4, examples of application of the method developed for two practical problems demonstrate its efficiency.The conclusions of the paper are summarized in Section 5.

Problem statement: general line of modeling and its bottleneck
Geotechnical activity such as mining, hydraulic fracturing and oil, gas and heat production leads to changes of stresses in a rock mass.Seismic and aseismic events are the result of these changes.They appear at flaws always present in rock mass, when the tensile or shear strength is reached at flaw location and surface.Then mutual displacement of two shores of the surface occurs in the form of either a jump (seismic event) with elastic energy release, or rapid, while smooth movement (aseismic event).An event itself leads to the change of stresses, which in turn may generate new events, and so on.Thus, chains of seismic events appear, and their totality is called seismicity.Clearly, to simulate seismicity, one needs to repeatedly evaluate stresses caused by two groups of factors: (1) "external", such as mining steps or HF propagation; and (2) "internal", due to activated flaws.
Therefore, as mentioned in Section 1, the simulation may be performed by complementing a conventional basic code of FEM, DEM, FDM or BEM, which accounts for the external factors, with subroutines, which account for the internal factors.The general line of modeling is described in detail in the paper (Linkov, 2013).The flowchart of Figure 1 gives its summary to delineate the bottleneck, which strongly influences the time expense for

Prescribing input data
In accordance with the general structure of the code, the input data consist of two groups.
The first one includes input data required by a basic code.It contains deterministic data on geometry of a region considered; physical properties of sub regions composing the region; and boundary and contact conditions at external surfaces, interfaces and surfaces of cracks.
It also includes data on initial conditions when time is explicitly present in constitutive or/ and in movement or/and conservation equations.For infinite regions, conditions at infinity are assigned.For HF problems, additional data on a fracturing fluid, pumping rate, fracture toughness, etc. are provided.These data are loaded in accordance with an instruction to a particular basic code.For uncertain input parameters to be quantified by comparing observed and simulated seismicity (Malovichko and Basson, 2014), an initial guess is taken (Malovichko, 2013).The corresponding subroutine is shown as BasicInput in Figure 1.
The second group consists of specific input data needed for seismicity simulation.Actually, it concerns with seeding flaws which, when activated by stresses, become the sources of seismic and aseismic events.Detailed description of these data may be found in the key-note lecture (Linkov, 2013).Note that the location and orientation of flaws may be set either randomly (Salamon, 1993) or by using the basic code to delineate local areas of stress concentration, where events are most expected (Malovichko and Basson, 2014), or by combining the two.The input data on mechanical properties of flaws (tensile and shear strength, shear rigidity and softening modulus) are quite uncertain, and they are roughly estimated from available data of seismic monitoring and/or recommendations of geotechnical engineers.Thus merely a reasonable minimal sets of input parameters are used when applying the modeling in practice (Salamon, 1993;Malovichko and Basson, 2014;Linkov et al., 2016).In the flowchart of Figure 1, the corresponding subroutine is called SeisInput.
Besides, exclusion of flaws, which appear activated at the in situ state (without mined areas or HF), is performed before accounting for the stresses induced by mining or HF.This is done by a subroutine called ExcInSitu (Figure 1), which may be considered as a part of prescribing seismic input data.Remaining flaws are re-numerated in the global numeration and their total number N is fixed.Only they are used in further calculations and their total number N is not changed.Later on, in the process of calculations, the totality N of flaws will be divided into two distinct arrays, one of which includes N na non-activated flaws, whereas the other consists of N a activated flaws.At the beginning, all the flaws are non-activated (N = N na ).The groups of activated and non-activated flaws are updated when some of the latter become activated by induced stresses and their changes.Then the number N na decreases, whereas the number N a increases with the sum N = N na þ N a unchanged.

Simulation of events and output processing
For a prescribed input data, a basic code of FEM, DEM, FDM or BEM calculates mechanical quantities (stresses, strains, displacements, fluid pressure, fluxes, etc.) in the area of interest.As mentioned in Section 1, it is assumed that special means are implemented to speed up solution of a resulting algebraic system in cases, when considering a large number of DOF.The basic code gives an information to see the connection between mechanical and seismic quantities and to improve interpretation of seismicity in well-established terms of Simulation of seismicity continuum mechanics.When simulating seismic and aseismic events, it provides tractions at non-activated flaws.
When using a conventional simulation of seismicity, the zero step of calculations, performed by the basic code (Figure 1), starts immediately after the mentioned exclusion of flaws, which appear activated at the in situ state.At the zero step, the basic code evaluates tractions corresponding to the initial geometry.The tractions, calculated by a basic code, are called "induced" to distinguish them from "seismic" tractions caused by activated flaws.Induced tractions are calculated by the subroutine InducedTrac, as shown in Figure 1.
The seismic tractions arise when a non-activated flaw becomes activated.It happens when the total tractions at its surface meet a strength condition.Then displacements of the flaw surfaces generate "seismic" stresses, which strongly influence the state of flaws located at its vicinity.At the beginning of calculations, all the flaws are non-activated, and there is no seismic tractions.However, they appear and should be accounted for after at least a single event is generated.Therefore, at each of time steps, it is necessary to calculate tractions generated by each of the activated flaws at locations and surfaces of each of nonactivated flaws, to add them to the induced tractions and to check whether the resulting total tractions meet a strength condition.If they do, a non-activated flaw becomes activated.This means that an event (seismic or aseismic) occurs.All its characteristics, such as its type, location, orientation of the displacement plane, time and stage of arising, energy and seismic moment, are written in special output arrays.In their turn, the newly activated flaws generate their own seismic tractions, and so on.In this way, chains of events appear at a time step considered.They are simulated in a properly organized cycle of the subroutine shown in gray as SeisTrac in Figure 1.
Clearly, in cycles performed by the subroutine SeisTrac, the influence of N a activated flaws (sources of events) on N na flaws, which are non-activated to the moment, is seen.Recall that the sum N a and N na equals to the total number N of flaws considered.Then with growing number N a of activated flaws, the repeated cycles become of excessive time cost.It is the bottleneck for simulation of seismicity, and the purpose of the discussion in the next section is to present a means to notably reduce the time expense needed to evaluate seismic tractions.
The cycles of the subroutine SeisTrac continue while new activated flaws (events) arise.When there is no new event, a next time step may start.
Then the subroutine, shown as StepGeom in Figure 1, accounts for the change of the geometry and physical quantities such as the fracturing pressure in HF problems.Again, the basic code evaluates induced tractions, which are again used by the subroutine SeisTrac to simulate chains of events.
After the last time step, the output arrays are used by the procedure SeisOutput (Figure 1) for statistical analysis of the simulated seismicity in accordance with similar analysis commonly applied to observed seismicity.Then the comparison of simulated and observed data becomes available with clear understanding mechanical features of each of the simulated events.

Speeding up evaluation of pair-wise influences
As mentioned in the previous section, the main time expense is connected with the subroutine SeisTrac, which evaluates tractions, induced by each of N a activated flaws at the location of each of N na non-activated flaws.Commonly the numbers N a and N na are of the same order, whereas their sum N a þ N na equals to the total number N of flaws considered.This implies that the number of arithmetic operations performed by the subroutine at a single time step, being proportional to the product N a Â N na , may be of EC 35,5 order N 2 .This leads to extreme time cost on multiple time steps when the number N of seeded flaws is of order used in practical calculations (hundreds, thousands and tens of thousands).Thus, there is need to improve the algorithm of checking influences of activated flaws on non-activated flaws.This is done by using the next considerations.
The pair-wise interactions of flaws have the same feature as those of sources in potential problems.Namely, mutual influence of N sources (flaws, particles and charges) decrease with growing distance R between them.In the problem of seismicity simulation, the induced stresses decrease quite fast, as (d/R) 3 , where d is the linear size of a flaw.This suggests following the line of the FMM (Rokhlin, 1983;Greengard and Rokhlin, 1987;Ying et al., 2004;Liu and Nushimura, 2006), which uses the feature discussed to reduce the number of arithmetic operations from the order N 2 to the order N.
In the FMM, the reduction is based on separating two groups of sources.The first group includes sources, which are relatively close to each other.Pair-wise interactions are accounted explicitly only for the particles of this group.As the number of sources in this group is less than the total number of sources, the time expense for accounting interactions between them is much less than that for the totality of all N sources.When N is large, this provides drastic reduction of the time expense because the remaining sources, which are in the second group, are accounted for through their integral influence.The latter is performed by specially designed algorithms of multipole expansions (in the classical FMM) or equivalent densities (in the kernel independent FMM).
Clearly, it is necessary to first distinguish the two groups.It is achieved by successive division of a region, containing sources, into collections of cells (boxes) of decreasing sizes.The division is continued until the number of sources in a cell exceeds a predefined number n l (n l !1), assigned from preliminary numerical tests (Rejwer et al., 2014).In this way, the region with sources is represented as a hierarchical tree, which consists of branches (cells with the number of sources greater than n l ) and leaves (cells with the number of sources not exceeding n l ) at successive levels of the division.Details concerning with the conventional building of the tree and distinguishing needed classes of cells may be found in the reviews by Ying et al. (2004) and Liu and Nushimura (2006).
Similar to the FMM, when simulating seismicity, it is also reasonable to separate strongly interacting (closely located) sources from those interacting weakly (located far enough).This opens the same possibility to explicitly account for pair-wise interactions only for closely located sources.Then similar to the FMM, the time expense drastically reduces.Similar to the FMM, the separation of strongly and weakly interacting sources may be performed by building the hierarchical tree.
Surely, there are differences between the problems to be taken into account.Specifically, in the problem considered, there is no need in multipole expansions (or equivalent and check surfaces), upward and downward traversing an hierarchical tree and, consequently, in corresponding translations.Meanwhile, there is need to agree the number of considered levels of an hierarchical tree and the sizes of objects at each level with sizes of flaws under consideration.Lists, defining neighbors and far regions of objects at hierarchical levels, should be adjusted to the problem.Improvements, which notably decrease the time expense, are to be suggested.Besides, even for conventional operations of the FMM, it is reasonable to present them in terms appropriate to the problem of seismicity simulation.The algorithm developed accounts for the differences and it is presented below in terminology of the problem considered.It is implemented in the computational code as the new subroutines TreeList and FastSeisTrac, shown in the flowchart of Figure 1.The first of the subroutines uses coordinates of flaws which are taken into account when simulating seismicity.It serves to build under the same loop both the hierarchical tree and influence lists.It is used once.

Simulation of seismicity
The second of the new subroutines, FastSeisTrac, now replaces the initial subroutine SeisTrac.(The latter is not used in the speeded algorithm.)The new subroutine FastSeisTrac uses influence lists prepared by TreeList.It serves to economical evaluation of seismic tractions repeatedly on time steps.

Building the hierarchical tree
We take a minimal spatial cube including all the N flaws.Denote the length of its side a 0 .The cube is called an object of level 0. It is divided into eight cubes with sides a 1 = a 0 /2.These are objects of Level 1.They are called children, whereas the starting cube is called their parent.It is assumed that a flaw belongs to a child object if its center belongs to this object.Thus, merely coordinates of central points of flaws are used for building the hierarchical tree.Consequently, the term "flaw" may be used instead of "point."Child objects, which do not contain flaws, are excluded from the tree.In its turn, each non-empty object is divided into eight cubes with sides a 2 = a 0 /4.The latter are objects of Level 2. They are children of the parent cube from which they are obtained by the division.Again, child objects, which do not contain flaws, are excluded from the tree.In its turn, each non-empty object is divided into eight cubes, now with the sides a 3 = a 0 /8.The division process is repeated K times.As a result, we obtain a hierarchical tree of octants with sides: On the l-th hierarchical level (l # K).In the conventional FMM, the total number K of levels becomes known at the end of building the hierarchical tree.Indeed, as mentioned, the building is performed by successive divisions of objects, which still contain the number of sources greater than a predefined maximal number n l .Since the number N of sources is finite, at some stage of divisions, there are no new objects containing more than n l sources.On this level, all the objects are leaves and the division is stopped.The number of this stage gives the number K of the last level of the hierarchical tree.For instance, when the sources are distrusted in the region uniformly, all the objects become leaves at the last level ent[log 2 (N/n e )].Clearly it depends on the choice of the maximal number n l of sources in a leaf.A proper choice of n l is established from preliminary tests (Rejwer et al., 2014).
In difference with the FMM, when using the hierarchical tree for seismicity simulation, the total number of levels K is an input parameter.It is estimated as follows.For the uniform distribution of flaws, the amount of flaws, located in a region V, is proportional to its volume.Consequently, as has been noted, we may roughly estimate the number of flaws in an object on a level l as (1/2 3l )N.However, the number of objects on a level rapidly grows with growing l.For the uniform distribution, it grows as a geometric progression with the factor 2 3 = 8; thus, on the fifth level, the number of objects reaches 262,144 objects.Each of them may contain sufficiently small activated flaws what would require evaluation tractions induced by them in neighboring cells.Clearly, not to have too many objects for evaluation of tractions induced by small flaws, the number K of levels should not be too big.On the other hand, in the problem considered, the side a 0 /2 K of an object at the last (K-th) level has to be less or at least comparable with the average linear size d av of seeded flaws.Therefore, the reasonable choice of the maximal number of levels is K = ent[log 2 (a 0 /d av )] þ 1. Commonly K is in the range from five to seven.Numerical examples in the next section show that the suggested choice appears to be optimal.
With the number K of levels fixed, we successively divide the objects onto octants to build the tree and to assign flaws to objects on a current hierarchical level.This process is EC 35,5 modified as compared with that conventionally used in the FMM.When coming to a successive level, we first use special re-numeration of the flaws in cells of this level.The aim and essence of the re-numeration are explained in details in the next section.
3.2 Special re-numeration of flaws Consider a parent-object at some level.For it, we have prescribed the total number M b of flaws belonging to it.The flaws are numbered in growing order from N 1 to N 2 , so that N 2 = N 1 þ M b À 1, and for each flaw, its number in the starting global numeration is known.The parent-object is divided into eight child-objects, numerated from 1 to 8. The M b flaws of the parent cell are analyzed to find the child, to which a flaw belongs.As a result, we attribute a child-object to each of M b flaws, and find the total number of flaws in each of the children.Denote the total number of flaws in the k-th child M k ; obviously, M b = P k = 1,.,8M k.If M k = 0, the corresponding child is empty and it is excluded from further analysis.Flaws of the first non-empty child k 1 with the total number of points M k1 obtain numbers from N 1 to N 1 þ M k1 À 1; points of the second non-empty child (if it exists) with the total number of points M k2 obtain the numbers and so on.Note that in the new numeration, the first element always has the number N 1 , and the last element always has the number N 2 .Hence the re-numeration does not influence the numeration of flaws in other objects on a considered level and on preceding levels.As a result, for each child, which becomes a parent on the next level, the situation is reproduced.The number of a flaw in the new numeration is linked with the number of its number in the initial (global) numeration, which defines the coordinates and properties of the flaw.The new numeration notably simplifies cycles with flaws belonging to an object.
Besides, in the process of building the hierarchical tree level after level, each new object successively obtains its global number.The global number of an object serves to save data on its level, integer coordinates on this level, the number of the parent object and the numbers (after re-numeration) of the first and the last flaws in it.

Lists for near region
In parallel with building the tree, for each object B, we determine its neighbor region.The latter consists of all objects on the same level as B, which have at least one common apex with B. We call the list of such objects Y list and denote it as L B Y .These lists are similar to U lists of the FMM (Ying et al., 2004); they differ from the latter by including objects merely of the considered level.
Figure 2 shows in gray the objects in the list L B Y in 2D (then the maximal number of objects in the list is 9).In 3D, the maximal number of objects is 27.It is evident that an object B belongs to its own Y list, the number of objects in Y list does not exceed 27 and if The last (reciprocity) relation allows us to diminish inspections of geometrical relation between pairs of objects as follows.Starting from the first object on the first level, we assign it to its Y list and run downward simultaneously with making the tree.It is sufficient to compare (and, when appropriate, add to already started lists) each newly created object B with those of preceding objects, which are on the same level as B.
After analyzing the last object on a considered level, Y lists of all the preceding objects become complete.As mentioned, building the tree and preparing lists for objects are performed simultaneously by the subroutine TreeList (the box with dashed contour line in the flowchart of Figure 1).

New subroutine FastSeisTrac
The key element, used to speed up evaluation of seismic tractions, is neglecting the influence of an activated flaw of linear size d on non-activated flaws located beyond a distance R, which significantly exceeds d.As the tractions induced by an activated flaw decrease as (d/R)3 , one may disregard its influence when: where d is a prescribed tolerance.Consider an activated flaw with the linear size d.For each level l of the tree, the flaw belongs to an object B l on this level (l = 0,. .., K).For l ! 2, the distance R between any flaw (including the considered), belonging to B l , and a flaw located beyond the Y list of B l is greater than the size a l of the cube B l : R > a l .Therefore, for such pairs of flaws, we have d/R < d/a l .Hence, by choosing the level l such that the size a l of its objects satisfies inequality: The condition [equation ( 2)] is certainly met.(In the case, when d ! a l ffiffiffi ffi d 3 p for all l ! 2, the level is set as l = 1).Then seismic tractions, generated by the considered activated flaw, may be accounted for merely those non-activated flaws, which belong to the Y list of B l .Since the latter contains at most 27 objects, the number of non-activated flaws in them is just a small part of the total number of non-activated flaws.We may estimate this part assuming, for certainty, that the distribution of flaws in the region considered is uniform.Then the number of flaws in a cell is proportional to its volume.A cell, belonging to the l-th level, has the side a 0 /2 l and its volume is (a 0 /2 l ) 3 .Consequently, the number of flaws in it is (1/2 3l )N.Hence, the number of traction evaluations for flaws in the cell is 1/2 3l -th part of evaluations performed when accounting for all N flaws seeded.For the maximal number of 27 cells in a Y list on the level l, the number of evaluations is at most 27(1/2 3l ) N .Therefore, the time expense for the evaluations becomes approximately 27/2 3l -th part of that when involving all non-activated flaws.For the second level, a rough estimation of the time reduction is 2 6 /27 % 2.4; for the third level, it is 2 9 /27 % 19.0; for the fourth level, it is 2 12 /27 % 150; for the fifth level, the reduction is 2 15 /27 % 1,200; and for the sixth level, the reduction is 2 18 /27 % 10,000.To properly comprehend these estimations, recall that activated flaws, whose influence on nonactivated flaws is checked, occur at various levels.Consequently, the reduction of the time expense for seismicity simulations is always less than that which would be if all activated flaws have appeared at the highest level K used.The examples of practical applications, given in Section 4, show that the actual reduction is of order 100.
It remains to properly assign the maximal level l, for which the sufficient tolerance condition [equation ( 3) is still met.Denote: Then the minimal factor m, satisfying the equation ( 3), is: As mentioned, the tolerance d is assigned.Then equation ( 5 4), a l = md with a l defined by equation ( 1), the level l is found as that for which the size a l of its objects is minimal among the sizes a k , satisfying the condition a k !md.Thus, combining equations ( 1) and ( 3), the level l of objects, containing the flaws influenced by an activated flaw of the size d, should not exceed l ¼ ent log 2 ffiffiffi ffi d . This estimation may be used when l # K. Recollect that the optimal maximal number K of levels is defined above as In practice, for m fixed in accordance with equation ( 5), it is sufficient to assign an appropriate level l (and corresponding object B l on it) to an activated flaw depending on its linear size d.This is done by the subroutine FastSeisTrac, developed by the authors, as a part of evaluation of seismic tractions.With the level l and the object B l known for an activated flaw, the seismic tractions, generated by this flaw, are calculated only for those non-activated flaws which are in the list L B l Y .This is performed by the subroutine FastSeisTrac (Figure 1), which replaces the conventional subroutine SeisTrac.As explained, the time expense notably decreases as compared with conventional accounting for all the non-activated flaws.Tests and examples, presented in the next section, show that with m = 4, the time expense for evaluation of seismic tractions decreases some two orders.What is also of essence, the output results are practically the same as those of conventional calculations.

Examples illustrating efficiency of the method developed
The efficiency of the means developed is illustrated by two examples of seismicity simulation, previously tackled by a conventional approach.The first of them concerns with the mining problem and the second with HF problem.In both cases, the basic code and seeding flaws were as follows.
4.1 Basic code and seeding flaws 4.1.1H-BEM code.The basic code used is the code of the BEM, founded on the 3D hypersingular boundary integral equations (H-BIE) for blocky systems with openings, faults, dykes, cracks and inclusions.The Poisson ratio is assumed the same, p = (close to 0.3) for all the blocks.In the problems under consideration, the tractions are continuous through a contact and the real H-BIE is (Jaworski et al., 2016):

Simulation of seismicity
where S is the total surface of blocks, the normal n to the surface is fixed arbitrary on contacts of blocks and Du = u þu -is the vector of displacement discontinuity.The superscript "þ" ("-") refers to the limit of a value from the side, with respect to which the normal n is outward (inward).The matrix J S is obtained by applying the traction operator T n to the matrix U of a fundamental solution: For the below-mentioned Kelvin's fundamental solution: q is the distance between the field x and integration y points; the matrix J S has the form: The hypersingular kernel is J H (x) = T n U S (x, y), where the matrix U S (x, y) = [ J S (y, x)] T defines the kernel of the potential of double-layer.Emphasize that the matrices J S and J H /m do not depend on the shear modulus.Equation ( 6) is solved by the boundary element method.The total boundary S is represented as a set of four-vertex plain elements S q .In particular cases, two of four vertices may coincide, then an element is a triangle.As shown in the paper by Jaworski et al. (2016), the integrals entering equation ( 6) are efficiently evaluated almost analytically for quite general approximations of densities.Thus, the basic code used is quite efficient.
4.1.2Seeding flaws.The seismic subroutine SeisInput, complementing the basic code, seeds rectangular flaws with centers randomly distributed in a parallelepiped with sizes three-to fivefold greater those of the region of mining or hydrofracture.The first group of seismic input data includes the total number of seeded flaws, their locations and sizes.The uniform random distribution is used for three coordinates of the flaw center, dip angle of a flaw (changing from 0 to p /2), strike angle (from 0-2p ), rotation angle in the dip plane (for a rectangular flaw changing from 0 to p ) and the aspect ratio b/a of a rectangular flaw with the smallest side a and the greatest b (from 1 to 3).For the size a, the exponential distribution is used.The mean size of seeded flaws is correlated with expected level of energy, estimated from observed data of seismic monitoring.The total number of seeded flaws depends on the volume of the parallelepiped and on the mean size of flaws.It is chosen in accordance with theoretical EC 35,5 estimations (Linkov, 2013) to have the density of flaws in the recommended range from 0.14 to 0.75.The second group of input data on seeded flaws includes their mechanical properties, such as shear rigidity, tensile strength, shear strength, shear softening modulus, which depend on a particular problem.They are specified below for the problems of mining and HF.Simulation of seismicity 4.2.3Mining conditions.The mining conditions are also those considered in the paper by Linkov et al. (2016).In particular, the mining depth is H = 870 m, the elasticity modulus of rock is E = 37 GPa, the Poisson's ratio is = 0.25 and the principal vertical stress is s 11 = À21.75MPa.Since the mined-out area is large, the pressure of the roof on the floor is taken into account.4.2.4Initial flaws.The previous simulations for this panel (Linkov et al., 2016) were performed under limitations on the time expense required for conventional accounting for pair-wise influence of flaws.They used seeding of 20,000 initial flaws.With the improved algorithm, it was reasonable to increase the number of seeded flaws to 75,000 to clearly distinguish the gain in time expense.The density of flaws and all other parameters of flaws were taken the same as previously done.A summary of the parameters used for simulation is given in Table I.
4.2.5 Tree and influence lists.The parameters of the hierarchical tree and influence lists, used in the new subroutines, are specified as follows.For the number N = 75,000 of randomly seeded flaws, the hierarchical tree built includes five levels (K = 5).The total number of objects in the structure is 12,066, the last fifth level contains 10,234 objects (cubes) of the size 15.625 m.Approximately, seven flaws are located in each cube of the last level.The accepted tolerance d = 0.0156 corresponds to m = 4.
4.2.6Comparison of conventional and improved algorithms.For most of the flaws, their sizes are quite small.Consequently, for them, when using the algorithm suggested, the corresponding level l was the last level (l = K), and Y lists were made for an object on the last level.Thus, for the flaws of small size, the number of evaluations was reduced from tens of thousands, required by the conventional approach, to first tens.Still, there are also flaws of greater sizes, for which the corresponding levels of influence are less than K.For objects on these levels, having notably greater volume (and consequently greater number of non- It can be seen that the time expense of the method developed is some 100-fold less than that of the conventional approach.What is of essence, a detailed comparison of the output results for the simulated seismicity shows just minor differences in the number of simulated events.The "lost" events present a small portion of the total number of generated events (less than 1 per cent on a time step).The difference is insignificant in view of the statistical nature of both input data on seeded flaws and output processing of the simulated seismicity.We conclude that the means developed provide notable reduction of the time expense under practically the same output results as those of the conventional approach.

Application to hydraulic fracturing
The second example refers to simulation of seismicity accompanying HF propagation.The problem, studied in the paper by Dobroskok and Linkov (2008), is revisited.We apply the accelerated algorithm and compare it with the conventional one.
4.3.1 Input data of the basic code.The input data are those used previously.Specifically, the center of a productive layer is located at the depth H = 500 m from the earth surface.The origin O of Cartesian system is placed in the center of the layer, the axes x 1 , x 2 and x 3 are directed, respectively, upward, to the east and north.These are the directions of principal stresses: s 11 = À7.5 MPa, s 22 = À5.0MPa and s 33 = À4MPa.The elasticity modulus of rock is E = 20 GPa, the Poisson ratio is = 0.25; thus the shear modulus is G = 0.5E/(1 À ) = 13.33 MPa.
A rectangular HF of the vertical height h = 50 m propagates in the horizontal direction x 2 .Thus, its plane is orthogonal to the direction of the minimal rock pressure 4 MPa, acting along the axis x 3 .The propagation is caused by the pressure p of a fracturing fluid pumped at the center O with the pumping rate q 0 = 2 m 3 /min.The net pressure Dp = p þ s 33 is positive.The fluid is Newtonian with the dynamic viscosity m = 3Â10 À8 MPa•s.The leakoff into rock formation follows the Carter dependence with the fluid loss coefficient C = 1.5s²Â10À4 m/Hs.The changes in time t of the net pressure Dp, the fracture length L and the fracture width w are described by the approximate solution:

Simulation of seismicity
The distributions of the totality of the simulated events and those which arise on time steps are quite similar to distributions obtained by Dobroskok and Linkov (2008).Specifically, the totality of events presents a cloud with the final fracture at its central part.The cloud spreading well beyond the final fracture, it hardly may serve to conclude on the sizes of the fracture confidently.In contrast, the distributions on time steps follow the propagating fracture front and they are about symmetric about it.Therefore, their tracing allows one to conclude on the initial and final positions of the front, and consequently it gives information on the final sizes of the fracture.4.3.4Comparison of conventional and improved algorithms.Because of radical difference in mechanical conditions, under which events occur in mining and hydraulic fracturing, the number of seismic events in the problem considered is quite small.It is about 60 events on a time step.Thus the time expense for FastSeisTrac is low; it does not exceed 1 second for a time step.As could be expected, the reduction of the time is not thus significant as in the previous example.Still the time expense became much less than that of non-accelerated calculations, which required about a minute for a step.The gain in the time is some 60-fold.The gain is reached without distortion of the output results.When using the new subroutine, it provides actually the same events as the conventional approach; in the considered ten steps at most, one to two events are lost what is insignificant for the statistical analysis.Finally after ten steps, the total number of events simulated by FastSeisTrac is 570 against 571 simulated by the conventional subroutine SeisTrac.Thus, the method developed provides actually the same results, whereas the time expense is significantly less.

Conclusions
The conclusions of the work are: (1) The analysis of conventional approaches to simulation of seismicity reveals that the most time costly stage is the evaluation of tractions induced by activated flaws on non-activated flaws.
(2) For reducing the time expense, it is reasonable to use the specific feature of pair-wise interaction of flaws: induced stresses decrease as O(1/R 3 ) the distance R between flaws grows.This suggests neglecting the influence of an activated flaw of the size d on non-activated flaws located at distances notably greater than d.
(3) A method developed uses the feature by: building an hierarchical tree for flaw locations; assigning to each activated flaw an appropriate hierarchical level depending on its size and tolerance accepted; assigning the influence list to an activated flaw on its level; and evaluation of tractions, induced by the flaw, merely for those non-activated flaws which are located in objects of the influence list.(4) Numerical implementation of the method developed has shown that it provides significant (two orders) reduction of the time expense as compared with a conventional approach, whereas the output results are practically the same.

Simulation of seismicity
Figure 1.Flowchart of seismicity simulation by complementing a basic code with seismic subroutines Figure 2. The object B of size 3 on Level 3, which contains an activated flaw of size d ) yields the corresponding m.In particular, for tolerance d = 0.01, we have m = 4.64; for d = 0.0156, m = 4. (In numerical examples presented in Section 4, we assume d = 0.0156 and consequently m = 4.) Since by equation (

4. 2
Application to mining 4.2.1 Initial geometry.The input geometry was those previously considered in the paper by Linkov et al. (2016).It corresponds to mining of the seam Fifth in Severnaya Mine of the Vorkuta coal district in April 2013 Figure 3.The panel of the longwall 512-z-V to be extracted is of the length 300 m.The mined-out area of the seam is shown in gray color.Since the seam is flat with small dip angle (less than 5 °), we use horizontal strips composed of square boundary elements (with side 50 m) to approximate the geometry.Still, to increase the resolution, smaller (30 m) square boundary elements are used for the mined part of Longwall 512-z-V.The total number of elements for zero step geometry equals 860.(The total number of unknowns is 2,580).4.2.2Geometry of mining steps.During April 2013, the mining front advanced about 100 m; thus, the rate of coal extraction was approximately 3 m per day.The area mined out in April is shown by the arrow.Seismicity simulations are performed for successive 10 advances of 10 m each.The extracted area for a ten-meter mining step is approximated by 5 Â 15 rectangular elements, which update the geometry.The final number of unknowns was 4,860.
Figure 3.The plan view of the seam Fifth in Severnaya Mine of the Vorkuta coal district used for approximation of geometry Dp D , length L D , time t D and opening w D are defined as Dp D = Dp/p n , t D = t/t n , w D = w/w n and L D = L/L n , respectively, with the Nordgren (1972) normalizing quantities:

Table I .
flaws), the reduction of the time expense is much less.The resulting reduction is illustrated in TableII. activated