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

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.

Rock column model in FLAC3D

An equivalent analytical solution is used for verification. The relationship between vertical stress increment Δσ_zz and vertical strain increment Δε_zz is:

where, α is the Biot coefficient; Δp 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, ρ_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.

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

2.1.1 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 HM-coupled scheme with hydraulic parameters as shown in Table 3.

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.

Displacement at the top of the model versus water level (full HM-coupling)

Figure 3 shows the results in terms of uplift, water pressure and stresses when water level has reached the final height of 100 m.

Vertical displacement, pore pressure and vertical stresses in rock column for final water level of 100 m (full coupling)

To verify these simulations against the corresponding analytical solution a Biot coefficient α = 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 parameters to actual mechanical conditions like stresses in case data/measurements are available.

Heave at the top of the rock column as function of water level (analytical vs. numerical solution)

(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.

Pore pressure and vertical displacement due to complete flooding of rock column (solid line = uplift, dashed line = pore pressure)

(3) Effective density scheme

The effective density (density under consideration of buoyancy) of rock after saturation is ρ_sat − ρ_w, which means the density of rock below the water level changes to ρ_s − (1 − n)ρ_w, because ρ_sat − ρ_w = (ρ_s + nρ_w) − ρ_w, where, ρ_sat is the saturated density. According to Archimedes' theory, the buoyancy force per unit rock volume (V_rock) is ρ_wgV_disp = ρ_wg(1 − n)V_rock, where, V_disp is the volume of displaced water. The overall force is (ρ_sg − ρ_wg(1 − n))V_rock, which also means the effective density of ρ_s − (1 − n)ρ_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) × 1000 = 1620 kg/m³.

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

Numerical calculated uplift using the effective density scheme

(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 σ_effective = σ_total − p. 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).

Uplift and stresses in rock column obtained by effective stress scheme

2.1.2 Water inflow from vertical boundaries

The following type of modelling considers water inflow (rebound) from sidewalls of the rock column according to Fig. 8.

Water level rebound after stop of pumping

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.

Uplift induced by water inflow from vertical boundaries

2.1.3 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.

Water inrush from bottom

Initial saturation of the model is zero. Inflow rate at bottom is fixed. Different inflow rates are investigated: 6 × 10⁻⁹ m³/s to 1.4 × 10⁻⁸ m³/s with an interval of 2 × 10⁻⁹ m³/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.

Uplift induced by water inflow from bottom boundary

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 top of the model and pore pressure at bottom. This relationship is shown in Fig. 11c.

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 × 10⁻⁸ m³/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.

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

Confined aquifer inside a rock column under hydrostatic pressure (water level = 100 m)

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.

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 × 100 m × 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.

2.5D profile model for parameter sensitivity study

3.1 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.

Effect of constitutive model of strata on uplift at model surface

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.

3.2 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.

Effect of geometrical size of mine goaf on surface uplift

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.

3.3 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.

Effect of mechanical and hydraulic properties on surface uplift (normalized to reference model and parameters)

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

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. (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. (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. (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. (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.