Speeded simulation of seismicity accompanying mining and hydrofracture

Abstract

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 a mean 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.

Citation

Rybarska-Rusinek, L., Rejwer, E. and Linkov, A. (2018), "Speeded simulation of seismicity accompanying mining and hydrofracture", Engineering Computations, Vol. 35 No. 5, pp. 1932-1949. https://doi.org/10.1108/EC-07-2017-0256

1. 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 finite-difference 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.

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.

2. 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 seismicity simulation. It is shown as a gray box. The flowchart also distinguishes new subroutines, developed by the authors to speed up simulation; they are shown as boxes with dashed contours. Particular parts of calculations performed are as follows.

2.1 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_nanon-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.

2.2 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 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 non-activated 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.

3. 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 order N². 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)³, 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² 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. 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.

3.1 Building the hierarchical tree

We take a minimal spatial cube including all the N flaws. Denote the length of its side a₀. The cube is called an object of level 0. It is divided into eight cubes with sides a₁ = a₀/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₂ = a₀/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₃ = a₀/8. The division process is repeated K times. As a result, we obtain a hierarchical tree of octants with sides:

(1) $$a_l=a_0/2^l$$

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₂(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³ = 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₀/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₂(a₀/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 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₁ to N₂, so that N₂ = N₁ + 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 = ∑_{k = 1,.,8} M_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₁ with the total number of points M_k1 obtain numbers from N₁ to N₁ + M_k1 − 1; points of the second non-empty child (if it exists) with the total number of points M_k2 obtain the numbers N₁ + M_k1 to N₁ + M_k1 + M_k2 − 1, and so on. Note that in the new numeration, the first element always has the number N₁, and the last element always has the number N₂. 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.

3.3 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 a cube C is in L^B_Y, then B is in L^C_Y. 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).

3.4 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)³, one may disregard its influence when:

(2) $${\left(d/R\right)}^3<\delta$$

where δ 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:

(3) $${(d/a_l)}^3<\delta$$

The condition [equation (2)] is certainly met. (In the case, when d ≥ a_l $\sqrt[3]{\delta }$ 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₀/2^l and its volume is (a₀/2^l)³. 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⁶/27 ≈ 2.4; for the third level, it is 2⁹/27 ≈ 19.0; for the fourth level, it is 2¹²/27 ≈ 150; for the fifth level, the reduction is 2¹⁵/27 ≈ 1,200; and for the sixth level, the reduction is 2¹⁸/27 ≈ 10,000. To properly comprehend these estimations, recall that activated flaws, whose influence on non-activated 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:

(4) $$m=a_l/d$$

Then the minimal factor m, satisfying the equation (3), is:

(5) $$m=1/\sqrt[3]{\delta }$$

As mentioned, the tolerance δ is assigned. Then equation (5) yields the corresponding m. In particular, for tolerance δ = 0.01, we have m = 4.64; for δ = 0.0156, m = 4. (In numerical examples presented in Section 4, we assume δ = 0.0156 and consequently m = 4.) Since by equation (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\left[{\log }_2\right(\sqrt[3]{\delta }{a}_0/d\left)\right]+1$. This estimation may be used when l ≤ K. Recollect that the optimal maximal number K of levels is defined above as K = ent[log₂(a₀/d_av)] + 1; therefore, if l > K, we set l = K.

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_Y^{B_l}$. 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.

4. 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.1 H-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):

(6) $${\displaystyle \underset{S}{\int }\left(\frac{1}{2{\mu }^{+}}-\frac{1}{2{\mu }^{-}}\right)}{J}^S(x,y)t_n\left(y\right)d{s}_y-{\displaystyle \underset{S}{\int }\frac{1}{2\mu }}{J}^H(x,y)\Delta u\left(y\right)d{s}_y=\frac{1}{2}\left(\frac{1}{2{\mu }^{+}}+\frac{1}{2{\mu }^{-}}\right)t_n\left(x\right),$$

where S is the total surface of blocks, the normal n to the surface is fixed arbitrary on contacts of blocks and Δu = 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:

$${\left(T_{n}U\right)}_{ij}=2\mu \left[\frac{\nu }{1-2\nu }\frac{\partial U_{kj}}{\partial x_k}{n}_i+\frac{1}{2}\left(\frac{\partial U_{ij}}{\partial x_k}+\frac{\partial U_{kj}}{\partial x_i}\right)n_k\right].$$

For the below-mentioned Kelvin’s fundamental solution:

$$U{\left(x,y\right)}_{ij}=\frac{1}{16\pi \mu (1-\nu )}\left[4\left(1-\nu \right)\frac{{\delta }_{ij}}{R}-\frac{{\partial }^{2}R}{\partial x_i\partial x_j}\right],$$

where $R=\sqrt{{\left(x_k-y_k\right)}^2}$ is the distance between the field x and integration y points; the matrix J^S has the form:

$$J_{ij}^S(x,y)=\frac{1}{8\pi (1-\nu )}[\frac{4\nu (1-\nu )}{1-2\nu }\frac{\partial }{\partial x_j}\frac{1}{R}{n}_i-\frac{\nu }{1-2\nu }\frac{{\partial }^{3}R}{\partial x_k^2\partial x_j}{n}_i+2(1-\nu )\left(\frac{\partial }{\partial x_k}\frac{1}{R}{\delta }_{ij}{n}_k+\frac{\partial }{\partial x_i}\frac{1}{R}{n}_j\right)-\frac{{\partial }^{3}R}{\partial x_k\partial x_i\partial x_j}{n}_k].$$

The hypersingular kernel is J^H(x) = T_nU^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/μ 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.2 Seeding flaws.

The seismic subroutine SeisInput, complementing the basic code, seeds rectangular flaws with centers randomly distributed in a parallelepiped with sizes three- to five-fold 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 π/2), strike angle (from 0-2π), rotation angle in the dip plane (for a rectangular flaw changing from 0 to π) 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 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.

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.2 Geometry 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.

4.2.3 Mining 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 σ₁₁ = −21.75 MPa. Since the mined-out area is large, the pressure of the roof on the floor is taken into account.

4.2.4 Initial 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 δ = 0.0156 corresponds to m = 4.

4.2.6 Comparison 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-activated flaws), the reduction of the time expense is much less. The resulting reduction is illustrated in Table II.

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.

4.3 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₁, x₂ and x₃ are directed, respectively, upward, to the east and north. These are the directions of principal stresses: σ₁₁ = −7.5 MPa, σ₂₂ = −5.0 MPa and σ₃₃ = −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₂. Thus, its plane is orthogonal to the direction of the minimal rock pressure 4 MPa, acting along the axis x₃. The propagation is caused by the pressure p of a fracturing fluid pumped at the center O with the pumping rate q₀ = 2 m³/min. The net pressure Δp = p + σ₃₃ is positive. The fluid is Newtonian with the dynamic viscosity μ = 3×10⁻⁸ MPa·s. The leak-off into rock formation follows the Carter dependence with the fluid loss coefficient C = 1.5⋅×10⁻⁴ m/√s. The changes in time t of the net pressure Δp, the fracture length L and the fracture width w are described by the approximate solution:

(7) $$\Delta p_D\left(t_D,x_{\text{2}}\right)=\frac{\text{160}}{\text{3\hspace{0.17em}π}}{w}_D\frac{h_D}{a_D}\frac{\text{2}}{\pi }\left[\text{1}+\text{1}\text{.515 exp}\left(-\text{1}\text{.15}{r}_D\right)\right]{\left(\text{1}-\frac{x_{\text{2}}}{L}\right)}^{\text{0}\text{.25}},L_D=0.48t_D^{0.552}$$

$$w_D=0.2\pi {\left[\frac{a_D/h{}_D}{1+1.515\hspace{0.17em}\exp \left(\mathrm{-}1.15r_D\right)}{L}_D\right]}^{0.25}$$

The dimensionless pressure Δp_D, length L_D, time t_D and opening w_D are defined as Δp_D = Δp/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:

$$p_n=\frac{3\pi }{32}{\left[\frac{2\mu G^2{q}_0^2}{{(1-\nu )}^2{C}^2{h}_0^4}\right]}^{1/3},L_n=\pi {\left[\frac{(1-\nu )\mu q_0^5}{256G{C}^8{h}_0^4}\right]}^{1/3},t_n={\pi }^2{\left[\frac{(1-\nu )\mu q_0^2}{32G{C}^5{h}_0}\right]}^{2/3},w_n={\pi }^2{\left[\frac{16(1-\nu )\mu q_0^2}{G{C}^2{h}_0}\right]}^{1/3}$$

Besides, h_D/a_D = h/a; r_D = a/b, where a and b are, respectively, the small and large sides of the rectangular HF: a = L, b = h when L < h (at the initial stage of the fracture propagation), whereas a = h, b = L when the fracture length L becomes greater than its height (L > h).

The basic code of the H-BEM uses the current net pressure and the current fracture length, defined by equation (7) at nodes of the mesh on the rectangular HF, at time steps of its propagation. The code evaluates stresses induced by the initial fracture and those arising due to the changes in the geometry and in the net pressure. The total number of elements for zero step geometry equals 10 (the total number of unknowns is 30). A ten-meter time step is approximated by ten rectangular elements, which update the geometry. The total number of steps is 10; thus the final number of unknowns is 300.

4.3.2 Input data on seeded flaws.

The volume with seeded flaws is taken as a parallelepiped with the size X₁ = 5h = 250 m in the vertical direction (x₁), X₂ = 8h = 400 m in the eastern direction (x₂) and X₃ = 2h = 100 m in the northern direction (x₃). Such a region of the total volume 10⁷ m³ is large enough to include flaws for simulating seismicity, induced by a HF of the height 50 m and the length 100 m, which appears after ten time steps of 10-m length each. The initial number of square flaws, seeded randomly in the volume, was 75,000. After exclusion of 6,493 of them by subroutine ExclusInSitu, the total number of flaws used in calculations is N = 68,507. The average size of flaws is d_av = 1 m.

4.3.3 Tree and influence lists.

For the number N = 68,507 of flaws, the hierarchical tree built includes six levels.

(K = 6). The total number of objects in the structure is 39,176, the last sixth level contains 33,280 objects (cubes) of the size 6.25 m. Approximately, two to three flaws are located in each cube of the last level. The accepted tolerance δ = 0.0156 corresponds to m = 4. Table III presents parameters used for simulation.

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.4 Comparison 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.

5. 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³) 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.

Figures

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

Figure 3.

The plan view of the seam Fifth in Severnaya Mine of the Vorkuta coal district used for approximation of geometry

Table I.

List of parameters used in simulation of seismicity in Severnaya Mine

Quantity	Value
Input data on mining conditions
Mining depth	H = 870 m
Elasticity modulus of rock	E = 37 GPa
Poisson’s ratio	ν = 0.25
Principal vertical stresses	σ11 = −21.75 MPa
Input data on seeded flaws
Initial cohesion	C0τ = 2.5 MPa
Average size of flaw	Dav = 2.56 m
Average distance	Lav = 7.9
Density of flaw	ξ = 0.32

Table II.

Number of simulated events and computational time without and with acceleration

	The zero step		Typical (sixth) step
Subroutine	Time of calculation (s)	No. of events	Time of calculation (s)	No. of events
SeisTrac	2013	10,168	105	631
FastSeisTrac	19	10,115	1	625

Table III.

List of parameters used in simulation of seismicity accompanying fracture propagation

Quantity	Value
Input data for the basic code
Depth of layer	H = 500 m
Elasticity modulus of rock	E = 20 GPa
Poisson’s ratio	ν = 0.25
Principal stresses	σ11 = −7.5 MPa
	σ22 = −5.0 MPa
	σ33 = −4.0 MPa
Input data on fluid
Dynamic viscosity	µ = 3×10−4MPa·s
Fluid loss coefficient	C = 1.5×10−4 m/√s
Constant inflow rate	q0 = 2 m/min
Input data on seeded flaws
Initial cohesion	C0τ = 0.1 MPa
Average size of flaw	dav = 1.0 m
Average distance	Lav = 5.11 m
Density of flaws	ξ = 0.20
Number of excluded flaws	6,493
Input data on fracture geometry
Fracture height	H = 50 m
Fracture total length	L = 100 m
Number of steps	10