A novel ground surface subsidence prediction model for sub-critical mining in the geological condition of a thick alluvium layer

Abstract

A substantial number of the coal mines in China are in the geological condition of thick alluvium layer. Under these circumstances, it does not make sense to predict ground surface subsidence and other deformations by using conventional prediction models. This paper presents a novel ground surface subsidence prediction model for sub-critical mining in the geological condition of thick alluvium layer. The geological composition and mechanical properties of thick alluvium is regarded as a random medium, as are the uniformly distributed loads on rock mass; however, the overburden of the rock mass in the bending zone is looked upon as a hard stratum controlling the ground surface subsidence. The different subsidence and displacement mechanisms for the rock mass and the thick alluvium layer are respectively considered and described in this model, which indicates satisfactory performances in a practical prediction case.

Keywords

ground surface subsidence / thick alluvium layer / sub-critical mining / prediction model

Introduction

The vast majority of Chinese coal reserves are distributed in North and Northwest China, where alluvium layers commonly cover the bedrock strata. Of particular note, the thickness of the alluvium layer is up to 100–300 m in the famous Loess Plateau and North Plain,. Thus it is critical that researchers and engineers predict the ground surface subsidence in the geological condition of thick alluvium layer. It is difficult to obtain satisfactory prediction results even though a number of subsidence prediction models have been developed over the past decades ( Knothe, 1957; Litwiniszyn, 1957; Awershen, 1959; Liu and Liao, 1965; National Coal Board Mining Department, 1975; Kratzsch, 1983; Peng, 1984; Ren et al., 1987; Lin et al., 1992; Liao, 1993; Cui et al., 2000; Sheorey et al., 2000; Álvarez-Fernández et al., 2005), . This is due to the fact that the overburden strata have been generally assumed to be a simple medium or a “comprehensive medium” in most conventional prediction models. This assumption is relatively reasonable when the overburden strata exhibit the same or similar mechanical behaviors. However, the alluvium layer, i.e., the third or quaternary stratum, is mainly composed of loess, alluvium, sand, gravel, etc. There are, in fact, significant differences between the mechanical behaviors of rock strata and alluvium layer so that the above assumption is too simple to reflect the displacement mechanism inside the overburden rock strata and the alluvium layer. It was not until 1999 that a ground surface subsidence prediction model was proposed ( Zhang et al., 1999) for the critical and super-critical mining in the geological condition of thick alluvium layer. However, this prediction model is not appropriate for sub-critical mining since the ground surface subsidence in different mining degrees obeys different laws. Therefore, an accurate prediction model for sub-critical mining in the condition of thick alluvium layer is urgently needed in China.

Because of the existence of a thick alluvium layer covering on the overburden rock, the law of the overburden rock displacement and ground surface subsidence depends not only on the properties of the overburden rock, but also on the properties of the thick alluvium layer. We will try to establish in this paper a novel ground surface subsidence prediction model suitable for the sub-critical mining in the geological condition of thick alluvium layer by investigating the displacement mechanisms in the overburden rock and the thick alluvium layer.

Basic assumptions for coal seam, overburden rock and alluvium layer

To establish a specific ground surface subsidence prediction model suitable for underground mining in the geological condition of the thick alluvium layer, we make the following assumptions for the coal seam, the overburden rock strata. and the thick alluvium layer:

The characteristics of the coal seam and the overburden rock strata are essentially horizontal, with the overburden rock strata exhibiting the similar mechanical behaviors.

In view of the mechanical behavior of weak compactness and looseness, the thick alluvium layer is regarded as random medium as are the uniformly distributed loads acting on the overburden rock mass.

Formation of a fissure arch and sub-critical mining

As is well known, the immediate roof and the main roof will successively collapse, and the fragmentized rock bulks will fall into the goaf when the underground coal seam is excavated. The displacement in the overburden rock manifests itself in the form of collapse, fracturing, bed separation, bending, and the consequential ground surface subsidence ( Kratzsch, 1983; Peng, 1984, Dai et al., 2002). Numerous fundamental investigations and simulation tests of similar materials demonstrate that a fissure arch will form above the excavation zone in the overburdened rock, and the bending rock strata will have a tendency to bridge over the goaf by arching between the solid abutments on each side with the advancing of the working panel ( Wang et al., 1997; Sansone and Ayres Da Silva, 1998; Wang et al., 2004). Figure 1 shows a schematic diagram regarding the formation process of a fissure arch induced by the underground mining ( Wang et al., 2004). When the excavation dimensions reach a certain size, the height of the fissure arch will no longer increase, but will expand horizontally. If the overburden rock strata are horizontal, the top trace of the fissure arch will show a straight line, separating the bending and the fractured zones ( Wang et al., 1997; Wang et al., 2004; Lai et al., 2006); thus, the trace of the fissure arch serves as both the separation and the boundary between the two zones. At this time, the rock mass in the bending zone (above the fractured zone) is primarily in the elastic state ( National Coal Board Mining Department, 1975; Kratzsch, 1983; Peng, 1984; He and Yang, 1991; Bai et al., 1995). However, with the goaf dimension continuously increasing, the bending rock strata are no longer able to arch over the goaf so that the maximum ground surface subsidence occurs. The sub-critical mining is typically referred to as the mining activity before the ground surface subsidence reaches maximum value. Hereafter, we will focus on the ground surface subsidence and other deformations in the sub-critical mining stage.

Subsidence of overburden rock in bending zone

A zone exists in the overburdened rock above the caved long wall face where large displacements have taken place. A major portion of ground surface subsidence occurs as a result of this zone ( Oravecz, 1986) which in actuality is the bending zone above the excavation zone. In the geological condition of the thick alluvium layer, the bending zone is referred to as the overburden rock beneath the thick alluvium layer. We can now investigate the subsidence of the overburden rock in the bending zone (see the zone above the fissure arch in Fig. 1), and consider it as a type of “excavation.” It is the “excavation” that leads to the ground surface subsidence through the thick alluvium layer.

The subsidence of the overburden rock in the bending zone w_b is actually the sum of the rock mass compression deformation on the excavation boundary and the deflection of rock stratum in the bending zone ( Su et al., 2003), namely,

(1) w_b=w_{\mathrm{1}}+w_{\mathrm{2}},

where w₁ is the rock mass compression deformation on the excavation boundary; and w₂ is the deflection of the rock stratum in the bending zone.

Rock mass compression deformation on the excavation boundary

In this section, we will examine a unit inside the rock mass as shown on the left side of Fig. 2. In the original stress state, the vertical stress {\sigma }_z^{\mathrm{0}} , and the lateral stresses {\sigma }_x^{\mathrm{0}} and {\sigma }_y^{\mathrm{0}} can be respectively expressed ( Sheorey, 1994; Fairhurst, 2003; Li, et al., 2009) as

(2) {\sigma }_z^{\mathrm{0}}={\rho }_{a}g{H}_a+{\rho }_{c}gz,

(3) {\sigma }_x^{\mathrm{0}}={\sigma }_y^{\mathrm{0}}={\displaystyle \frac{\mu }{\mathrm{1}-\mu }}{\sigma }_z^{\mathrm{0}},

where H_a and ρ_a are, respectively, the thickness and the density of the alluvium layer; ρ_c is the average density of the rock mass; z is the vertical distance between the unit and the interface of the rock mass and alluvium layer; g is the acceleration of gravity; and μ is the Poisson’s ratio of the rock mass.

After the excavation of the coal seam, the equilibrium state of the original stress will be demolished and the stresses will be redistributed, i.e., the abutment pressure zones will form in the coal wall near the excavation boundary. The vertical stresses are much greater than the original ones in the abutment pressure zones. The peak abutment pressure zone is generally located near the excavation boundary ( Jacobi, 1966; Kulakov, 1975; Yavuz, 2004; Lai et al., 2006), and the concentration coefficient k stays constant in certain geological and mining conditions, generally k=1.5 ‒ 4.0 ( Yang, 1988; Xie et al., 1999).

In light of previous research, if there is no alluvium layer on the overburden rock, the ground surface subsidence on the excavation boundary can be expressed ( Yang, 1988) as

(4) w_{\mathrm{1}}\mathrm{\prime }={\displaystyle \frac{k-\mathrm{1}}{\mathrm{2}E}}{\rho }_{c}g{H}_c^{\mathrm{2}},

where E is the Young’s modulus of the rock mass; H_c and ρ_c are, respectively, the thickness and the average density of the rock mass.

Similarly, if there is extra stress (k‒1) ρ_agH_a from thick alluvium layer on the overburden rock mass (see Fig. 2), the compression deformation of the rock mass on the excavation boundary due to the extra stress can be expressed as

(4′) w_{\mathrm{1}}\mathrm{\prime }\mathrm{\prime }={\displaystyle \frac{k-\mathrm{1}}{E}}{\rho }_{a}g{H}_a{H}_c,

where E is the Young’s modulus of the rock mass; H_a and ρ_a are, respectively, the thickness and the density of the alluvium layer, and H_c is the thickness of the rock mass.

On the basis of Eqs. (4) and (4′), we can obtain the subsidence on the interface of the rock mass and alluvium layer above the excavation boundary in the case of a thick alluvium layer covering on the rock mass as

(5) w_{\mathrm{1}}=w_{\mathrm{1}}\mathrm{\prime }+w_{\mathrm{1}}\mathrm{\prime }\mathrm{\prime }={\displaystyle \frac{k-\mathrm{1}}{E}}\left({\rho }_{a}g{H}_a{H}_c+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\rho }_{c}g{H}_c^{\mathrm{2}}\right).

Deflection of rock stratum in bending zone

As mentioned above, the top trace of the fissure arch approximates a straight-line segment, which separates the bending zone from the fractured zone as shown in Fig. 3. The stresses within the fissure arch are elevated, while the stresses below are diminished ( Evans, 1941; Corlett and Emery, 1956; Panek, 1964; Wright, 1973). Therefore, it is feasible to estimate the supporting stresses from the fractured zone on the bending zone to be zero in the sub-critical mining stage. The stress state of rock strata in the bending zone is illustrated in Fig. 4.

Obviously, as compared with the thick alluvium layer, all rock strata belong in a hard stratum, which plays a key role in controlling ground surface subsidence. Accordingly, a thick hard rock stratum can be theoretically looked upon as a thin plate and the corresponding deformation theory for a thin plate can be used to study the ground surface subsidence ( Wang and Wu, 1995; Wu et al., 1997).

Furthermore, the goaf dimension is much greater than the thickness of the rock stratum in the bending zone and the deflection of the rock stratum is much less than its thickness, so it can thus be simplified as a rectangular elastic thin plate until it fractures and collapses ( Berry, 1960; Salamon, 1963; Li et al., 2006) as shown in Fig. 4 and Fig. 5.

In this section, we will examine the coordinate system o_{\mathrm{1}}\xi \eta in Fig. 5. The dimensions of the rock plate are, respectively, L₁ and L₂ in the directions of strike and inclination; the thickness of the rock plate, namely, the thickness of the bending zone, is h. In terms of our assumption, the load intensity on the surface of the rock plate from the thick alluvium layer is considered as the uniformly distributed loads, i.e., q₀ = ρ_agH_a. Thus, the displacement differential equations is as follows (Timoshenko and Woinowsky-Krieger, 1959; Xu, 1992; Reddy, 2006):

(6) {\nabla }^{\mathrm{2}}{\nabla }^{\mathrm{2}}{w}_{\mathrm{2}}\left(\xi ,\eta \right)={\displaystyle \frac{q}{D}},

where {\nabla }^{\mathrm{2}} is the Laplace Operator, {\nabla }^{\mathrm{2}}={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial {\xi }^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial {\eta }^{\mathrm{2}}}} ; D is the bending rigidity of the rock plate, D={\displaystyle \frac{E{h}^{\mathrm{3}}}{\mathrm{12}\left(\mathrm{1}-{\mu }^{\mathrm{2}}\right)}} ; E is Young’s modulus of the rock plate; μ is Poisson’s ratio of the rock plate; and q is the total load acting on the rock plate. The total load consists of two parts: the uniformly distributed overburden alluvium load q₀ and the self-weight of the rock plate ρ_cgh, i.e., q= q₀+ρ_cgh.

Under the simple supporting conditions, i.e., when ξ=0 or L₁, then w₂ = 0 and {\displaystyle \frac{{\partial }^{\mathrm{2}}{w}_{\mathrm{2}}}{\partial {\xi }^{\mathrm{2}}}}=\mathrm{0} ; and when η=0 or L₂, then w₂=0 and {\displaystyle \frac{{\partial }^{\mathrm{2}}{w}_{\mathrm{2}}}{\partial {\eta }^{\mathrm{2}}}}=\mathrm{0} , we may obtain the deflection of the rock plate according to the Navier solution:

(7) w_{\mathrm{2}}=f\left(\xi ,\eta \right)={\displaystyle \frac{\mathrm{16}q}{{\mathrm{\pi }}^{\mathrm{6}}D}}{\displaystyle \sum _{i=\mathrm{1},\mathrm{3},\mathrm{5}}^{\infty }{\displaystyle \sum _{j=\mathrm{1},\mathrm{3},\mathrm{5}}^{\infty }{\displaystyle \frac{\sin {\displaystyle \frac{i\mathrm{\pi }\xi }{L_{\mathrm{1}}}}\sin {\displaystyle \frac{j\mathrm{\pi }\eta }{L_{\mathrm{2}}}}}{ij{\left({\displaystyle \frac{i^{\mathrm{2}}}{L_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{j^{\mathrm{2}}}{L_{\mathrm{2}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}}}.}}

This is a rapidly converging series, and a satisfactory approximation is obtained by taking only the first term of the series, namely,

(7′) w_{\mathrm{2}}\approx f\left(\xi ,\eta \right)={\displaystyle \frac{\mathrm{16}q}{{\mathrm{\pi }}^{\mathrm{6}}D}}{\left({\displaystyle \frac{L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}}{L_{\mathrm{1}}^{\mathrm{2}}+L_{\mathrm{2}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}.

This result is about 2.5 percent in error (Timoshenko, Woinowsky-Krieger, 1959; Xu, 1992; Reddy, 2006). For simplification, let

(8) w_{m_{\mathrm{2}}}={\displaystyle \frac{\mathrm{16}q}{{\mathrm{\pi }}^{\mathrm{6}}D}}{\left({\displaystyle \frac{L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}}{L_{\mathrm{1}}^{\mathrm{2}}+L_{\mathrm{2}}^{\mathrm{2}}}}\right)}^{\mathrm{2}},

then

(9) w_{\mathrm{2}}=f\left(\xi ,\eta \right)=w_{m_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}.

From Eq. (9), it may be seen that maximum deflection occurs at the center of the plate \left({\displaystyle \frac{L_{\mathrm{1}}}{\mathrm{2}}},{\displaystyle \frac{L_{\mathrm{2}}}{\mathrm{2}}}\right) .

According to the classical theory , the applicability of the elastic thin plate is as follows:

1) The deflection of the rock plate w₂ is far less than the thickness of thin plate h;

2) The ratio of the plate thickness h to the minimum dimension of the thin plate ranges from 1/5 to 1/100.

In practical underground mining engineering, the applicability can be liberalized to some extent ( Wu and Wang, 1994; Ma, 2001; Chang and Wang, 2003; Li et al., 2006, Wang et al., 2008).

Total subsidence on the interface between rock mass and alluvium layer

According to the Eqs. (1), (7′), and (9), the total subsidence in the bending zone w_b is easily obtained:

(10) w_b=w_b\left(\xi ,\eta \right)=w_{\mathrm{1}}+w_{m_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}},

Because there is almost no separation or crack in the bending zone, it is reasonable to assume the subsidence of the rock in the bending zone as the subsidence on the interface between the rock mass and alluvium layer.

Derivation of the ground surface subsidence

Ground surface subsidence basin

As previously discussed, the ground surface subsidence is mainly induced by the subsidence of the bending zone above the excavation zone, and the subsidence on the interface between the rock mass and alluvium layer can be considered as a kind of “excavation,” which will lead to ground surface subsidence through the thick alluvium layer as shown in Fig. 6. That is to say, the subsidence on the ground surface is directly attributed to the subsidence on the interface between the rock mass and the alluvium layer.

To derive the ground surface subsidence basin, two coordinate systems must be established: oxy and o'ξη (Fig. 7). The oxy coordinate system is on the ground surface with x-axis and y-axis respectively pointing to the directions of the strike and inclination; meanwhile, the o'ξη coordinate system is on the interface between the rock mass and the alluvium with ξ-axis and η-axis respectively pointing to the directions of the strike and inclination.

Based on the stochastic medium theory, the probabilistic integral method is well suited to describe the ground surface subsidence and other deformations due to the excavation under the thick alluvium layer, which is regarded as random medium. According to the probabilistic integral method, we have the ground surface subsidence basin ( Litwiniszyn, 1974; Baxevanis and Plexousakis, 2002):

(11) w=w\left(x,y\right)={\displaystyle \underset{\sigma }{\iint }{w}_b\left(\xi ,\eta \right)w_e\left(x-\xi ,y-\eta \right)\mathrm{d}\sigma },

where w_b (ξ, η) is the excavation under the thick alluvium layer; and w_e (x‒ξ, y‒η) is the micro-unit mining subsidence basin, which is expressed by

(12) w_e\left(x-\xi ,y-\eta \right)={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}\left[{\left(x-\xi \right)}^{\mathrm{2}}+{\left(y-\eta \right)}^{\mathrm{2}}\right]},

where r is the major influence radius of the thick alluvium layer.

In Eq. (10), we have obtained the subsidence on the interface between the rock mass and thick alluvium layer w_b (ξ, η), which is considered as a type of “excavation.” By substituting Eqs. (10) and (12) for Eq. (11), the ground surface subsidence basin expression could be considered as

w=w\left(x,y\right)

(13) w=w\left(x,y\right)={\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\left[w_{\mathrm{1}}+w_{m_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{1}}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}\right]{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}}}{\displaystyle {\displaystyle {\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}\left[{\left(x-\xi \right)}^{\mathrm{2}}+{\left(y-\eta \right)}^{\mathrm{2}}\right]}\mathrm{d}\xi \mathrm{d}\eta ,}}

i.e.,

(14) w={\displaystyle \frac{w_{\mathrm{1}}}{r^{\mathrm{2}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi }{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}}}\mathrm{d}\eta +{\displaystyle \frac{w_{m_{\mathrm{2}}}}{r^{\mathrm{2}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi {\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta .

It is clearly impossible to retrieve the explicit primary functions in Eq. (14). However, it is very convenient to calculate the integral quantities by means of numerical integration.

Ground surface subsidence and other deformations in principal profiles

Ground surface subsidence in principal profiles

According to Eq. (14), a ground surface subsidence formula can easily be obtained in principal profiles. For the strike principal profile:

(15) w\left(x\right)={w\left(x,y\right)|}_{y=\frac{L_{\mathrm{2}}}{\mathrm{2}}}={\displaystyle \frac{w_{\mathrm{1}}}{r^{\mathrm{2}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{2}}}{\mathrm{2}}-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta {\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi }}+{\displaystyle \frac{w_{m\mathrm{2}}}{r^{\mathrm{2}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{2}}}{\mathrm{2}}-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta {\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi .

Now, let C_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{2}}}{\mathrm{2}}-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta } and C_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{2}}}{\mathrm{2}}-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta .

Hence,

(16) w\left(x\right)={\displaystyle \frac{C_{\mathrm{1}}{w}_{\mathrm{1}}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi }+{\displaystyle \frac{C_{\mathrm{2}}{w}_{m_{\mathrm{2}}}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi .

Obviously, \mathrm{0}\le C_{\mathrm{1}},C_{\mathrm{2}}\le \mathrm{1} .

Similarly, we can also obtain the ground surface subsidence in the inclination principal profile:

(17) w\left(y\right)={w\left(x,y\right)|}_{x=\frac{L_{\mathrm{1}}}{\mathrm{2}}}={\displaystyle \frac{D_{\mathrm{1}}{w}_{\mathrm{1}}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}{\mathrm{e}}^{\begin{array}{l}-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}\hfill \end{array}}}\mathrm{d}\eta +{\displaystyle \frac{D_{\mathrm{2}}{w}_{m_{\mathrm{2}}}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_y}}}{e}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta ,

where D_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{1}}}{\mathrm{2}}-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi } and D_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(\frac{L_{\mathrm{1}}}{\mathrm{2}}-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi , \mathrm{0}\le D_{\mathrm{1}},D_{\mathrm{2}}\le \mathrm{1} .

Other deformations in principal profiles

According to the basic theory of ground movement, i.e., functional relations exist between ground surface subsidence and the other deformations, thus the other deformations can be derived on the basis of the ground surface subsidence expression. As is well known, the slope and curvature are, respectively, the first order and the second order derivatives of the corresponding subsidence, thus it is not difficult to obtain the slope and curvature expressions in the strike principal profile:

(18) I\left(x\right)={\displaystyle \frac{\partial w\left(x\right)}{\partial x}}={\displaystyle \frac{C_{\mathrm{1}}{w}_{\mathrm{1}}}{r}}\left({\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{x}^{\mathrm{2}}}-{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-L_{\mathrm{1}}\right)}^{\mathrm{2}}}\right)-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{C}_{\mathrm{2}}{w}_{m_{\mathrm{2}}}}{r^{\mathrm{3}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\left(x-\xi \right)\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\mathrm{d}\xi ,

(19) k\left(x\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}w\left(x\right)}{\partial x^{\mathrm{2}}}}=-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{C}_{\mathrm{1}}{w}_{\mathrm{1}}}{r^{\mathrm{3}}}}\left[x{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{x}^{\mathrm{2}}}-\left(x-L_{\mathrm{1}}\right){\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-L_{\mathrm{1}}\right)}^{\mathrm{2}}}\right]-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{C}_{\mathrm{2}}{w}_{m_{\mathrm{2}}}}{r^{\mathrm{3}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{1}}}\sin {\displaystyle \frac{\mathrm{\pi }\xi }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(x-\xi \right)}^{\mathrm{2}}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{\left(x-\xi \right)}^{\mathrm{2}}}{r^{\mathrm{2}}}}\right]\mathrm{d}\xi .

In view of Awershen’s assumption that horizontal displacement and torsion are respectively proportional to the slope and curvature ( Awershen, 1959; Liu and Liao, 1965; National Coal Board Mining Department, 1975), it is relatively simple to write down their expressions:

(20) U\left(x\right)=b\cdot I\left(x\right),

(21) ϵ\left(x\right)=b\cdot k\left(x\right),

where b is the coefficient of lateral displacement.

Similarly, the expressions of slope, curvature, horizontal displacement and torsion in the inclination principal profile can simply be defined as

(22) I\left(y\right)={\displaystyle \frac{\partial w\left(y\right)}{\partial y}}={\displaystyle \frac{D_{\mathrm{1}}{w}_{\mathrm{1}}}{r}}\left[{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{y}^{\mathrm{2}}}-{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-L_{\mathrm{2}}\right)}^{\mathrm{2}}}\right]-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{D}_{\mathrm{2}}{w}_m}{r^{\mathrm{3}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\left(y-\eta \right)\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{1}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}}\mathrm{d}\eta ,

(23) k\left(y\right)={\displaystyle \frac{\partial w^{\mathrm{2}}\left(y\right)}{\partial y^{\mathrm{2}}}}=-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{D}_{\mathrm{1}}{w}_{\mathrm{1}}}{r^{\mathrm{3}}}}\left[y{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{y}^{\mathrm{2}}}-\left(y-L_{\mathrm{2}}\right){\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-L_{\mathrm{2}}\right)}^{\mathrm{2}}}\right]-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{D}_{\mathrm{2}}{w}_m}{r^{\mathrm{3}}}}{\displaystyle {\int }_{\mathrm{0}}^{L_{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\pi }\eta }{L_{\mathrm{2}}}}}{\mathrm{e}}^{-\frac{\mathrm{\pi }}{r^{\mathrm{2}}}{\left(y-\eta \right)}^{\mathrm{2}}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{2}\mathrm{\pi }{\left(y-\eta \right)}^{\mathrm{2}}}{r^{\mathrm{2}}}}\right]\mathrm{d}\eta ,

(24) U\left(y\right)=b\cdot I\left(y\right),

(25) ϵ\left(y\right)=b\cdot k\left(y\right).

Similar to the conventional prediction models, the related parameters in this prediction model, such as the coefficient of lateral displacement, the major influence radius, the concentration pressure coefficient, the mining dimensions, the extracted seam thickness, the mining depth, and the thickness of alluvium layer, are determined by the actual observed data of the local mining fields. In addition, the mechanical parameters for the rock mass, including the average density, the Poisson’s ratio, and the Young’s modulus, depend on the properties of the actual overburdened rock mass.

Finally, as previously mentioned, it is impossible to retrieve the explicit primary functions in the above equations. In general, the numerical integration methods are used to calculate the predicted ground surface subsidence and other deformations. Fortunately, numerous software systems have now been developed which can provide an integrated development environment for various numerical computations. Matlab programming was used to calculate the predicted results stated in this paper.

A prediction case

To evaluate the above prediction model with actual observed data, a ground surface subsidence observation line was established at Dongpang Coal Mine, which is located north of Xingtai, Hebei Province in China, at a distance of about 19 km (see Fig. 8).

For the purposes of surveying the ground surface subsidence and other deformations on the strike principal profile, we established the ground surface subsidence observation line with 23 measured points (including 2 controlling points and 21 observation points) on the strike principal of working panel 2107 (see Fig. 9). In view of the actual mining depth, the separation distances between the observation points were kept around 25 meters, and the 2 controlling points (R₁, R₂) were beyond the influence area of the underground mining. The observation precision was determined based on the China Coal Mines survey handbook, so that the actual observations were carried out precisely.

The land surface in this area is relatively plain at an elevation of about 95 m above the mean Huanghai Sea level. The main strata belong to the Ordovician system, and the related geological and mining parameters are as follows:

The average inclination angle of the coal seam α=6.5°, which approximates horizontal mining; the extracted seam averages 4.0 m in thickness, i.e., m= 4.0 m; the average mining depth H = 257 m, the average thickness of alluvium layer H_a= 140 m (including clay, sandy clay, and gravel), the density of alluvium layer d_a=17,000 N/m³, the overburden rock strata mainly consist of shale and limestone, the average thickness of rock mass H_c=117 m, the average density of rock mass d_a =26,000 N/m³, the Poisson’s ratio: μ=0.22, the Young’s modulus E=24.0 GPa, which is supposed to be half the value of the rock samples in the laboratory ( Lin, 1984; Boreek and Chudek, 1985; Mohammad et al., 1997); the mining dimensions are 139 m along strike and 160 m along dip, i.e., L₁=139 m, L₂=160 m; the lateral movement factor b=0.44, the tangent of the main influence angle tanβ=1.75, the concentration pressure coefficient: k=3.5; the thickness of the bending zone, i.e., the thickness of the plate h, was determined by using empirical formula ( Mining Institute of China, 1981; Bai, et al. 1995; Du and Weng, 1997; Zhang and Shen, 2004):

(26) h=H_{\mathrm{c}}-{\displaystyle \frac{m}{a_{\mathrm{1}}m+a_{\mathrm{2}}}},

where H_c is the thickness of the rock mass; m is the extracted seam thickness; a₁ and a₂ are respectively the coefficients depending on the rock strata lithology. According to the actual parameters of local mine a₁=0.008 and a₂=0.0148, we can obtain h= 31.5 m.

Since there are 18 observed points inside the ground surface subsidence basin, it is significant to compare the predicted and actual observed results of these observed points. The predicted subsidence and horizontal displacement values, as well as the corresponding actual observed values, are respectively shown in Figs.10(a) and 10(b) with ξ-axis being the horizontal axis.

To evaluate the prediction accuracy, we make the assumptions that there are no observation errors and the prediction errors are approximately equal for all points. Thus, the mean square errors of ground surface subsidence and horizontal displacement between the predicted and the actual observed values can be calculated by using the Bessel formula:

(27) m_w=\pm \sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n{\left(W{p}_i-W{o}_i\right)}^{\mathrm{2}}}}{n-\mathrm{1}}}},

(27′) m_u=\pm \sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n{\left(U{p}_i-U{o}_i\right)}^{\mathrm{2}}}}{n-\mathrm{1}}}},

where subscript i is the mark of the actual observed points; n is the total number of the actual observed points, now n =18; Wp_i and Wo_i are, respectively, values of the predicted and actual observed subsidence; and Up_i and Uo_i are, respectively, values of the predicted and actual observed horizontal displacement.

Finally, we can yield m_w≈±0.150 m, and m_u≈±0.141 m. Obviously, the predicted values for the subsidence are in good agreement with the corresponding actual observed values. The horizontal displacement prediction mean square errors are a somewhat large, which is probably due to the movement factor b of the local mine field or other influence factors. In view of the complexity of the geological and mining conditions, the predicted results can be acceptable.

Conclusions

1) Because there are tremendous differences between the mechanical behaviors of the rock mass and alluvium layer, the mechanisms of displacements and deformation for rock mass and alluvium layer are greatly distinguishable from one another, thus it is not reasonable to apply the conventional prediction models to ground surface subsidence in the geological condition of thick alluvium layer.

2) The different mechanisms of the displacements and deformation between rock mass and alluvium layer have been respectively considered and described in this prediction model, making it possible to significantly improve the prediction accuracy of the ground surface subsidence and other deformations.

3) The proposed prediction model is only suitable for sub-critical mining subsidence in the geological condition of thick alluvium layer, thus that it cannot replace the conventional prediction models in other geological and mining conditions.

References

[1] Álvarez-Fernández M I, González-Nicieza C, Menéndez-Díaz A, Álvarez-Vigil A E (2005). Generalization of the n-k influence function to predict mining subsidence. Eng Geol, 80(1–2): 1–36

[2] Awershen C R (1959). Rock Movement due to Underground Mining in Coal Mine. Beijing: Coal Industry Press, 24–32 (in Chinese)

[3] Bai M, Kendorski F, Van Roosendaal D (1995). Chinese and North American high-extraction underground coal mining strata behavior and water protection experience and guidelines. Proceedings of the 14th International Conference on Ground Control in Mining, Morgantown, 209–217

[4] Baxevanis Th, Plexousakis M (2009). Estimation of base settlement from the surface subsidence profile: two-Dimensional plain field of displacements. Int J Numer Anal Methods Geomech, 33(8): 1109–1121

[5] Berry D S (1960). An elastic treatment of ground movement due to mining - I. Isotropic ground. J Mech Phys Solids, 8(4): 280–292

[6] Boreek M, Chudek M (1985). Rock Mass Mechanics in Mining. Beijing: Coal Industry Press (in Chinese)

[7] Chang Z Q, Wang J Z (2003). The space problems research under thick loose layer bend sinking. Journal Mine Surveying, 9: 36–38 (in Chinese)

[8] Corlett A V, Emery C L (1959). Prestress and stress redistribution in rocks around a mine opening. Bull Can Min Metall., 372–383

[9] Cui X M, Miao X X, Wang J A, Yang S, Liu H, Song Y, Liu H, Hu X (2000). Improved prediction of differential subsidence caused by underground mining. Int J Rock Mech Min Sci, 37(4): 615–627

[10] Dai H Y, Wang J Z, Cai M F (2002). Seam dip angle based mining subsidence model and its application. Int J Rock Mech Min Sci, 39(1): 115–123

[11] Du S G, Weng X H (1997). Relationship between dip angle of coal seam and heights of Three Zones. Journal of Engineering Geology, 5(3): 211–217 (in Chinese)

[12] Evans W H (1941). The strength of undermined strata. Trans Inst Min Metall, 50: 475–500

[13] Fairhurst C (2003). Stress estimation in rock: a brief history and review. Int J Rock Mech Min Sci, 40(7-8): 957–973

[14] He G Q, Yang L; Research on Mine Subsidence (1991). Xuzhou: China University of Mining Technology Press (in Chinese)

[15] Jacobi O (1966). The increase of roof flaking in longwall faces as a result of working under pillar edges and of abutment pressure of adjacent workings. Int J Rock Mech Min Sci Geomech Abstr, 3(3): 221–230

[16] Knothe S (1957). Observations of surface movements and their theoretical interpretation, Proceedings of the European congresson ground movement, University of Leeds, UK, 210–215

[17] Kratzsch H (1983). Mining Subsidence Engineering. Berlin: Springer-Verlag

[18] Kulakov G S (1975). The stress distribution in the abutment pressure zone according to the results of field measurements. Int J Rock Mech Min Sci Geomech Abstr, 12(10): 231–234

[19] Lai X P, Cai M F, Ren F H, Xie M W, Esaki T (2006). Assessment of rock mass characteristics and the excavation disturbed zone in the Lingxin Coal Mine beneath the Xitian river, China. Int J Rock Mech Min Sci, 43(4): 572–581

[20] Li G, Mizuta Y, Ishida T, Li H, Nakama S, Sato T (2009). Stress field determination from local stress measurements by numerical modeling. Int J Rock Mech Min Sci, 46(1): 138–147

[21] Li W X, Liang X L, Zhao S T, Mei S H (2006). The elastic-thin-plate bending method for predicting ground subsidence due to underground mining. Engineering Mechanics., 23(8): 177–181 (in Chinese)

[22] Liao C P (1993). Fuzzy infuence function method for calculating mine subsidence in a horizontal seam. Geotech Geol Eng, 11(4): 235–247

[23] Lin S, Reddish D J, Whittaker B N (1992). An integrated analytical model of subsidence induced by level seam extractions. Geotech Geol Eng, 10(3): 203–221

[24] Lin Y M (1984). Simulation Research on Experimental Rock Mechanics. Beijing: Coal Industry Press (in Chinese)

[25] Litwiniszyn J (1957). The theories and model research of movements of ground masses. Proceedings of the European Congress on Ground Movement, University of Leeds, 206–209

[26] Litwiniszyn J (1974). Stochastic Methods in the Mechanics of Granular Bodies. Wien: Springer-Verlag

[27] Liu B S, Liao G H (1965). Basic Law Governing Ground Movement in Coal Mines. Beijing: Coal Industry Press, 55–60 (in Chinese)

[28] Ma D (2001). Discussions about the relations of three-dimension elastic theory, classical thin plate theory and thick-walled construction theory. Journal of Liaoning Institute of Technology, 21(1): 48–50 (in Chinese)

[29] Mining Institute of China (1981). Rock Strata and Ground Surface Displacement in Coal Mines, Beijing: Coal Industrial press (in Chinese)

[30] Mohammad N, Reddish D J, Stace L R (1997). The relation between in situ and laboratory rock properties used in numerical modelling. Int J Rock Mech Min Sci, 34(2): 289–297

[31] National Coal Board Mining Department (1975). UK Subsidence Engineer’s Handbook, NCB Mining Department, London

[32] Oravecz K I (1986). Improved prediction of surface subsidence using the influence function approach. In: Sangorm Symposium: The Effect of Underground Mining on Surface, 73–80

[33] Panek, L A (1964). Design of bolting systems to reinforce bedded mine roof. U.S. Bureau of Mines Report of Investigations5155, 16

[34] Peng S S, Chiang H S (1984). Long Wall Mining, New York: Wiley, 707–708

[35] Reddy J N (2006). Theory and Analysis of Elastic Plates and Shells (2nd edition). Boca Raton: CRC Press

[36] Ren G, Reddish D J, Whittaker B N, Whittaker B N (1987) Mining subsidence and displacement prediction using influence function methods. Min Sci Technol, 5(1): 89–104

[37] Salamon M D G (1963). Elastic analysis of displacements and stresses induced by the mining of seam or reef deposits–Part 1. Fundamental principles and basic solutions as derived from idealized models. J S Afr Inst Min Metall, 64: 128–149

[38] Sansone E C, Ayres Da Silva L A (1998). Numerical modeling of the pressure arch in underground mines. Int J Rock Mech Min Sci, 35(4–5): 436

[39] Sheorey P R (1994). A theory for in situ stresses in isotropic and transversely isotropic rock. Int J Rock Mech Min Sci Geomech Abstr, 31(1): 23–34

[40] Sheorey P R, Lou J P, Singh K B, Singh S K (2000). Ground subsidence observations and a modified influence function method for complete subsidence prediction. Int J Rock Mech Min Sci, 37(5): 801–818

[41] Su Z J, Liu W S, Yang L (2003). A mechanical model of subsidence due to mining and its application. Mech Eng, 25(2): 45–47

[42] Timoshenko S, Woinowsky-Krieger S (1959). Theory of Plates and Shells. New York: Mcgraw-Hill Book Company, Inc.

[43] Wang C G, Song Z Q, Shi Y K (2004). Study on the relation between stratum movement and subsidence of flat seam mining. Rock and Soil Mechanics, 25(8): 1343–1346 (in Chinese)

[44] Wang J A, Shang X C, Ma H T (2008). Investigation of catastrophic ground collapse in Xingtai gypsum mines in China. Int J Rock Mech Min Sci, 45(8): 1480–1499

[45] Wang J Z, Li Y S, Zhou X (1997). Ground movement caused by mining under thick alluvium. Journal of China Coal Society, 22(1): 90–93 (in Chinese)

[46] Wang J Z, Wu L X (1995). Study of strata and surface movement and mining under village houses in the mountainous area in Furong coal field. Journal of China University of Mining and technology, 24(2): 52–58 (in Chinese)

[47] Wright F D (1973). Roof control through beam action and arching. In: Cummins A B, Given I A, eds. SME Mining Engineering Handbook, New York: Society of Mining Engineers, 13: 80–96

[48] Wu L X, Qian M G, Wang J Z (1997). The influence of a thick hard rock stratum on underground mining subsidence. Int J Rock Mech Min Sci, 34(2): 341–344

[49] Wu L X, Wang J Z (1994). Study of deformation model of a controlling holding plate when large area is extracted continuously. Journal of China Coal Society, 19(3): 233–242 (in Chinese)

[50] Xie H P, Chen Z H, Wang J C (1999). Three-dimensional numerical analysis of deformation and failure during top coal caving. Int J Rock Mech Min Sci, 36(5): 651–658

[51] Xu Z L (1992). Mechanics of Elasticity (II). 3rd edition. Beijing: Higher Education Press (in Chinese)

[52] Yang L (1988). The new knowledge about the mechanism and the law of displacement in rock and ground due to coal mining. Journal of Fuxin Mining Institute, 7(1): 9–18 (in Chinese)

[53] Yavuz H (2004). An estimation method for cover pressure re-establishment distance and pressure distribution in the goaf of long wall coal mines. Int J Rock Mech Min Sci, 41(2): 193–205

[54] Zhang J C, Shen B H (2004). Coal mining under aquifers in China: a case study. Int J Rock Mech Min Sci, 41(4): 629–639

[55] Zhang X D, Zhao Y H, Liu S J (1999). A new method of calculating surface subsidence and deformations under thick alluvial soil. Chinese Journal of Nonferrous Metals, 9(2): 435–440 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg