Underground extraction of coal mining could lead to serious environment problems because of surface movement and subsidence. The damage to environment could continue for a long time. A number of computational methods have been developed for simulating ground movement and deformation caused by underground mining (Braeuner 1973; Peng and Chyan 1981; Alejano et al. 1999; Zhao et al. 2004). The most widely used method is subsidence engineering handbook (SEH) for predicting subsidence based on the graphical solutions (Anon. 1975). The SEH method offers fast and easy calculations where empirical predictions are based on actual measured data. However, this method is only applicable for simple coal panel’s geometry and for solving two-dimensional (2D) problem only. The profile functions method have been successfully applied to predict subsidence profiles induced by working in horizontal or moderately dipping coal seams (Torano et al. 2003), but this method cannot be accounted for coal seams with irregular panel. The complex cases can be modeled by analytical methods such as finite element analysis which offers some possibilities for subsidence prediction by Najjar and Zaman (1993), Yang et al. (1993), but these are invariably fraught with some difficulties. Moreover, influence functions method is considered as a powerful application for three-dimensional (3D) subsidence prediction for all shapes of extraction panels (Sheorey et al. 2000). However, the use of influence functions is very time consuming and calibration is difficult. A sensitivity analysis method is proposed for the GIS-based mapping of the ground subsidence hazard near abandoned underground coal mines (Oh et al. 2011).

With the development of computer technologies, the computational method is widely accepted. Most research has focused on final subsidence at the center line above mining operations to simulate and assess the surface structural damage without considering the dynamic extraction process. However, the occurrence of ground movement and deformation caused by underground extraction has a relationship with its location and mining sequence obviously. It is defined as a four dimensional (4D) problem and each of sub-components has a different effect on structures subject to subsidence. It is suggested that the location of mined panels is the most important parameter influencing coal mining subsidence (Suh et al. 2016).

The influence factors of subsidence work alternately. For example, horizontal strain is a major component in cracking, tensile failure of concrete dam structure, and leakage of water from a reservoir when it is under mined. Since it is common for a structure to be subjected to compression and extension strain in different directions, its response is not a simple problem of predicting to a particular numerical value of subsidence. It should be subjected to alternating stages of strain owing to the sequence of extraction, and a spatial–temporal (4 dimension) effect must therefore be considered.

In a 4D numerical simulation of a mining area with multiple extraction panels in a complicated mining sequence, it is difficult to obtain the component distribution of ground subsidence with the current existing prediction methods. Therefore, it is not possible to assess structural damage accurately. Numerical methods that could simulate the movement considering spatially and temporal processes are desirable for the reliable design of the mining layout in order to minimize structural damage. In considering the development of the computer-based prediction methodology, it is important to predict 4D ground subsidence components at any point with any shape of excavation so as to cover a wide range of mining geometry as well as to provide automation, intellection, and visualization in the simulation of subsidence process. The recent development of geographical information systems (GIS) comprises a technology designed to support integrative modeling, to conduct interactive spatial analysis and for understanding various processes. In case of a progressive ground subsidence simulation during undermining, GIS would be effective and efficient in computing such 4D ground movement and deformation if a prediction system could be properly developed. It also should be addressed that the GIS approach could result in an uncertain issue once it relies on just a small amount of information (Longoni et al. 2016).

In our present research, a GIS based computational methodology is developed for calculating the distribution of ground movement and deformation at arbitrary surface points (Djamaluddin et al. 2012). The GIS based numerical simulation is validated by two study cases. First is the prediction of 21-year of ground subsidence due to complex underground mining geometry in Japan. Second is the assessment of dynamic ground movement due to underground mining sequence in China. The tight coupling model, the component object model (COM) using GIS, is used to develop a GIS based tool for ground movement prediction system. Therefore in this paper, first, the fundamental calculation of ground movement over time using the Knothe time model and stochastic prediction model is briefly introduced. Second, a GIS-based prediction methodology for four-dimensional analysis is described. Third, two case studies are given for the validations the GIS-based computational model. Using GIS COM protocol, a computational implementation of the developed tool is discussed for prediction of 4D ground subsidence from underground mining sequences. Finally, the computational system software of the developed tool within GIS is presented.

2.1 Dynamic calculation method for ground movement

In an underground excavation progress, the excavated zone greatly influences both the rock mass and the ground surface. Hence, the final ground movement that takes place at a given point depends on time and the creep characteristic of the rock mass undergoes. The final value of the subsidence mainly depends on the size and position of the mining extraction. An analytical method was suggested to describe the creep characteristic (Knothe 1953). Furthermore, Knothe assumed that the rate of subsidence of a point is proportional to the difference between the possible final subsidence Send of this point, owing to extracting a portion of the instantaneous movement of this point at the time that S_dyn is considered (Fig. 1), which is

where c is the factor of proportionality, which is related to the physical and mechanical properties of the overburden rock and soil; \(s_{t}\) is the time-delay factor; \(t\) is the time.

Illustration of dynamic subsidence analysis of a surface point P by extraction sequence in relation to the time factor

The coefficient c is a parameter to describe the influence of mining and geological conditions on the movement process in time. By using the Knothe (1957) model, the ground subsidence at any point could be calculated when the mining ongoing and it is also possible to obtain the other time dependent variation of movement.

The Knothe time model is adopted in this research for the dynamic ground movement prediction, and the key process is to determine the time factor coefficient c. The most significant factors influencing the ground movement were summarized as followings (Whittaker and Reddish 1989):

• Distance from the working face;

• The thickness of overburden;

• The rate of advance of the working face;

• The method of the working panel; and

• The geological structure.

This is not an easy problem for determining the above influences and calculating the factor of proportionality. However, with a series of measured time-subsidence curves in the ground surface, the coefficient c, describing the ground movement rate, could be determined (Berry 1977; Burns 1981). The dynamic ground movement at a time of a given point could be formulated in fractions if its final subsidence according to Kratzsch (1983).

To obtain the coefficient c, the time-movement curves of the measured data from in situ investigations of subsidence should be draw at first. Second, by knowing the final subsidence rate, the progressive subsidence data, and their interval time, the coefficient c can be calculated by back-analysis of Knothe formulae. Finally, the typical coefficient c in the corresponding mining regions could be derived by analyzing several numbers of the measured stations.

2.2 Stochastic medium theory and spatial subsidence prediction

In all prediction methods, the assumption of rock mass behavior generally represents continuous or discontinuous materials. Two concepts have been devised with special regard to the mining subsidence process and seem to be adaptable especially to discontinuous materials. The first concept is the rock mass as a stochastic medium and second concept is the gap-diffusion as reported by (Braeuner 1973). Because the overlying strata behave in a complex manner, and the movement of the rock mass is governed by a number of known and unknown factors, a stochastic medium theory is a widely accepted model for the prediction of 3D ground movement (Litwiniszyn 1957). Based on stochastic medium theory, a series of solutions for subsidence calculations in different geological and extraction conditions have been obtained in China and Japan. The stochastic solutions have been adopted in mining practice and constructions of underground space to solve the excavation problem under railways, rivers and buildings.

A stochastic medium model is a tool for estimating probability distributions of potential outcomes by allowing for random variation in one or more inputs over time. To calculate movement of a surface point P using the stochastic model, an excavation panel can be divided into infinitesimal areas. According to the principle of integrated subsidence effect of the excavation panels, the consequence would be equal to the sum of the effects caused by those infinitesimal areas. Based on the stochastic medium theory, the occurrence of a rock-mass movement over the extraction element may be a random event that takes place with a certain probability. A unit with infinitesimal width, length and thickness (\(\partial w\partial l\partial m\)) in an extraction panel is called the extraction element. The vertical displacement at any point in the movement slice is defined as the basic influence function (\(S_{e}\)). An event in which surface movements take place in an infinitesimal area, \(dA = dxdy\), at horizon z, with point \(P(x,y,z)\) at its center, is equivalent to the simultaneous occurrence of two events composed of a movement in the horizontal strip \(dx\) and the horizontal strip \(dy\) through point P (Fig. 2). Fundamentally, the probability can be written separately for these two events by \(C(x^{2} )\) \(dx\) and \(C(y^{2} )\) \(dy\), respectively, where C is the subsidence slice function. The probability for a simultaneous occurrence of these two events is

Illustration of the probability prediction of ground movement at a point as a result of the extraction element within a given extraction panel

For calculating the other components of ground movement, such as vertical displacement, horizontal displacement, curvature, slope and strain, calculation procedures were also used. The vertical displacement or subsidence at point P is given:

The slope at point P along direction \(\phi\) is given:

The curvature at point P along direction \(\phi\) is given:

The horizontal displacement at point P along direction \(\phi\) is given:

The horizontal strain at point P along direction \(\phi\) is given:

where \(S_{\hbox{max} }\) is the maximum possibility subsidence; \(e_{n}\) is subsidence influence, depending on the extraction area (\(C_{x} \times C_{y}\));\(m\) is the coal-seam thickness; \(a\) is the subsidence factor; \(\alpha\) is the angle of dip; \(l\) is the panel length along strike; \(w\) is the panel width along dip;\(r = H/\tan \gamma\), \(r_{1} = H_{1} /\tan \gamma_{r}\), \(r_{2} = H_{2} /\tan \gamma_{d}\), radius of the circle of influence; \(\gamma_{r}\) is the angle of draw to the rise,\(\gamma_{d}\) is the angle of draw to the dip; and \(H\) is the depth along strike, \(H_{1}\) is the depth along the boundary of the rise side, \(H_{2}\) is the depth along the boundary of the dip side.

3.1 A strategy to integrate GIS and models

Based on the classic calculation method by Knothe function and Stochastic medium theory, a new spatial model for subsidence of underground mining could be developed based on GIS. GIS provides a wonderful platform for dealing with spatial data and graphical output. However, the general purpose of GIS, which provides only a basic tool, cannot be employed to model specific problems. How to integrate the subsidence-prediction modeling to GIS is a question to be solved. It is an information-integration problem, a little like combining one GIS to another for data-transfer purposes. For example, analytical subsidence-calculation points that represent ground movements for application to a simulation model can be designed to link to GIS automated directly. At the same time, a modeling study undertaking in GIS supplies a basis for simplification of the interaction between the different users involved through the establishment of a common data structure that can be visualized using the same GIS-based system. Joining the variety of data, models, and tools into a robust system of GIS is a research topic that is approached ranging from so called “loose” integration to “tight” integration.

A strategy based on COM technology coupling with GIS was applied to overcome the problems of complex geometric modeling of mining panel and data integration in present study. The coupling strategy includes integrated data management services of GIS, and automated exchange of data becomes possible through a standardized interface using COM method. COM is a standard, which enhances software interoperability by allowing different GIS components, possibly written in different programming languages, to communicate directly (Matthew and Michael 2002). The integrated model is accomplished within the GIS system to achieve an effective and efficient calculation method.

3.2 Integrated prediction model within GIS

By integrating the effect of all extraction elements in an excavation panel as shown in Fig. 3, all the subsidence components (dS _p ) related to a corresponding surface point (S _p ), as illustrated in the figure by a 3D view of grid points (calculation points) and a 3D polygon (extraction panel), may be calculated. For the extraction panel, with reference to the vector-based polygon (Fig. 4), the spatial data of panel geometry, mining sequence, ground movement parameters, and excavation depth can be stored in the 3D polygon. A feature table is used to relate the subsidence parameters in the panel data set of polygon. In the 3D polygon attribute table, ‘PolygonZM’ is the shape of 3D polygon attributes, and ‘ID’ is the extraction sequence. ‘AngDip’ is dip inclination of seam panel, ‘UpwardAng’ is the upward angle of the panel from the east, and ‘Thick’ is the extraction thickness. All of these are related to the geometrical coordinates. The corresponding parameters of subsidence are represented by ‘SubFac’ (subsidence factor), ‘HoMoFac’ (horizontal movement factor), ‘TiFac’ (time-delayed subsidence factor), ‘UpTan’ (tangent of draw angle in rise side), and ‘DownTan’ (tangent of draw angle in dip side). The values of depth vertices are stored in each 3D polygon that give spatial geometry in x, y, and z. In other words, a 3D polygon has spatial geometry that can be used to identify the main strike direction and dip inclination of each panel. Together with a spatial model in GIS, a triangulated irregular network (TIN) model is adopted to identify the strike direction and dip inclination of an excavation panel.

Probability distribution of dynamic movement at a grid point as a result of given 3D polygon extraction panels

Example of 3D polygon panels with related spatial data geometry and features table of data set attributes

Subsidence at the surface points is calculated from the 3D polygon panels as geometrical excavation areas, and each panel of the polygon is referenced to the global coordinate system in GIS. It is assumed that the inclination of a 3D polygon has an upward angle direction (\(\varphi\)) with a reference from the east. To obtain distribution of surface subsidence from inclined panels, the global coordinate panel (X, Y) is transformed to local coordinate (X″, Y″) in which the upward angle direction of each panel is set to be the same as the east direction. The radius of subsidence influence circle of each panel is assumed to be the downward and upward part of the panel, and the main direction of the radius influence circle is set to be the same as the upward angle of the panel. An example of a polygon panel before and after coordinate transformation is shown in Fig. 5a. The coordinate transform polygon panel vertex (x″, y″) that gives the depth can be performed by simple equations. Let global coordinate X, Y and its transform coordinate X″, Y″ be derived as follows:

where (X, Y) are for global coordinates, and (X″, Y″) are for local coordinates.

a Coordinate transformation of a 3D polygon panel. b Obtaining the radius of subsidence influence circles on the transformed panel

The main direction angle of the upward panel is set as the subsidence-influence circle direction. The main direction of the upward angle (\(\varphi\)) and the main dip inclination of the panel (\(\alpha\)) could be obtained from the depth value of the vertices on a panel in the polygon area. Referring to the coordinate transform, it is able to get the minimum and maximum values of the panel vertices in the x-coordinate. The dip inclination of the panel is important for the calculation of the subsidence-influence circle area of each panel. It should be determined at the first step. Referencing known x-coordinate values after transformation, the radii of the subsidence-influence panel circles could be calculated by the given depths. Figure 5(b) is an example of a panel after coordinate transformation, and a section plan (X″, Z) that shows the vertices values of the inclined panel’s depth. The radius distances at the upward (r ₁ ) and downward (r ₂ ) parts of the panel are effected by the inclined dip angle α, and the solution equations could be derived as follows:

where \(H_{i}\) is the depth of the panel vertex; \(X_{\hbox{min} }\), \(X_{\hbox{max} }\) are the minimum and maximum values of the x-coordinate of panel vertices for local coordinates; \(X_{i}\) is the x-coordinate of the panel vertex; \(\alpha\) dip inclination of the panel; \(\gamma_{r}\), \(\gamma_{d}\) are the angles of draw for the upward and downward parts of the panel, respectively.

3.3 Computational method algorithm

An flow chart is shown in Fig. 6 to demonstrate the algorism of prediction analysis using GIS functions to calculate a surface subsidence. In the computational process of ground movement, the panel on the global coordinate system is transformed into the local coordinate system at first. The transformation coordinate panel is performed to get the radius and angle direction of the main subsidence zone. Then, the prediction of subsidence distribution is carried out in the stochastic-prediction procedures. Finally, the calculated ground movement results are summarized according to the working panel number.

Computational processes for 3D dynamic movement prediction using GIS

In order to predict the ground movement using the GIS functions, the following data, such as calculation points, extraction panels and subsidence parameters should be prepared at first. Meanwhile, the input parameters for ground movement calculations could be separated into two categories. The first one is the generation of surface grid points for providing calculation points, including the distance of each grid point in the x and y directions, in which the number of points along the grid lines, the number of grid lines, point intervals along the grid lines, and the grid-line directions need to be established. The second one is entering the calculation data of ground movement, including the extraction panel sequence ID, seam upward angle, inclination of dip, subsidence factor, excavation thickness, horizontal-movement factor, angle of draw and time factor. Since the ground movement is calculated only with the above entering data, the calculation results of each panel would be uncertain once the input data are uncertain.

3.4 Computational processes using spatial model

The ground movement computation is not limited within the GIS. It also could be performed outside the GIS. In this case, the GIS system would be used only as a spatial-related database of ground movement for storing, displaying, and updating the inputted data. The main advantage of this approach, using an existing external subsidence-prediction model, is to save time in programming the model algorithms into the GIS. A disadvantage of this method is the complication caused by the conversion of complex geometrical data to and from external models. Prediction models calculate the movements of a surface point for an extraction panel in three dimensions. Because an individual mine may consists of multiple panels and complex geometry, the use of the 3D prediction model for obtaining the spatial distribution of subsidence is very time-consuming without the GIS, as each mining panel has to be calculated separately. To overcome the data conversion problem of complex geometrical, the calculation model of ground movement can be established within the GIS. All modules related to the GIS spatial-analysis function are shown in Fig. 7. Every GIS function is implemented by a GIS component. All of the ground movement-related parameters such as spatial geometry, surface-point calculation data, and the subsidence parameters could be obtained from the functions of the GIS data module. The surface subsidence as well as vertical displacement, horizontal displacement, slope, curvature, and horizontal strain are calculated with the stochastic prediction module. Finally, the surface subsidence grid-point calculation and its subsidence-prediction kriging interpolation can be obtained by a function of the GIS spatial analysis. Because a GIS component is implemented in the calculation model, the three-dimension problem can be computed effectively.

Processes for dynamic subsidence prediction

3.5 Computational workflow

A three-step process is applicable in order to implement a subsidence-prediction model in the GIS environment. In step 1, it includes deciding on an appropriate working projection, establishing spatial extents of the study area, and assembling variously used spatial data from the study area so that the spatial component can be overlapped correctly. This step combines the appropriate subsidence data together into a GIS spatial database such as a mining-panel map, a seam layer, a surface-elevation layer, and spatial information from overburden strata. Step 2 is to make the spatial geometry of the mining panel, in which the ground movement calculation parameters are stored. Usually the mining-panel map is digitized manually, transformed from an image into GIS vector data. In order to construct a geometry representing the mining panel to store the calculation parameters for the ground movement computational model, a 3D polygon is adopted and each geometry panel is required for the depth of panel vertices (co-ordinate in the z direction). The prediction value of subsidence is stored into 3D polygon once it is arrived at correctly, and it is combined with the integration of data in step 3.

4.1 Subsidence resulting from 21 years of mining beneath a reservoir

According to the results of previous studies, including a sequential calculations of ground movement and a series of excavation process simulations between 1944 and 1967, the following parameters are obtained:

• subsidence factor (a) = 0.85,

• horizontal movement factor (b) = 0.21,

• tangent of draw angle (γ) = 1.428,

The time-delayed subsidence is represented by the time factor. It shows that the final value completed in 3 years (first year = 0.83, second year = 0.90, third year = 1). The subsidence continued even after coal mining was completed, and the maximum value was 3.27 m up to 1967. The progressive subsidence simulation results of past mining stages including 1950, 1961 and 1967, are shown in Fig. 8a–c. For these spatial simulations, all calculation points were interpolated by using the kriging interpolation method. The simulated maximum vertical displacement is 1364 mm in 1944, and becomes to 1994 mm after 6 years of mining (in 1950). The magnitude is 2505, 2576 and 3012 at 1955, 1958 and 1961, respectively. At the same time, the zero subsidence contours are indicated by the dotted line. This line encloses the subsidence area and identifies the beginning of subsidence and the limit of ground movements that may cause damage to surface structures. It is found that a surface subsidence increased progressively about from 1956 to 1967 at the reservoir area. Eight surveys were carried out during the twenty one year mining period to monitor subsidence data around the reservoir every year.

a Simulation of subsidence distribution for mining up to 1950. b Simulation of subsidence distribution for mining up to 1961. c Simulation of subsidence distribution for mining up to 1967

4.2 Subsidence from a coal-mining sequence beneath railway and road lines

A coal mine in China has been studied. From September 20, 1999, to July 20, 2000, the longwall face of no. 1315 is excavated. The average mining depth is 269 m below the ground surface. The overburden consists mainly of clay, sand, and sandstone of varying thicknesses. The surface topography of the coal field area is almost flat. The coal seam thickness is 5.6 m. The average angle of dip of the coal seam is 4°, and the average mining rate of the longwall face is 120 m/month. The strike and dip length of the working face are 1160 m and 150 m, respectively. The monitoring points were arranged along a railway line and a small road line with a total of 60 measurement positions. The layout of coal panel and the measurement points is illustrated in Fig. 9. The time factor of the subsidence can be calculated by the Knothe time model. The calculation of the time coefficient (c) is 13 year⁻¹, which is determined from the measured points. The time factors for the relevant duration of influence of the coal panels are:

Plan layouts of measuring points and mining extraction progress

• First month: \(\left( {1 - e^{{ - 13 \cdot \frac{1}{12}}} } \right)\) = 0.662

• Second month: \(\left( {1 - e^{{ - 13 \cdot \frac{2}{12}}} } \right)\) = 0.885

• Third month: \(\left( {1 - e^{{ - 13 \cdot \frac{3}{12}}} } \right)\) = 1.000

The general subsidence parameters could be determined by the observations and borehole data. As a result, subsidence factor is 0.8, the horizontal displacement factor 0.2, and the tangents of both the draw angle at dip side are 2.6 and of the angle at rise side is 2.2. The progressive surface subsidence owing to the mining sequence can be predicted in terms of monthly sequences. Figure 10a–c illustrate the comparison of measured and calculated subsidence point values at the road line in the months of January, February, and March, respectively.

a Comparison of measured and predicted subsidence curves with mining progress in January. b Comparison of measured and predicted subsidence curves with mining progress in February. c Comparison of measured and predicted subsidence curves with mining progress in March

5.1 Structure of the integrated system

An integrated system based on GIS platform has been developed to predict subsidence due to underground mining in which all calculations and data processing have been implemented in a computer program as an extension tool in ArcGIS, called as mining subsidence damage assessment system (MSDAS-GIS). Figure 11 illustrates the whole system structure in which a GIS component is used to fulfill the analysis of spatial–temporal. All the subsidence and environmental impact analysis related GIS data can be managed and analyzed effectively as the same as the ordinary GIS software. Spatial analyst and 3D analyst are used to analyze 3D-polygon panel geometry and to provide more accurate possibility of input parameters. All data for the subsidence calculations are in GIS vector data, and the final calculation results can be transformed into GIS raster data. At the same time, the composite algorithms and iteration procedures of the dynamic subsidence analysis problem can also be implemented perfectly.

The main system of MSDAS-GIS toolbar in GIS platform based on ArcGIS of ArcMap technology

In this research a tight coupling method based on COM technology has been used to overcome the problems of GIS-based model integration. Figure 12 illustrates the coupling models of GIS in which a COM method is used to communicate between the models and the GIS components. The subsidence—analysis—related GIS data can be managed and analyzed effectively in the same manner as ordinary GIS software. In order to achieve a more accurate possibility of input parameters, the tools of 3D analyst and spatial analyst are used to analyze the spatial geometry of mining panel. With the proposed method, the composite algorithms and iteration procedures of the dynamic subsidence analysis problem can also be implemented successfully. The concept of a GIS-based component model is that developed system can be assembled with reusable model components. Dynamic models of impact assessment and prediction systems often have common components that usually represent a specific aspect or functionality of the system. GIS software is built based on common modular of shared components. These components are interchangeable which can be replaced by similar components depending on how the user wants to customize the model. After some specified recoding or modification, the related functional parts can be built with standard interfaces and become reusable. These reusable components form a resource that is available for creating customized simulation models. Different components may come from a same parent model and different models may provide similar components, which have different internal structures and calculations but perform the same functionality. In the system of MSDAS-GIS, the component model is employed for fulfilling all the GIS functions.

Illustration of object inside a GIS ArcMap is accessed by COM method

5.2 Computational process and system interface

In this part, the whole computational processes of mining subsidence analysis are expressed in more detail. The subsidence data modeling is established as an integral part of the GIS, and the developed system provides a way for subsidence prediction as well as environmental impact assessment analysis. Figure 12 shows the main interface of MSDAS-GIS toolbar within GIS. It could be divided into three categories that have relationship with the implementation calculation process of mining subsidence. The system provides a set of supporting tools that is using the complex data preparation while enabling the utilization of the spatial and 3D analyst supported by GIS.

In order to perform a subsidence prediction analysis in GIS, there are four steps named as general setting, pre-processing, subsidence calculation and post-processing. The interface of the pre-processing of subsidence calculation is given in Fig. 13a. This form is used to prepare a 3D GIS-polygon dataset in which all the subsidence parameters and geometrical parameters and extraction sequence information are included. MSDAS-GIS provides tools for creating a new dataset of 3D GIS-polygon, identifying its polygon depth as well as panel vertices number. There are five steps in subsidence parameters preparation that is necessary to perform the parameter-related 3D GIS-polygon. In order to calculate subsidence of a grid point, a program was developed to generate surface grid points for providing calculation points, including the distance of each grid point in x and y directions, where the number of point along grid line, number grid line, point interval along grid line, and grid line direction are needed to be established. The form interface for generating a grid point is shown in Fig. 13b. Figure 13c shows the form interface for calculating surface subsidence due mining sequence. Due to panel sequence related subsidence data and multiple calculation panel, it is difficult to manage all of these data without 3D-polygon-based datasets, thus in the calculation system, a polygon dataset is used to store all panel datasets as a spatial geometry.

a The form interface for “Pre-processing” subsidence prediction. b The form interface for “Generating grid” calculation point. c Form interface for “Calculation core” subsidence prediction based on stochastic method. d The form interface for “Post-processing” subsidence prediction

In this 3D-polygon panel dataset, a feature table is used to relate the subsidence parameters and spatial geometry. Developed system provides form interface in the calculation to make an easy calculation phase process and calculation of 3D subsidence components in arbitrary direction.

The form interface for post-processing subsidence prediction is given in Fig. 13d. All the calculation results can be stored into GIS point-grid. Each calculation point such as vertical displacement, slope, curvature, horizontal displacement and horizontal strain can be transformed into a GIS raster by a surface interpolation as well. Simulation ground movements for each calculation number are possibly shown in this interface.

In subsidence calculation profiles, the calculated components and phase numbers of working panels will be shown on “MapView” on the main interface, as shown in Fig. 14, and the ground movement’s components changing will be shown by two dimensional (2D)-graph.

An example of the calculating process for dynamic movements showing the simulation of five components subsidence movements along the centerline of longwall mining

This paper has demonstrated a GIS based developed tool to predict 4D subsidence from underground mining sequence. The strategy of close coupling between a prediction model and GIS using COM technology is used. Using a GIS-based prediction methodology, it is possible to calculate spatial and temporal ground-surface movements induced by multi-seam longwall mining. Differential-subsidence characteristics, such as progressive vertical displacement, slope, curvature, horizontal displacement, and horizontal strain, can be computed using the GIS developed tool.

Subsidence in the study area of Japan is simulated. The calculated result has been used to evaluate subsidence-induced damage resulting from mining beneath a reservoir. The progressive ground subsidence resulting from underground mining sequence in the study area of China is also simulated. The calculated results have been used to compare observed subsidence values at road and railway lines. It should be addressed that the GIS related analysis result has an uncertain issue because it relies on a small amount of information. The accuracy of the simulation results depends on the original GIS data.

As the future study, the GIS-based developed tool should be improved to consider the effect of steeply mining, heterogeneity of rock mass, and influence of fault and surface topography conditions. Moreover, a multi-layer approach should be performed to determine diverse factors such as geological structure and ground water level, which can improve the GIS subsidence prediction model.