Simulation Possibilites of the Post-Mining Goafs Impact on the Deformations Induced by Next Underground Mining Operations with Use of the Cellular Automata Method

Underground mining operations in the area of a rock mass affected by previous exploitation may cause additional deformations to appear on the surface. The size of these deformations can be significant, and their character is often non-linear. The nature of these deformations cannot be justified solely by the impact of current mining operations. At the same time, the predictive method of S. Knothe, widely used in Poland, does not explicitly include these types of phenomena. In the area of intensive and long-term mining exploitation, such as the Upper Silesian Coal Basin, the practical possibility of simulating this occurrence may be helpful in the planning of new mining exploitation under construction objects. Today we are usually limited to numerical modelling methods like finite difference method (FDM). This one base on the principle of mechanical similarity. The theoretical usefulness of method (and its similar) has already been proven many times. The main impediment to its practical application is the lack of recognition of the rock mass in terms of its mechanical properties. The presented method is a new approach to the possibility of modelling the subject phenomenon. The method has not been used in practical forecasting mining area deformation caused by underground deposits mining. It’s characterized by a huge potential for further development.


Introduction
Underground mining exploitation in the area of the Upper Silesian Coal Basin has been carried out continuously for several hundred years. In the first period mining works were carried out at small depths. The development of mining over the decades has contributed to the increase in the average depth of exploitation. This fact caused that currently mining operations take place in conditions of the rock mass strongly affected by previously extracted mining works (Kowalski et al. 2010;Strzałkowski 2010). This favours the so-called activating old goafs, and thus leads to the resumption of exploitation proceeds after the period of their earlier stabilisation (Głowacki et al. 2013;Pilecki 2011). Activation of old goafs has a large impact on the course of the rock mass deformation process. Observed influences of conducted mining operations show significant differences in relation to forecasts that do not take into account the presence of old mining workings. An example may be the noticeable influence of reactivation of old goaf, documented by measurements, in mining exploitation conditions at coal mine ''Marcel'' on the ''W'' measurement line (Fig. 1).
What draws attention here is the much higher value of subsidence measured than it would result from reprognosis, on the section between points 48-86. The additional impacts observed in this section are too large and too uneven so that they can only arise from the impact of the analysed exploitation. These influences are the result of activation of the previously extracted deposits and unfinished inflows of previous exploitation. In the past, before the measurements were done, an intensive exploitation of quite shallow longwalls was done under the analysed section of the measurement line. Deformations measured on this line perfectly illustrate the impact of mining exploitation in the conditions of a rock mass disturbed by earlier mining operations.
There are more examples that could be shown here (Niedojadło and Palka 2003;Niedojadło et al. 2006;Mielimąka 2009). Unfortunately, it is impossible to include and describe all of them here, even briefly. What makes them particularly special is the increase in observed subsidence on the land surface disproportionate to the thickness of the exploited deck. Usually, the phenomenon of adding up of mining influences is non-linear. It is noteworthy that each of examples is different and there is no single way to influence the goaf impact on the next distribution of rock mass deformation. Each time, the available information on the structure of the rock mass, its non-uniformity and the multiplicity of previous mining operations should be taken into account, including the time factor. The mutual positioning of the exploitation parcels and the direction of their exploitation is not without significance.
The finite difference method (FDM) gives considerable possibilities for modelling the described phenomenon (Wesołowski 2013). It has already been shown that it can be used to simulate the impact of old goafs on the actual displacement distribution process as well as a number of factors that cause non-linear addition of rock mass displacements caused by underground mining. However, a serious impediment to the practical application of the method is the lack of sufficient recognition of the rock mass mechanical properties. In addition, it happens that a rock sample subjected to strength laboratory tests under real conditions behaves differently.
A completely different approach to the problem of simulating the impact of heterogeneity of the rock mass structure and the possibility of simulating the process of its displacement distribution in a non-linear way is the application of cellular automata theory (Sikora 2019). The method base on the principle of geometric similarity and allows to describe in an extremely simple way the influence of the heterogeneity of the mining rock on the distribution of its displacements caused by underground mining operations. In the scope of mining area protection and broadly defined rock mechanics, the method has not been used in practice so far. Also, there were used Author's software solutions due to the lack of existing ones that can be used. The article attempts to determine possibilities of simulation of the old goafs impact on the size and distribution of deformation of the rock mass in case of new mining operations with the use of cellular automaton method. The advantage of the cellular automata method in the context of the subject of research is its characteristics, which is based on the principle of geometric similarity. The parameters of the method correlate with the parameters of the S. Knothe method, which allows the selection of parameters, among others from the results of geodetic observations. At the same time, the method allows for taking into account in a natural way the heterogeneity of the rock mass structure and the properties of nonlinear summation of mining influences.
The paper presents the assumptions of the method and shows its operation on a representative abstract example. First numerical calculations were made using the proven FDM method. The reference calculation model was based on the work of Kołodziejczyk and Wesołowski (2010) describing how to select deformation parameters in relation to typical values of parameters of Knothe's theory. Then, by fittingusing the least squares method-a simulation was performed using a finite cellular automaton. The comparison showed compliance in qualitative and quantitative terms. Classic Knothe theory does not directly capture the influence of many factors, among others associated with the heterogeneity of the rock mass structure and the phenomenon of non-linear summation of influences caused by the impact of earlier operation. Obtaining similar results by the geometric-integral method may be difficult.

Cellular Automata Method
The theory of cellular automata (CA) is completely different from the previously used methods of forecasting rock mass continuous deformations caused by underground mining. This method is a new approach in this area. The basic considered computational model is the finite state automaton (Sikora 2010). There are also a number of other concepts for using cellular theory and automata in this issue (Niemiec 1985;Saavedra-Rosas and Jarosz 2012). It is an abstract, mathematical and iterative model of the behaviour of a dynamic system, the structures of which are described by the network of cells and their states, transitions and rules of these transitions. The essence of the method is the possibility to naturally take into account the variables of the properties of the medium, which in effect allows to describe and analyse the impact of heterogeneity of the rock mass and nonlinear properties of the rock mass on the surface subsidence. The method adopts the principle of geometric similarity. Originally the characteristics of the method is more consistent with the characteristics of the loose centre of Litwiniszyn (1954). Its parameters are simple to physical interpretation. They can be determined based on the results of observations of the effects of previous underground exploitation. The method does not require special mathematical knowledge to apply it and is extremely effective in comparison to other numerical methods.
The precursor of the idea of application the method in the field of rock mass mechanics and protection of mining areas is Niemiec (1985). In turn, Sikora (2010Sikora ( , 2016Sikora ( , 2019) developed a mathematical characteristic allowing for its practical application. The works of Białek and Sikora (2012) and Mielimąka and Sikora (2017) confirmed the great potential of the method development. The theoretical assumptions have been pre-verified on examples of actual exploitation.

The Computational Model as a Cellular Automaton
The computational model treated as a finite cellular automaton is constructed by discretizing a fragment of a rock mass in the form of a regular cellular grid (Sikora 2010). By definition, all cells have the same shape and adhere closely to each other. The basic parameters of the model are the dimensions of a single cell reproduced in reality. In the case of a twodimensional model, they are marked as S k -width and W k -height (Fig. 2). The model simulates the process of filling the postmining void with rock rubble. Model assumptions can be greatly simplified by assuming that only the force of gravity affects the fulfilling of emptiness (Sikora 2016). The model works by automatically executing a specific algorithm (the so-called transition function) on each cell separately (Sikora 2019). In the process of ''outflow'' of the volume of the post-extracted emptiness assigned to a given cell, in the two-dimensional model, the most often considered is a three-celled neighbourhood (space around a given cell within which data is exchanged). In the distribution the parameter of the main transition P is distinguished. Its value determines the part of the base cell ''subsidence'' that will be passed to the cell that is directly above it. The rest is usually passed in equal proportions to the remaining cells from the neighbourhood (Fig. 3).
In a programming environment, the task is carried out by a programming loop (Sikora 2019). As a part of the loop operation, each cell-row by row and column by column-is subjected to the same transition function that determines the range and proportions of the subsidence transition from a given base cell. An example showing the temporary position of the transition in the automaton grid is shown in Fig. 4.
The simulation consists in a separate deformation distribution for each elementary cell initiating the initial post-operation emptiness (cell in row No. 1 and  (Sikora 2016) column No. 4 in the Fig. 4). The distribution simulation from subsequent cells is carried out in accordance with the direction of progress of mining works. The results are summed up in the same grid or separated one called collection grid. The same transition function used each time causes the change of the state of specific cells in relation to the cell evolved at a given moment.
Numerous observations of the effects of the method led to develop the initial characteristics allowing its practical application (Eq. 1) (Sikora 2010).
where: I max -maximal inclination in case of a full subsidence trough, ag-maximal subsidence in case of a full subsidence trough (a-mining coefficient of S. Knothe's theory, g-exploitation thickness), Hdepth of mining exploitation, A-adjustment parameter due to the value of the main transition parameter P.
In case of the geometrical-integral Knothe's theory (1984), a number of geological properties of the rock mass defining its ability to create subsidence trough with a given slope are captured in the form of a single parameter tgb. In the cellular automata method, the parameter of the maximum slope a I is defined analogously to the tgb parameter (Eq. 2) (Sikora 2016).
In contrast to the tgb parameter, the a I parameter is a variable depending on the depth of mining exploitation H. This confirms the results of Rogusz's research (1977) on the variability of the tgb parameter depending on the depth of mining operations.
The maximum slope parameter a I can be determined from the criterion of matching the theoretical subsidence trough to the trough resulting from geodetic measurements.
The distribution function can be non-linear. There are different ways to solve this problem. One of them, modelled on a similar solution by Litwiniszyn (1954), is the dependence of the transition function on the slope value in a given cell. Simulation of the distribution of deformations for the entire seam, carried out cell by cell, causes the deformation state of the collective table to change constantly. In practice, the effects of non-linear summation of inflows are obtained by making the value of the main transition P dependent on the current value of the slope I in a given base cell at row m and column n (Białek and Sikora 2012) according to the dependence.
where: P p -value of the main transition for the undisturbed distribution function, z d -parameter defining the acceptable variation of the main transition P p , a d -parameter of a nonlinear distribution function depending on the inclination value I already saved in a given cell about the index m,n; a d 2 R þ . The selection of parameters delinealizing the distribution is most often based on the percentage assessment of the subsidence trough asymmetry and iterative selection of a d and z d parameters.
Also, it can be assumed that this distribution process will occur with a certain loss. The assumption correlates with observations that indicate that in the case of the full trough maximum subsidence is smaller than the longwall's operational height. The degree of delinearization and the process of accumulation of part of the subsidence in the computational model can be additionally dependent on distance of the elementary cell imaging the extracted volume (the progress of the exploitation front) from the initial cell (start-up zone). In this way can be simulated volume changes and in particular the observed (in reality) phenomenon of gradual shallow of the subsidence trough bottom. With regard to the structure of the computational model and the characteristics of the method, it can be stated that post-mining goafs can be mapped in at least two ways. The first consists in performing numerical calculations for subsequent longwalls-one by onetaking into account the impact of previously registered displacements on the distribution of current deformations. The simulation involves evaluating all cells that make up the entire longwall one by one. The results are added up in a summary table. During iterative distribution of subsidence (equivalent of post-mined void) from a single cell to cells in a cellular neighbourhood, some of the subsidence previously recorded, which are derived from distribution simulations for the preceding longwalls, will be additionally added. Some difficulty in using this solution is to determine the part of the subsidence from a given cell (in summary table) that will affect the current distribution. The second approach is mapping previously extracted longwall as a separate one, assuming that its thickness will be a certain part of the original thickness. In this way, the assumption is made that the inflows resulting from its reactivation will be a part of the primary inflows. The issue has been described, among others by Szpetkowski and Pytlarz (1967) who, as a result of observation of the effects of exploitation, estimated that this impact is about 12% of the original impact. A lot has changed since then-including exploitation is carried out at a greater depth, however, the author's model experiments show that in the case of the cellular automata method, a value of 10-12% is appropriate in the absence of additional information. This solution is much less complicated in terms of implementation in the computational model and due to the practical possibility of using the method. Especially when mapping more longwalls in one automaton grid. The main advantage of the solution is the lack of the necessity of inputting new parameters difficult (or impossible) to be determined under real operating conditions.

Computational Model
In order to verify the theoretical assumptions regarding the possibility of taking into account the impact of post-mining goafs using the cellular automata method, first a representative computational model was built. The example concerns the exploitation of two abstract longwalls. The first stage of calculation is to determine the deformations after extracting the first one. In the next stage, the possibility of taking into account the impact of deformation caused by mining out the first longwall on deformations caused by extraction the second longwall were demonstrated.
The first was made simulation using the FDM method, providing a reference for the cellular automata method.

Computational Model Structure
To carry out simulated mining exploitation and determine its effect on surface deformations, a model was constructed in a flat state of deformation with dimensions of 2400 m 9 1100 m. At the depth of 600 m (depth of the seam floor), the 2-m-thick deck designated for exploitation was modelled. The assumed field length in this seam will be 1050 m.
At a distance of 150 m below deck No. 1, a 3-m height longwall No. 2 is mapped. The assumed field length in this seam is 700 m, while the edge of exploitation in seam 2 has been shifted by a distance of 150 m relative to the operational edge in the seam 1. Mutual location of exploitation fields in both seams will also allow to determine the impact of old goafs and influence of the edges of these goafs on the deformation of the model's surface. The geometric scheme of the model is shown in Fig. 5.

Finite Difference Method
Numerical modelling, carried out for the purpose of this work, was done with a use of the ITASca Consulting Group FLAC programme (Fast Lagrangian Analysis of Continua) based on the finite difference method (FDM) (FLAC 1992).
The FLAC program is used to build numerical rock mass models and to simulate the behaviour of ground and rock centres that experience plastic flow or brittle fracture after reaching the point of ductility or strength limit. This program is especially recommended for solving rock engineering problems. It is used both to assess the behaviour of the rock mass in the area of mining excavations, and to simulate phenomena occurring in significantly larger areas. The FLAC program enables modelling rock mass discontinuities, simulating the behaviour of rock mass and structures under dynamic constraints (Wesołowski 2013). The FLAC program also allows for mapping of plastic changes in the rock mass during the simulation of mining exploitation (FLAC 1992; Kołodziejczyk and Wesołowski 2010).

Numeric Model in the FLAC Program
The structural model presented in Fig. 5 is the basis for the construction of the numerical model in the FLAC finite difference program. To carry out numerical calculations, the model's space was divided by a rectangular net with dimensions of 10 by10 m (floor and roof of the seam) and in the case of seams 10 m 9 2 m (seam 1) and 10 m 9 3 m (seam 2).
Building the network, it was assumed that nodal points, located on the vertical side edges of the disc, can move freely in the vertical direction, and in the horizontal direction their displacements are equal to zero. The joints located on the basis of the model disk can move freely in the horizontal direction. The vertical value of displacements of these points was determined as zero. Other nodal points belonging to the model are free to move in any direction of the X-Z plane.
When determining the boundary conditions, it was assumed that the value of the primary stresses in the rock mass comes only from the gravitational forces. In the case of a non-tectonic rock mass, this assumption is sufficient to determine the initial conditions of simulated exploitation (Tajduś 2010).
As was shown (Wesołowski 2013), taking into account the layered structure of the rock mass model enables obtaining the most accurate determination of most deformation indices. Taking into account the layered structure of the rock mass, it was assumed that the ceiling and the floor of seams will be described with homogeneous, alternating layers of sandstone and shale with the same thickness of 40 m. The last surface layer with a thickness of 20 m will be a cover layer.
In terms of mathematics, it was assumed that all layers building the rock mass model were described by the ubiquitous joint centre. This model is an anisotropic plastic centre containing the areas of weakness of a particular orientation. In this model, Coulomb-Mohr's endurance (plasticity) condition was implemented (Labuz and Zang 2012). The plasticization can occur both within the areas of weakness and rock mass. Strength and deformation parameters of the layers were adopted on the basis of literature (Prusek and Bock 2008;Tajduś 2013).
Defining the values of the parameters of the weakening planes, the case described in the Sainsbury paper was used (Sainsbury et al. 2008). The range of variability of material parameters of rock layers of the ubiquitous joint model adopted for calculations are presented in Table 1 (Kołodziejczyk and Wesołowski 2010).
The simulation of exploitation will consist in cyclical removal of individual finite element mesh zones in the order corresponding to the extraction of longwalls. At the same time, contact elements (separation plane) are inserted between the roof and floor layer, which prevents mutual penetration of the ceiling and floor of the deck. This way of mining exploitation simulation eliminates the necessity of introducing additional parameters attributed to infarct zones. The simulation of the rock mass displacement corresponds to the quasi-static model, in which the effects caused by the operation (removal) of the next element of the seam are immediately revealed. There are no sticky elements and no time concept. Due to the plastic properties of the centre, the order of removal of individual mesh zones is important. The adopted model assumptions can be considered average in terms of the method in relation to the observed typical mining influences on the land surface (Kołodziejczyk and Wesołowski 2010).

Building a Model as a Cellular Automaton
In this case, for the numerical calculations, a model was built in which were assumed the values of the basic parameters correlating with the typical values in the case of geometric and integral methods (in the range of Polish mining conditions). The value of the maximum slope parameter a I = 2.0 was assumed. Based on the mathematical characterization of the method (Eq. 1 and 2) the height of the cell W k and the value of the main transition P were selected iteratively, assuming constant cell width S k = 10 m. The adopted parameter values can be considered average.
Non-linear effects were implemented based on the results of a computer simulation carried out for the ubiquitous joint model. Parameters of the delinearizing function (Eq. 3) and the degree of accumulation of subsidence in the grid of the automaton were selected iteratively to match the asymmetry effect obtained with FDM method.
Post-mining goafs were included in this case as a separate deck. It was assumed that the inflows resulting from its reactivation will be a part of the original influences (Szpetkowski and Pytlarz 1967). Finally, it was assumed that the thickness of the longwall No.1 will be 12% of the initial one, assuming the mining operation with the roof collapse method. The simulation was carried out assuming unchanged parameters in relation to simulations of vertical displacements caused by the original exploitation in the seam No. 1. Nonlinear transition function with the The angle of inclination of planes a 1. Mining operations simulated by the Finite Difference Method under extracted field in seam 1 results in the creation of a subsidence trough, for which the maximum subsidence value is s max-= 3003 mm. Assuming such a value of maximum subsidence and a 3.0 m operation thickness, the determined value of the exploitation coefficient should be a = 1,001. With regard to the comparative model, the determined value of the maximum subsidence is higher by approx. 16% to the maximum subsidence obtained in the conditions of the rock mass not damaged by earlier mining exploitation (Fig. 7). The reason for such a significant increasement of subsidence is the volume changes of rocks and goafs lying above the seam 2. 2. There is also a noticeable increase in the range of the resulting subsidence on the surface of the model. Additional subsidence covers the whole goaf area within longwall No. 1. With reference to the comparative model, the range of (additional) subsidence increased by approx. 380 m. In this part of the trough, the influence of the edges of the finished exploitation in seam 1 is clearly marked (ordinate x = 1960 m). In case of starting edges, the increasement of subsidence (in order to the comparative model) was small and doesn't exceed 40 mm. 3. The results obtained by numerical calculations using the cellular automata method are very similar (Fig. 7). As a result of numerical calculations, a subsidence trough was created, in which the maximum subsidence amounted to s max-= 2992 mm. In the area of the starting edge, the nature of the deformation in the case of results obtained with both methods is very similar. Certain differences in the quantitative range can be seen in the region of the other edge, in particular over the operational edge in seam 2 and in the middle of the goaf in deck 1. In this case, the differences in subsidence reaches around 25%. Nevertheless, it should be recognized that the similarity is high. The coefficient of determination determining the degree of matching of both subsidence profiles is R 2 = 0.98.

Summary
The article presents the results of numerical analysis of the impact of so-called goaf and rib area on current underground mining exploitation. Using the FLAC3D finite difference program and the cellular automata method (and proprietary software CA3D created for the purposes of conducting research), a simplified model was created which was an idealized example of mining exploitation under goaf zone. The adopted model assumptions can be considered average in terms of the method in relation to the observed typical mining influences on the land surface (in conditions of Polish mining operations). In order to determine the impact of reactivation effect on the model's surface deformation, an additional comparative simulation was performed, where the post-mined zone was omitted. Based on the results of numerical calculations, the following conclusions were made: Considering the above conclusions, it should be stated that the main goal of the work has been achieved. The article demonstrates that using the cellular automata method it is possible to simulate the impact of post-mining goafs on the distribution of rock mass subsidences caused by underground mining operations. The possibilities relate to the simulation of the effect of a disproportionate increase of subsidences on the land surface in relation to the thickness of the extracted longwall and the volume changes observed in reality revealing in the significant asymmetry of the subsidence trough. These effects are not captured by commonly used for forecasting geometric-integral methods.
At the same time, attention should be paid to the simplicity of the method. The selection of parameter values can be done basing on results of geodetic observations or on the basis of average values for a given mining region, which closely correlate with the parameters of S. Knothe's theory. In the manner presented in the article, the cellular automata method can be treated as a supplement to elements that are not explicitly included in geometrical-integral methods.
The method assumptions presented in the article can be used in the future to assess specific cases of actual underground mining operations in regions where intensive mining has previously been carried out.
Funding Not applicable.

Compliance with Ethical Standards
Conflict of interest The author declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.