Numerical simulation of flooding induced uplift for abandoned coal mines: simulation schemes and parameter sensitivity

Numerical simulation approaches have been widely applied to study mining induced subsidence, and they are potential methods to study the flooding induced uplift for abandoned mines. This paper gives an overview about different numerical approaches to simulate uplift induced by flooding abandoned underground mines, including three different hydraulic conditions, considering both unconfined and confined water conditions. Four basic simulation schemes using 1-dimensional rock column models verified by analytical solutions demonstrate these procedures. The results reveal that flooding induced uplift is mainly related to the pore pressure in the mine goaf. The parameter study documents that height and stiffness of the mine goaf have the strongest influence on maximum surface uplift.


Introduction
Mining induced subsidence has already been studied for decades. Several approaches are proposed, such as empirical methods (including graphical method and profile function method) (Asadi et al. 2004;Hood et al. 1983), influence function methods (Cui et al. 2013;Hejmanowski 2015;Hu and Lian 2015;Luo 2015) and physical models (Ghabraie et al. 2015;Wang et al. 2021;Ye et al. 2018;Zhang et al. 2020). With the ongoing development of computational hard-and software, in recent years numerical simulation techniques are applied to predict and analyze subsidence (Cheng et al. 2018;Fathi Salmi et al. 2017;Li et al. 2020;Sepehri et al. 2017).
Especially in recent years, increasing attention is paid to the ground surface uplift induced by the flooding of abandoned mines. There are also some empirical methods to study the flooding induced uplift (Bekendam and Pottgens 1995;Fenk 2009), but they have site-specific character and their use is restricted to areas having identical geological, mining and hydrological conditions. Also, they can describe only simple geometrical conditions. At the same time, there is a lack of widely-used influence function methods and physical models to handle the uplift process. The satellite remote sensing technique is recently applied to monitor the uplift process (Bateson et al. 2015;Caro Cuenca et al. 2013;Devleeschouwer et al. 2008). But the accuracy of satellite monitoring data is influenced by numerous factors and restricted in general. For instance, reflected radar signals are easily affected by the rural and agricultural environment (Ashrafianfar 2014) or the atmosphere (Schäfer 2012). Besides, data processing from interferogram to DEM (digital elevation model) can also produce errors (Walter 2012). Therefore, the use of satellite radar data is only recommended to determine the long-term trend of surface movement. Hence, numerical simulation is a feasible method to study the uplift process and to provide predictions.
Numerical simulations to study flooding induced uplift for abandoned mines involve the selection of appropriate simulation tools and schemes. Recent publications show first applications of numerical tools to simulate the uplift & Jian Zhao jian.zhao@doktorand.tu-freiberg.de process (Dudek et al. 2020;Todd et al. 2019;Zhao and Konietzky 2020). This paper provides an overview about different numerical simulation approaches including full coupling, partial coupling, effective density and effective stress method to simulate flooding induced uplift for abandoned coal mines. Parameter studies document the general influence of several parameters on the uplift process. Also several hydraulic conditions like initializing rising water level inside the model, water inflow from outer vertical boundaries and water injection from the bottom boundary are considered.
2 Numerical simulation schemes to simulate flooding induced uplift The different simulation schemes are demonstrated by using a rock column with height of 100 m and quadratic cross section with edge length of 10 m (Fig. 1a). The basic model parameters are shown in Table 1. To verify the pure mechanical response, the settlement of the rock column driven by gravity is documented in Fig. 1b. Please note that the bottom boundary is fixed and the vertical boundaries are fixed only in the horizontal direction (roller boundaries). The code FLAC 3D (Itasca 2017) is used for the numerical simulations. An equivalent analytical solution is used for verification. The relationship between vertical stress increment Dr zz and vertical strain increment De zz is: where, a is the Biot coefficient; Dp is the pore pressure increment; K is the bulk modulus; G is the shear modulus. For a dry model, the pore pressure increment is zero, then the above given equation can be expressed as: where, q s is the solid dry density; H m is the model height; Z is the position for displacement determination. By integrating the equation given above, the vertical displacement of the model induced by self-weight is given by the following equation: When Z = H m , the vertical displacement at the model top u top is calculated as: Using the parameters given above, the analytical solution as well as the numerical result give u top = 4.9 mm.

Uplift induced by unconfined (phreatic) aquifer
The rock column shown in Fig. 1a is used considering a rising water level. The implementation can be performed in three different ways: (i) initializing rising water level in model zones, (ii) prescribing rising water level along the sidewalls (boundaries) of the model, (iii) injecting water from the model bottom (source).

Initializing rising water level in model zones
Water level induced uplift simulation can be performed by full or partial (two-step) coupling as well as effective density and effective stress schemes. Table 2 describes the key mechanisms for these different simulation schemes. Whenever it is necessary to study the water flow in a rockmass, full coupling or two-step coupling is the only choice. However, if simulation of the fluid flow process itself is not the main task, the other schemes are recommended due to their much lower computational effort.
(1) Full coupling The rock column model is used to study the fully HMcoupled scheme with hydraulic parameters as shown in Table 3.  Density q s (kg/m 3 ) 2500 Numerical simulation of flooding induced uplift for abandoned coal mines: simulation schemes… 1239 The simulations are performed in such a way, that the water level is rising step-by-step up to 100 m. The obtained vertical deformation (heave) at the top of the rock column is shown in Fig. 2. Figure 3 shows the results in terms of uplift, water pressure and stresses when water level has reached the final height of 100 m.
To verify these simulations against the corresponding analytical solution a Biot coefficient a = 1 is used. According to Eq. (1), the relationship between stress, pore pressure and strain is given as follows: When the position for determination of vertical displacement Z varies between zero and actual height H w , Eq. (5) could be expressed as: where, H w is the actual water level height. After integration, the following equation is obtained: If Z is above the actual water level, the increments for vertical stress and pore pressure are both zero, which means the dry part does not contribute to the heave or with other words, the vertical displacement in the dry part corresponds to the maximum displacement in the saturated part, given by the following equation: If the model is fully saturated (H w = H m ), the vertical displacement at the top of model is: Inserting of the parameters into the above equation gives u top = 1.7 mm, which is consistent with the numerical simulation result of 1.7 mm, as shown in Fig. 3.
Also for other water levels, analytical and numerical solutions give identical results as documented in Fig. 4. Considering Fig. 4 and Eq. (9), for a 1-dimensional rock column model with a phreatic (unconfined) water level or aquifer, the relationship between uplift and water level is a quadratic function. Note, that the presented HM-coupling is restricted in that sense, that permeability and porosity are kept constant during the coupling. However, it is simple to include functions which relate changing hydraulic Effective density pure mechanical simulation distinguish between density above and below the actual water level Effective stress pure mechanical calculation impose water pore pressure perturbation according to hydrostatic level and apply effective stress concept according to actual water level, meanwhile consider the saturated density below the actual water level parameters to actual mechanical conditions like stresses in case data/measurements are available.
(2) Two-step coupling To minimize the computational effort, the coupling can be restricted to the final position of the rising water level only. The rising water level as function of time is calculated, but the mechanical response in terms of effective stresses only for the final position of the water level (Fig. 5). This figure demonstrates that uplift and pore pressure for the final stage are identical to the analytical solution as well as for the full HM-coupling.
(3) Effective density scheme The effective density (density under consideration of buoyancy) of rock after saturation is q sat -q w , which means the density of rock below the water level changes to q s -(1n)q w , because q sat -q w = (q s ? nq w ) -q w , where, q sat is the saturated density. According to Archimedes' theory, the buoyancy force per unit rock volume (V rock ) is q w gV disp = q w g(1n)V rock , where, V disp is the volume of displaced water. The overall force is (q sg -q w g(1n))V rock , which also means the effective density of q s -(1n)q w for the rock below the water level. For the case study used in this research, the effective density below the water table is 2500 -(1 -0.12) 9 1000 = 1620 kg/m 3 .
In this scheme, fluid flow and pore pressure variations are not considered, but the calculated uplift is the same as obtained by using the coupling methods: 1.7 mm (Fig. 6).
(4) Effective stress scheme The effective stress method (or pore pressure perturbation) is based on Terzaghi's theory. When considering variable pore pressure as perturbation, the effective stress is r effective = r totalp. Uplift and stresses obtained by this method are shown in Fig. 7. Uplift obtained by applying the effective stress scheme is also 1.7 mm. Pore pressure and stresses are all consistent with the results obtained by using the full coupling methods (Figs. 3, 7).

Water inflow from vertical boundaries
The following type of modelling considers water inflow (rebound) from sidewalls of the rock column according to Fig. 8. Initially the rock column is dry. A rising water level is assumed at the outer model boundaries. Water level is rising stepwise, realized by changing the hydraulic boundary conditions, which leads to a water flow from the boundaries into the rock column. Figure 9a documents that this procedure delivers the same results as the other procedures described above.

Water inflow from bottom boundary via water injection
The third approach considers a rising water level originating from water inflow through the bottom of the model as shown in Fig. 10. Initial saturation of the model is zero. Inflow rate at bottom is fixed. Different inflow rates are investigated: 6 9 10 -9 m 3 /s to 1.4 9 10 -8 m 3 /s with an interval of 2 9 10 -9 m 3 /s. Figure 11a reveals that vertical deformation at the surface of the model increases with rising water level inside the model. At the same water level, the uplift increases with increasing inflow rate. Relationships between maximum uplift at the top of the model, pore pressure at bottom of the model and inflow rate are illustrated in Fig. 11b. Figure 11b demonstrates that vertical deformation and pore pressure are both positively correlated with inflow rate, so there must be a certain function between uplift at  This linear relationship can be explained via analytical formulas. If pore pressure applied at the bottom of the model is larger than hydrostatic pore pressure (here 1 MPa), maximum uplift will be larger than induced by hydrostatic pore pressure. According to Eq. (5), the vertical displacement at any position could be written as: where, P b is pore pressure at the bottom of the model; zero pore pressure indicates the water level position (phreatic surface). When the water level reaches the model top surface (H w = H m ), vertical displacement is a function related to pore pressure at model bottom and position height, as given by the following equation: When Z = H m , vertical displacement at the top surface is: When inflow rate is 1 9 10 -8 m 3 /s and water level has reached the model top, pore pressure at bottom is about 10.78 MPa, substituting data into Eq. (12), u top-= 20.71 mm. The uplift obtained from the numerical simulation is about 18.72 mm as shown in Fig. 11c. This deviation can be minimized by finer meshing. In fact, if initial pore pressure at bottom is fixed at 10.78 MPa, and 0 MPa at top, final uplift displacement of the surface is very consistent with the analytical result.
If the model is drained at top boundary, pore pressure at top is always zero, then top surface will not continue to move upwards, hence, final uplift at top surface is also consistent with the analytical solution. This model is suitable for a phreatic (unconfined) aquifer, in which top surface is a permeable boundary. However, if the model is undrained at the top boundary, a non-zero pore pressure at the surface will develop (refer to Eq. (10)), so that displacement at top surface will continue to increase. This model is suitable for a confined aquifer, in which top surface is an impermeable boundary and may be connected to an aquitard or aquifuge.

Uplift induced by confined aquifer
Verification of the numerical simulation approach for a confined aquifer under hydrostatic pressure conditions is documented in this part. When a confined aquifer is located at a depth between 0 and 50 m and a hydrostatic water level (100 m) up to the top of the rock column model is assumed, numerically simulated uplift and pore pressure distribution are shown in Fig. 12.
For a confined aquifer, the uplift under hydrostatic water level is given as follows: where, P a and P b are pore pressures at the top and bottom of the confined aquifer; H c is the height of the confined aquifer layer; H w is the water level height.
For H c = 50 m and H w = H m , the analytical solution for the uplift is about 1.40 mm, which is close to the simulation result (Fig. 12). What's more, Eq. (13) reveals that the uplift is linearly related to the rising water level for the 1-dimensional rock column model connected to a confined aquifer.

Parameter study
To investigate the effect of a few key parameters on the uplift, numerical 2.5D profile models (two elements in direction perpendicular to the x-z plane) are used. The model size is 8000 m 9 100 m 9 1000 m, as shown in Fig. 13. The equivalent mine goaf is 2000 m long and 50 m high, located 800 m below surface. Zone edge length is 50 m. Bottom is fixed and roller boundaries are applied at the sidewalls.

Effect of constitutive model for strata
The effect of the constitutive model for the surrounding rock mass (strata) is studied by considering the elastic and the Mohr-Coulomb model. The equivalent mine goaf is always considered as elastic material in a simplified manner to represent the reduced stiffness and density (Zhao and Konietzky 2020). The mechanical parameters are shown in Table 4.
The hydraulic properties for wet mine goaf and surrounding rock strata are shown in Table 5. Confined water conditions are assumed, that means mine water only exists in the mine goaf and corresponding water pressure is hydrostatic related to the model surface. The effective stress approach according to Table 2 is applied. The influence of the constitutive model of rock strata on uplift is shown in Fig. 14. Figure 14 reveals that the effect of the constitutive model on uplift is limited because this process is nearly pure poro-elastic. However, in general we can conclude, that plastic response of the rockmass leads to reduced uplift compared with pure elastic response.

Effect of size and inclination of mine goaf
To document the effect of size (horizontal extension and height) as well as inclination of the mine goaf, the parameters given in Table 6 and Fig. 15d are used. Figure 15a, e indicate that maximum uplift reaches a peak value, when goaf length reaches a critical value (full mining size). Uplift is positively correlated to goaf height according to Fig. 15b, e. The sensitivity of uplift to goaf height is far greater than it is to any other geometrical parameters of the mine goaf (Fig. 15e). Besides, Fig. 15 also reveals that uplift is negatively correlated to mine dip angle, while positively to goaf depth. Figure 15d reveals that in case of inclined goafs the position of maximum heave is shifted, connected with an asymmetric uplift pattern.

Effect of mechanical and hydraulic properties
To investigate the effect of mechanical and hydraulic properties, the parameters shown in Table 7 were  considered. During the parameter study only one parameter was changed while the other parameters were kept constant according to Tables 4 and 5. Permeability was not considered, because it only has influence on the evolution in time, and the simulation scheme of effective stress does not need to consider permeability. Figure 16 illustrates that the equivalent elastic modulus of the wet mine goaf is more important for uplift than other mechanical and hydraulic properties. Otherwise, as porosity of wet mine goaf increases, wet density under full saturation gradually increases and consequently uplift shows a slight drop. This can be explained also via Eq. (13). Hence, uplift has a negative correlation with elastic modulus of wet mine goaf and porosity.
Apart from these two properties of the wet mine goaf, the influence of mechanical properties of the surrounding strata on uplift are weak. In detail, according to Fig. 16, there is nearly no effect of friction angle of strata and tensile strength on uplift. However, for extreme small elastic modulus of the surrounding rockmass around the mine goaf, flooding induced uplift is small (Fig. 16). The reason is that the larger uplift occurring in the roof zone of the mine goaf cannot be transferred to the ground surface through the overlying rockmass due to the low stiffness of the strata. For the cohesion of strata: when cohesion  reaches higher values, uplift becomes close or identical to the elastic solution (Fig. 14).

Conclusions
This study explains different simulation schemes for flooding induced uplift for abandoned mines via 1-dimensional numerical rock column models. A parameter study using 2.5D numerical models is used to investigate the effect of several parameters on the uplift pattern and corresponding magnitude. The main conclusions are as follows: (1) There are four schemes for simulating flooding induced uplift: full and partial coupling as well as effective density and effective stress method. When only focusing on ground surface movement and water level rising is known, it is not necessary to consider full or partial HM-coupling. Considering effective density or effective stress is sufficient and saves computational time and reduces simulation complexity.
(2) The investigated three different approaches to simulate the flooding process (rising water level is induced by initializing rising water level inside the goaf zones, inflow of water from vertical outer boundaries or water injecting from bottom boundary), provide identical results. (3) Under confined aquifer conditions, the flooding induced uplift is linearly related to the rising water level. While for unconfined aquifers, the uplift is represented by a quadratic polynomial related to the water level. (4) Compared with other impact factors, height of mine goaf and equivalent elastic modulus of the goaf have the strongest influence on maximum surface uplift. To obtain accurate simulation results, these two parameters should be calibrated accurately when carrying out numerical simulations of flooding induced uplift using large-scale complex geological models.