Agnieszk Mlinowsk,Ryszrd Hejmnowski,Huyng Di

a AGH University of Science and Technology,Al.Mickiewicza 30,30059 Krakow,Poland

b China University of Mining and Technology,Beijing,D11 Xueyuan Road,100083 Beijing,China

Keywords:Mining operation Modeling of surface deformations Stochastic models Probability integration method Influence function Salt hard coal mining copper ore mining

ABSTRACT Mathematical modeling of surface deformations caused by underground mining operation is commonly carried out with use of empirical,numerical or stochastic models.One of the most frequently applied model for prediction of ground deformation in many countries is Knothe model.The model developed by Knothe belongs to the stochastic methods and is based on the influence function.In China a prediction method named Probability Integration Method (PIF)was established by Liu Baochen and Liao Guohua based on the stochastic medium theory.Modified version of that model allows to predict ground movements caused by mining operation in extremely complex technical and geological conditions.That model is commonly applied for coal,metal ore and salt deposits.The article presents several modifications of the mathematical model used in China and Poland.This model is very widespread in the world,therefore the generalizations proposed in the article can be implemented for the purposes of prediction surface deformations for various types of deposits in many countries.The presented generalizations were then tested on specific examples of coal mining,copper ore mining and rock salt deposit.The obtained results indicate high efficiency of methods based on the influence function in complex geological and mining conditions.

1.Introduction

Mining-induced ground movements can be predicted with several geomechanical and stochastic models.One of the most frequently used models applied across the world is the influence function method,i.e.geometric integral methods.These methods support effective prediction of mining-induced ground movements without involving much labor and calculation power.These methods stemmed from the first calculation model proposed in Germany at the beginning of the 20th century.At that time the calculation algorithms were used for predicting vertical strains to any point of the surface,which were a function elementary area of the deposit extracted in a given region.Initially these methods were very simple,though with time they were generalized to integral calculation models [1,2].In the 1950s and early 1960s these geometrical-integral methods were transformed into a form of compact Ehrhardt-Sauer method,which was based on Knothe-Budryk model [3-5].The Knothe-Budryk model was developed in the successive decades in Poland into numerous applications for different types of mining operations like coal mining,copper ore mining,and salt mining[6-14].Underground coal mining condition in central and eastern Europe are similar due to the fact that Knothe theory was also implemented for subsidence modeling in Germany and Slovakia[15-20].The mining conditions in western Europe are more complicated due to the steep seam mining,but also for that area influence function methods were adjusted [21,22].The Knothe theory was also adjusted to the mining conditions in China,Vietnam,USA,Australia,Turkey,India,and Iran[23-35].

The theory of predicting ground movements applying the basic Knothe model and influence functions has remained unchanged since its beginning in the 1950s,although the future adaptations to various types of mineral deposits are an interesting solution.They required adaptations and generalizations which considerably exceeded the original,without minimizing its author’s merits.These new elements were discussed at greater length in the paper as they grow in significance now.

2.Original Knothe-Budryk calculation model

2.1.Influence function and modeling components of vertical and horizontal strains

As mentioned above,the research was based on a stochastic model relying on Knothe influence function[4,5].According to this model,the exploitation of an elementary volume of deposit (dV)and the horizontal plane generates an elementary subsidence at point A.This subsidence can be described as:

where dsAis the elementary subsidence at point A;f the influence function;dP the elementary area of the chamber;and x and y the coordinates of point A.

Employing the principle of superposition,it has been assumed that vertical displacement/subsidence at point A is a sum of elementary subsidence coming from all elementary volumes in the exploitation field area:

where sAis the vertical subsidence at point A.

In Knothe model,the following generalized influence function for a flat state of deformation is assumed as following:

where smaxis the maximum/final subsidence;and h the parameter of the influence range.smaxcan be calculated by the following equation:

where α is the extraction ratio depending on the way the post

mining void has been filled out;and g the thickness of exploitation.Following the parametrization of the influence function,i.e.assuming the limited range of influence of the exploitation and introduction of the parameter R,Eq.(3)can also have this form:

where R is the radius of the influence range,also called the range of the major influence.R is connected with the strength characteristic of the rock mass through Eq.(6):

where tanβ is the parameter corresponding to strength properties of rock mass;and H the mining depth.

Making use of the assumption that the horizontal displacements at point A are proportional to the derivative of vertical deformations,Budryk (after Avershyn)implemented Eqs.(7)and(8),by which the horizontal displacements and horizontal strain could be predicted at any point (taking point A as an example)on the surface,into the Knothe model [36]:

where uAand vAare the horizontal dispalcements at point A;εAthe horizontal strain at point A;andthe coefficient of horizontal displacements.

The maximum strain at the calculation point is predicted for a planned mining to be realized in a given time.The strain,calculated for a given moment,can be predicted from Eq.(9):

There are usually two values of extreme perpendicular-oriented strains which are determined at the given point A (Fig.1).The major angle,i.e.left angle,is measured in reference to the positive semi axis of local coordinates system,at which the extreme strain occurs.The angle of the maximum strain is calculated based on Eq.(10):

where φεmaxis the left angle of the maximum horizontal strain εmax.

The above model is congruent with Litwiniszyn stochastic media theory,especially the assumed non-compressibility and compaction of the medium and superposition of influence [37].In practice these assumptions turn out to be purely theoretical and do not correspond to the real geological condition in the rock mass.The postulated stochastic characteristics of the medium was not considered in numerous modifications of the Knothe model.

2.2.Modeling of deformations in a function of time

In the 1950s Knothe published a differential equation on the basis of which a dependence could be derived for modeling the subsidence in time,addressing the issue of manifestation of mining impacts on the surface in a function of time [4,5]:

where c is the coefficient of time depending on local geological conditions (0.001 ＜c ＜30);sFthe final potential subsidence at a point due to the ongoing extractions;and s(t)the momentary actual subsidence at a point due to the ongoing extractions.

The following Eq.(12),shown in Fig.2 as a solution of Eq.(11)with the boundary conditions as s(t=0)=0 and also instantaneous depletion of the deposit,allows determining the subsidence at any time after stopping or terminating the extraction:

As for the constantly advancing extraction,the solution of Eq.(11)is slightly more complex:

The time of subsidence that the mining impacts on after the activity is terminated may range from a few to hundreds of years,depending on the c value.For hard coal extraction the time when the impacts are visible ranges from 0.2 to 2.5 years whereas for fluidal deposits and salt domes it may be tens or hundreds of years.

Fig.1.Modeled distribution of horizontal strains in the given point A.

Fig.2.Standardized subsidence over time with different time coefficients calculated by Eq.(12).

3.Adaptation of Knothe's theory for modeling continuous deformations in various types of deposits

Depending on the type of deposit,its underground exploitation evokes surface deformations varying in time and space.This is mainly connected with the geological building of the deposit,its geomechanical properties,and mainly the mode of extraction and liquidation of post-mining voids.Depending on the dynamics with which the impacts appear and their time and space course,the parametrization of the influence function seems to be a key element.The parameters of the model can be evaluated when the displacements are known,e.g.measured with geodesic surveying methods over the terminated extraction conducted in similar geological and mining conditions.Then the Knothe’s theory model can be adjusted to local conditions.

3.1.Calculation algorithms

The original Knothe model was based on the influence function of a parametrized Gauss function.Owing to the integral form of vertical strains model (Eq.(2)),Eq.(5)was no longer efficient.Therefore,calculations were performed for schematized rectangular extraction lots and tabularized integral solutions.With time,as the informatics infrastructure developed,numerical algorithms were available,and some of them are worthy of special attention.In 1989 Drze˛z´la suggested a polar integration algorithm(along the radius vector and angle),which facilitated the modeling for extraction in almost every shape and schematized the fragments of a circle section [6].Such numerical integration was also implemented in computer programs by Piwowarski and Zych [10].Białek proposed other modifications of the original Knothe model by introduce of additional parameters connected with propagation of influence to the basic influence function [9].And Hejmanowski introduced a numerical algorithm based on rasterizaion of the modelled deposit [7].

3.1.1.Hejmanowski’s numerical algorithm

The principle of the algorithm is substituting the integral with summing of elementary extraction impacts evoked by elementary volume element (EVE):

where N is the number (quantity)of exploitation elements;i the number of a subsequent individual exploitation element;j the number of the (individual)calculation point;t the time difference between the exploitation date of every EVE and the date of the prediction;V(t)the volume of the post-exploitation elementary subsidence bowl at the time t,due to mining of the layer volume in the ith element of exploitation;di,jthe (horizontal)distance between the point calculation and the centre of the exploitation element;and fz(di,j,Hi)=the influence function.

The deposit is divided into an arbitrary number of such elements (Fig.3).EVE has a shape of square base cuboids,uniform for the whole modeled reservoir.The horizontal dimension(square)is adjusted to the distance of the reservoir from the calculation level,accounting for the accuracy.The EVE generates an elementary subsidence trough.At the same time,the EVE is a carrier of local information about the deposit and the way in which it is exploited.It may contain information about the local depth,local thickness,extraction coefficient,time of extraction,pressure (fluidal deposits),etc.

Based on this approach,the algorithm facilitates modeling for arbitrary exploitation lots and any spatial range of deposition.This model can be equally used both for mineral deposits having the form of deposit formations and for veins.The influence of a given fragment of reservoir on the rock mass displacements can be modelled by introducing the surface information of the planned mining area within the EVE as compared to the total EVE area (Fig.4).Numbers corresponding to the percentage of mined area of a given EVE are introduced in particular cells of the model.Thanks to this the volume of an arbitrary EVE,i.e.its influence on rock mass displacement,can be calculated as following:

where VEis the volume of EVE to be mined;E the ratio of surface of mining area to total surface in a given EVE;L the length of the edge of EVE;and g the local thickness of deposit layer to be mined.

The algorithm allows modeling the development of mining in time,even though slightly different functions of time and detailed solutions are applied for various types of deposits.The twoparameter Schober-Sroka function of time is a very good solution especially for deposits where the post-extraction voids tend to contract slightly slower [17].

Schober-Sroka function,whose properties were discussed in detail in relation with the cause of deformation,i.e.the convergence of working with the subsidence basin.The volume of elementary basin needed in the forecast of displacements (Eq.(8))is expressed with the following equation [8,17]:

where θ(t)is the function of time and is expressed as following:

where ξ is the convergence coefficient characterizing delaying properties of roof rocks in the workings area;and ϑ the coefficient characterizing delaying properties of the caprock in displacement propagation.

These two time parameters are the same as the coefficient of time (c)introduced by Knothe:

Fig.3.Division of the deposit into EVE and generation of elementary subsidence troughs.

Fig.4.Modeling of non-mined out pillar.

For the traditional hard coal and ore extraction,the coefficient ξ has the following values:(1)15 ＜ξ ＜50 for hard coal (goaf system);(2)0.2 ＜ξ ＜1.5 for copper ores (room and pillar);(3)0.4 ＜ξ ＜1.0 for ore mining;and(4)0.001 ＜ξ ＜0.1 for the salt mining and oil/natural gas deposits.

Over years this model was adjusted to the specific characteristics of manifestation of mining impacts on the copper ore,hard coal and even salt deposits.

3.1.2.Numerical algorithm applied in China by huayang dai

The probability integral method is used to calculate the surface movement and deformation caused by strip mining in China [23-25].Based on the following assumptions,the ground deformation can be simulated:

For the subsidence,

For the tilt,

For the curvature,

For the horizontal movement,

And for the horizontal deformation,

Table 1 Prediction parameters of probability integral method.

where Wcmis the maximum subsidence of surface at full mining;Ucmthe maximum horizontal movement of surface at full mining;r the main influence radius;θ0the mining influence propagation angle;and D the mining area(taking the inflection point offset into account).

3.2.Modeling continuous deformations for hard coal deposits in China

Based on the probability integral method,the ground movements and deformation caused by strip mining in China were calculated [23-25].In order to adjust the prediction model to the geological and mining condition,the local parameters were established (Table 1).

The prediction software is the Mining Subsidence Analysis System(MSAS)developed by China University of Mining and Technology (Beijing)[24,38,39].

3.2.1.Geological and mining conditions

The coal strata of the mining area No.8 of the Wutongzhuang coal mine in Fengfeng mining field belongs to Carboniferous and Permian [40].At present,coal seam #2 is designed and manufactured,with a dip angle of 3° to 18° (10° as the average),and with an average coal thickness of 3.3 m.The surface elevation is 150 to 257 m,the alluvium thickness is 56.3 to 138.7 m,and the submersible water level is 10 m below the surface.In this area,a coordinated mining technology is adopted for coal pillar mining in village W.By February 2017,the 803 and 807 working faces have been mined out(Fig.5).It is necessary to analyze the results of predicted surface movement and village mining impact on 805 and 811 working faces.The basic parameters of the working faces are given in Table 2.

3.2.2.Analysis of the prediction results

The presented result of the modeling relates to the ground movements caused by underground coal mining of longwall panels 805 and 811 (Fig.6).The following values were obtained:maximum subsidence of surface of 1981 mm,maximum tilt along south direction of 8.3 mm/m,maximum tensile deformation of 6.4 mm/m,and maximum tilt along east direction of 10.7 mm/m.

In the village area,the maximum subsidence was 1800 mm,maximum tilt along north and south direction was 7.5 mm/m and maximum tensile deformation value was 3.0 mm/m (Fig.6).

According to the assessment in the village area 45 houses will be damaged at grade Ⅱ,and 49 damaged at grade III,while others will be damaged at grade I.The distribution of predicted house damage grade is shown in Fig.6.

According to the prediction results,grade I damage means the simple maintenance,grade Ⅱand grade III mainly adopt maintenance measures.The most hazarded areas are classified as Ⅳcategory.In order to ensure the safety in build-up area,appropriate hazard managements need to be taken into consideration.So,the predicted ground movements can support planning how to strengthen the buildings and also that the monitoring policy should be.

Fig.5.Mining area No.8 of the Wutongzhuang coal mine (after [40]).

3.3.Modeling continuous deformations for copper ore deposits

Modeling vertical strains and deformations for copper ore deposits was carried out for a deposit in southwest Poland at a depth exceeding 950 m,the thickness of the seam oscillated from 3 to 5 m and tilt of 5°-8°.

The prediction was calculated with application of the MODEZ software.The calculation algorithm was based on division mining panels into EVE,which generates definite impacts.

The presented calculation model was implemented into MODEZ software.That IT system allows modeling ground movement and deformation at indicated periods of time.The calculations were performed for the next 50 years of exploitation,determining extreme horizontal deformations in time (Eq.(8))and maximum surface subsidence.

The calculation model is parametrized with properly selected or estimated parameters,which characterize geomechanical properties of the rock masses.During parametrization of the model for concession purposes attention should be paid to the variability of parameters in the mining areas.Moreover,it was necessary to take local geological conditions into account.The following parameters of modified Knothe model were used in this presented example:(1)tanβ-the tangent of radius of major influence range in Knothe’s theory,tanβ=1.35 ÷ 1.55;(2)α-the extraction coefficient,depending on the mining system.The following values of extraction coefficient were assumed:α=0.5 (exploitation of the deposit with roof deflection,the production from this part totaled to ca.80% of resources);α=0.2 (roof deflection and backfilling);α=0.1 (cutting out of the deposit into big stoops);(3)u-the parameter describing deviation of subsidence trough in view of the tilt of the rock mass,u=0.6;(4)ξ-the parameter of function of time describing delay of convergence of a selected element of extraction field,ξ=0.5-1.0/year;(5)η-the time parameter describing delayed properties of caprock,η=3.0/year.

With these values of time parameters above,the global coefficient of time(Eq.(19))varied in the following range:(1)c-the general parameter of time,c=0.4-0.75/year;(2)rs.-the radius of the influence range on the level of roof of the exploited deposit,rs.=50 m;(3)B-the coefficient of proportionality between vertical and horizontal deformations,B=0.32R.

The continuous deformations could be predicted based on adjusted model and estimated parameters.The subsidence caused by copper ore mining in Poland does not exceed 4 m.The time at which such impacts occur on the terrain surface is about 10 years.In the analyzed area the maximum of predicted surface subsidence is 2.75 m(Fig.7).Three predicted terrain categories estimated base on predicted horizontal strains are presented in Fig.7:(1)I terrain category-horizontal strains form 0.3 to 1.5 mm/m(marked in yellow);(2)Ⅱterrain category-horizontal strains form 1.51 to 3 mm/m (marked in orange);and (3)III terrain category-horizontal strains form 3.01 to 6 mm/m (marked in red).

3.4.Modeling continuous deformations for Wieliczka salt domes

Adapting the theoretical model for simulation of continuous deformation in the salt deposit conditions was a complex problem.The theoretical model was developed for a loose or brittle rock masses,not for the rheological one (as salt deposit).For the salt deposit the presented modification of the original model does not meet the basic assumption about requirements of continuity of rock masses.The first works on the adaptation of the model to the simulation of salt rock mass deformation dated back to the early 1980 s.Schober and Sroka proposed a two-parameter function of time (Eq.(18)),thanks to which the rate at which the salt caverns converge can be taken into account [17].This function was also used for rock mass deformation modelling presented in the paper.The solution firstly applied in Germany for salt caverns deformation modelling for Asse salt mine,was implemented for Polish salt mines.The adapted Knothe model for salt deposit was applied for rock mass modelling in historical,over 700 years old,Wieliczka salt mine.Since the 1990 s the mine has not been active,but the progressing convergence of over 2000 chambers generates a constantly increasing rock masses displacements.The Wieliczka town is in the center of these displacements therefore old workings need to be constantly protected,and on the other hand,deformations should be predicted and surface movements should be monitored.

When predicting salt rock masses deformation,the following parameters were applied:(1)tanβ=0.72-1.0;(2)α=1.0;for backfilled or protected chambers otherwise,α is determined on the basis of different solutions;(3)u=0.6;(4)ξ=(0.0007÷0.002)/year;and (5)η=3.0/year.

With thus assumed values of the parameter of time,the global coefficient of time(Eq.(19))was changed in the following way:(1)c=(0.0007-0.002)/year;(2)rs.=50 m;and (3)B=0.22R.

The subsidence of terrain in the next 100 years was assessed on the basis of the results of predictions for the given parameters and assumed geometry of underground salt chambers.The predicted maximum subsidence of terrain caused by convergence of underground workings is 1.4 m (Fig.8).The dynamics of subsidence occurrence will be very slow.The annual increase of subsidenceshould not exceed 7 cm in zones where the surface movements are the most significant.

Table 2 Basic information of the working faces.

Fig.6.Predicted modeling results of the ground movements (after [40]).

Fig.7.Predicted subsidence in next 50 years for underground copper ore mine in meters (after [12]).

Fig.8.Predicted subsidence in next 100 years for Wieliczka salt mine in meters(after [12]).

4.Concluding remarks

Geological and mining conditions of various types of deposits in Poland and China are complex,making the reliable prediction of subsidence and deformations very difficult.Original stochastic models (e.g.Knothe model)were developed for quite different mining conditions(shallow deposits and relatively slow extraction rate).The contemporary conditions are completely different due to the high concentration of extraction,depth over 1000 m and variable tilt of deposits.New methods of predicting surface deformations have been worked out,generalized for the exploitation of various types of deposits.These models also have been parameterized adequately to the local geological and mining conditions.

As shown in the paper,the advancement of works on the development and modification of the primary method and model of mining induced deformations in Poland and China allows modeling ground movements caused by the impact of coal,copper and salt underground mining.These methods stemmed from Knothe model.However,the preliminary assumption done by Knothe is no more reliable due to the changing mining condition.A new finding presented in the paper is the novel modification of the influence function method which brings much more accurate ground movement modeling for different types of minerals.Results of this research revealed that function methods are time-and cost-effective and reliable.Modified function methods could be applied for prediction of ground movements caused by every type of mineral reservoir exploitation.

Results of the presented research benefits for authorities and mining companies which are faced with the problem of induced ground movements.Especially,the presented solution may be of interest communities which are intensely build-up and where the ground movements are significant threat to buildings and infrastructures.The presented prediction model may serve as a novel tool supporting safety management on the areas hazarded by ground movements.

Acknowledgements

This paper is funded by the national key project ‘‘The Belt and Road”talent recruitment project named:Comparison of Mining Subsidence Research in China and Poland (No.G2017001).

Part of the research was financed from the Grant for Statutory Research AGH-University of Science and Technology in Krakow,Poland No.16.16.150.545.