# Time Function Model of Surface Subsidence Based on Inversion Analysis in Deep Soil Strata

Authors: Jihuan Han; Chenchen Hu; Jiuqun Zou
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

As a common geological disaster, surface subsidence caused by mining underground resources has always been a hot and difficult topic in the civil engineering field. Aimed at the shortcomings of existing time function models in predicting mining subsidence in deep soil strata, a more accurate and reasonable time function model, called the composite function model, was established based on an inverted analysis of measured data. The results showed that the composite function model could describe the whole subsidence process of a deep soil surface and agreed well with the measured data. The model parameters were calculated by specific formulas, which improved the reliability of the subsidence prediction results under different mining conditions. The new model provided important guiding significance for preventing subsidence geological disasters and determining the coal mining time under the buildings, the railways, and the water bodies in deep soil strata.


---

## Body

## 1. Introduction

Energy and mineral resources play a key role in national economic development around the world [1]. With the implementation of the deep coal mining strategy, since 2002, China has built 71 shafts with soil thicknesses over 400 m. Mining subsidence not only causes great waste of land resources but also seriously affects the safety of underground structures and surface buildings [2, 3]. The geological hazards caused by mining subsidence are particularly serious in deep soil strata [4, 5]. Therefore, a time function model must be established to describe the temporal and spatial variations of surface subsidence. At present, the commonly used time function models include the Knothe model [6], Kowalski model [7], Sroka–Schober model [8], revised Hrries model, [9] and improved Knothe model [10]. However, neither the Knothe model nor the Kowalski model can describe the whole process of surface subsidence. The parameter values in the Sroka–Schober model and the revised Hrries model are difficult to determine, resulting in more uncontrollable factors. The improved Knothe model is inaccurate in predicting the deep soil surface subsidence. Considering the shortcomings of the abovementioned models, a new time function model was proposed by inverting and analysing the measured data of deep soil surface subsidence. Compared with the improved Knothe model, the new model is more accurate and reasonable in predicting the deep soil surface subsidence.


## 2. Typical Time Function Models

### 2.1. Regular Analysis

When the coal seam is mined out, the overlying rock strata in goaf will move due to the destruction of stress balance, and its moving range will increase with the expansion of the goaf area. Once the working face has a certain advance distance, coal mining affects the surface, causing the surface to move and form a subsidence basin. According to the subsidence process of the actual surface point, some researchers believe that the subsidence, subsidence velocity, and subsidence acceleration are all 0 at the initial time [6–11]. The regular laws of subsidence, subsidence velocity, and subsidence acceleration with time can be described as follows: subsidence: 0⟶+max; subsidence velocity: 0⟶+max⟶0; and subsidence acceleration: 0⟶+max⟶0⟶-max⟶0.Surface subsidence is a complex process with time and space. The whole process can be regarded as an approximately inverse S-shaped curve with time, as shown in Figure1 [6–12]. Thus, the whole process can be roughly divided into four stages on the basis of velocity variation: the initial stage (I), accelerated stage (II), decelerated stage (III), and stable stage (IV). At stage I, the subsidence velocity is very low and the subsidence increases slowly from 0. At stage II, the subsidence velocity increases to the extremum with the subsidence increasing rapidly. At stage III, the subsidence velocity decreases while the subsidence continues to increase. At stage IV, the subsidence velocity decreases to 0 and the subsidence increases to a stable value. When the coal seam is mined under the conditions of thick soil and thin rock strata, due to the gravity of the soil strata, the overlying rock strata will have large bending deformation in a short time. In other words, the duration of stage I is relatively short in deep soil strata.Figure 1
The whole process of surface subsidence.


### 2.2. Model Analysis

#### 2.2.1. Knothe Model

Assuming that the subsidence velocity at a certain moment is directly proportional to the difference between the dynamic subsidence value at this moment and the maximum subsidence value, Knothe [6] obtained a differential function as follows:(1)dwtdt=cw0−wt,where wt is the instantaneous subsidence value of a point on the surface at time t, w0 is the maximum subsidence value of the point, and c is a time coefficient related to lithology.With the initial conditionw(0) = 0, (1) can be derived as(2)wt=w01−e−ct.The Knothe time function can be expressed as(3)Tt=1−e−ct,where Tt is the dimensionless instantaneous subsidence of a point on the surface at time t, Tt=wt/w0.


#### 2.2.2. Kowalski Model

Considering the phenomenon of surface subsidence lagging behind underground mining, Kowalski [7] rewrote the Knothe model as(4)Tt=qt1−Ae−ct−τop,qt=0,t<τop,1,t≥τop,where A is a model parameter and τop is the time from initial mining to the appearance of surface subsidence.Obviously, the Kowalski model still cannot directly reflect the variation law of subsidence velocity.


#### 2.2.3. Sroka–Schober Model

Considering the convergence and compaction of porous rock strata, Sroka et al. [8] proposed a time function with two parameters, which can be defined by(5)Tt=1+ξf−ξe−ft−ff−ξe−ξt,where f is the relative convergence rate of rock strata, f = 0.01 a−1 means that the annual convergence is 1% of the original volume, ξ is the time coefficient of the overlying rock strata, when the roof is managed by caving method, and the value of ξ ranges from 20 to 70 a−1.The Sroka–Schober model is used for the prediction of rock mass subsidence, while it can also calculate surface subsidence by adjusting the model parameter values. For example, whenξ = 20 a−1, f = 2 a−1, the dimensionless subsidence curve is calculated by the Sroka–Schober model, as shown in Figure 2. It shows that the Sroka–Schober model is superior to the Knothe model, especially at the initial stage. However, it is difficult to obtain the actual values of f and ξ.Figure 2
Dimensionless subsidence curves calculated by the Sroka–Schober model and Knothe Model.


#### 2.2.4. Revised Hrries Model

In view of the deficiencies of the Gompertz model and Logistic model in predicting foundation subsidence, Liu et al. [9] established a time function by revising the Harris model, which can be derived as(6)Tt=1−11+w0mtn,where m and n are the model parameters with positive values.The revised Harris model can describe the whole process of surface subsidence. However, the parametersm and n solved by fitting cannot reflect the stratigraphic characteristics. That is, it is difficult to determine which geological and mining factors control these parameters, and this reduces the reliability of prediction results. In addition, (6) contains the parameter w0, which indicates that the model values are affected by the maximum subsidence value. This is not a time function in the complete sense.


#### 2.2.5. Improved Knothe Model

Regarding the surface subsidence curve as a symmetrical distribution, Chang and Wang [10] improved the Knothe model as follows:(7)Tt=T1t=0.5e−cτ−t−e−cτ0<t≤τ,T2t=0.52−e−ct−τ−e−cττ<t≤T0,where T0 is the total subsidence time and τ is the time required for the subsidence velocity to reach the maximum, τ=0.5T0.In fact, the surface subsidence curve is an asymmetrical distribution in deep soil strata. Taking the mining of the 1308 working face of Guotun Coal Mine in Shandong Province as an example, Z1, Z2, and Z3 are the surface observation points along the strike. The measured subsidence [13] and the theoretical subsidence calculated by (7) are shown in Figure 3. Obviously, there is a large deviation between the theoretical and measured values, especially at the accelerated and decelerated stages. Therefore, the improved Knothe model cannot accurately and reasonably predict deep soil surface subsidence.Figure 3
Comparison between the measured and theoretical subsidence.


## 2.1. Regular Analysis

When the coal seam is mined out, the overlying rock strata in goaf will move due to the destruction of stress balance, and its moving range will increase with the expansion of the goaf area. Once the working face has a certain advance distance, coal mining affects the surface, causing the surface to move and form a subsidence basin. According to the subsidence process of the actual surface point, some researchers believe that the subsidence, subsidence velocity, and subsidence acceleration are all 0 at the initial time [6–11]. The regular laws of subsidence, subsidence velocity, and subsidence acceleration with time can be described as follows: subsidence: 0⟶+max; subsidence velocity: 0⟶+max⟶0; and subsidence acceleration: 0⟶+max⟶0⟶-max⟶0.Surface subsidence is a complex process with time and space. The whole process can be regarded as an approximately inverse S-shaped curve with time, as shown in Figure1 [6–12]. Thus, the whole process can be roughly divided into four stages on the basis of velocity variation: the initial stage (I), accelerated stage (II), decelerated stage (III), and stable stage (IV). At stage I, the subsidence velocity is very low and the subsidence increases slowly from 0. At stage II, the subsidence velocity increases to the extremum with the subsidence increasing rapidly. At stage III, the subsidence velocity decreases while the subsidence continues to increase. At stage IV, the subsidence velocity decreases to 0 and the subsidence increases to a stable value. When the coal seam is mined under the conditions of thick soil and thin rock strata, due to the gravity of the soil strata, the overlying rock strata will have large bending deformation in a short time. In other words, the duration of stage I is relatively short in deep soil strata.Figure 1
The whole process of surface subsidence.


## 2.2. Model Analysis

### 2.2.1. Knothe Model

Assuming that the subsidence velocity at a certain moment is directly proportional to the difference between the dynamic subsidence value at this moment and the maximum subsidence value, Knothe [6] obtained a differential function as follows:(1)dwtdt=cw0−wt,where wt is the instantaneous subsidence value of a point on the surface at time t, w0 is the maximum subsidence value of the point, and c is a time coefficient related to lithology.With the initial conditionw(0) = 0, (1) can be derived as(2)wt=w01−e−ct.The Knothe time function can be expressed as(3)Tt=1−e−ct,where Tt is the dimensionless instantaneous subsidence of a point on the surface at time t, Tt=wt/w0.


### 2.2.2. Kowalski Model

Considering the phenomenon of surface subsidence lagging behind underground mining, Kowalski [7] rewrote the Knothe model as(4)Tt=qt1−Ae−ct−τop,qt=0,t<τop,1,t≥τop,where A is a model parameter and τop is the time from initial mining to the appearance of surface subsidence.Obviously, the Kowalski model still cannot directly reflect the variation law of subsidence velocity.


### 2.2.3. Sroka–Schober Model

Considering the convergence and compaction of porous rock strata, Sroka et al. [8] proposed a time function with two parameters, which can be defined by(5)Tt=1+ξf−ξe−ft−ff−ξe−ξt,where f is the relative convergence rate of rock strata, f = 0.01 a−1 means that the annual convergence is 1% of the original volume, ξ is the time coefficient of the overlying rock strata, when the roof is managed by caving method, and the value of ξ ranges from 20 to 70 a−1.The Sroka–Schober model is used for the prediction of rock mass subsidence, while it can also calculate surface subsidence by adjusting the model parameter values. For example, whenξ = 20 a−1, f = 2 a−1, the dimensionless subsidence curve is calculated by the Sroka–Schober model, as shown in Figure 2. It shows that the Sroka–Schober model is superior to the Knothe model, especially at the initial stage. However, it is difficult to obtain the actual values of f and ξ.Figure 2
Dimensionless subsidence curves calculated by the Sroka–Schober model and Knothe Model.


### 2.2.4. Revised Hrries Model

In view of the deficiencies of the Gompertz model and Logistic model in predicting foundation subsidence, Liu et al. [9] established a time function by revising the Harris model, which can be derived as(6)Tt=1−11+w0mtn,where m and n are the model parameters with positive values.The revised Harris model can describe the whole process of surface subsidence. However, the parametersm and n solved by fitting cannot reflect the stratigraphic characteristics. That is, it is difficult to determine which geological and mining factors control these parameters, and this reduces the reliability of prediction results. In addition, (6) contains the parameter w0, which indicates that the model values are affected by the maximum subsidence value. This is not a time function in the complete sense.


### 2.2.5. Improved Knothe Model

Regarding the surface subsidence curve as a symmetrical distribution, Chang and Wang [10] improved the Knothe model as follows:(7)Tt=T1t=0.5e−cτ−t−e−cτ0<t≤τ,T2t=0.52−e−ct−τ−e−cττ<t≤T0,where T0 is the total subsidence time and τ is the time required for the subsidence velocity to reach the maximum, τ=0.5T0.In fact, the surface subsidence curve is an asymmetrical distribution in deep soil strata. Taking the mining of the 1308 working face of Guotun Coal Mine in Shandong Province as an example, Z1, Z2, and Z3 are the surface observation points along the strike. The measured subsidence [13] and the theoretical subsidence calculated by (7) are shown in Figure 3. Obviously, there is a large deviation between the theoretical and measured values, especially at the accelerated and decelerated stages. Therefore, the improved Knothe model cannot accurately and reasonably predict deep soil surface subsidence.Figure 3
Comparison between the measured and theoretical subsidence.


## 2.2.1. Knothe Model

Assuming that the subsidence velocity at a certain moment is directly proportional to the difference between the dynamic subsidence value at this moment and the maximum subsidence value, Knothe [6] obtained a differential function as follows:(1)dwtdt=cw0−wt,where wt is the instantaneous subsidence value of a point on the surface at time t, w0 is the maximum subsidence value of the point, and c is a time coefficient related to lithology.With the initial conditionw(0) = 0, (1) can be derived as(2)wt=w01−e−ct.The Knothe time function can be expressed as(3)Tt=1−e−ct,where Tt is the dimensionless instantaneous subsidence of a point on the surface at time t, Tt=wt/w0.


## 2.2.2. Kowalski Model

Considering the phenomenon of surface subsidence lagging behind underground mining, Kowalski [7] rewrote the Knothe model as(4)Tt=qt1−Ae−ct−τop,qt=0,t<τop,1,t≥τop,where A is a model parameter and τop is the time from initial mining to the appearance of surface subsidence.Obviously, the Kowalski model still cannot directly reflect the variation law of subsidence velocity.


## 2.2.3. Sroka–Schober Model

Considering the convergence and compaction of porous rock strata, Sroka et al. [8] proposed a time function with two parameters, which can be defined by(5)Tt=1+ξf−ξe−ft−ff−ξe−ξt,where f is the relative convergence rate of rock strata, f = 0.01 a−1 means that the annual convergence is 1% of the original volume, ξ is the time coefficient of the overlying rock strata, when the roof is managed by caving method, and the value of ξ ranges from 20 to 70 a−1.The Sroka–Schober model is used for the prediction of rock mass subsidence, while it can also calculate surface subsidence by adjusting the model parameter values. For example, whenξ = 20 a−1, f = 2 a−1, the dimensionless subsidence curve is calculated by the Sroka–Schober model, as shown in Figure 2. It shows that the Sroka–Schober model is superior to the Knothe model, especially at the initial stage. However, it is difficult to obtain the actual values of f and ξ.Figure 2
Dimensionless subsidence curves calculated by the Sroka–Schober model and Knothe Model.


## 2.2.4. Revised Hrries Model

In view of the deficiencies of the Gompertz model and Logistic model in predicting foundation subsidence, Liu et al. [9] established a time function by revising the Harris model, which can be derived as(6)Tt=1−11+w0mtn,where m and n are the model parameters with positive values.The revised Harris model can describe the whole process of surface subsidence. However, the parametersm and n solved by fitting cannot reflect the stratigraphic characteristics. That is, it is difficult to determine which geological and mining factors control these parameters, and this reduces the reliability of prediction results. In addition, (6) contains the parameter w0, which indicates that the model values are affected by the maximum subsidence value. This is not a time function in the complete sense.


## 2.2.5. Improved Knothe Model

Regarding the surface subsidence curve as a symmetrical distribution, Chang and Wang [10] improved the Knothe model as follows:(7)Tt=T1t=0.5e−cτ−t−e−cτ0<t≤τ,T2t=0.52−e−ct−τ−e−cττ<t≤T0,where T0 is the total subsidence time and τ is the time required for the subsidence velocity to reach the maximum, τ=0.5T0.In fact, the surface subsidence curve is an asymmetrical distribution in deep soil strata. Taking the mining of the 1308 working face of Guotun Coal Mine in Shandong Province as an example, Z1, Z2, and Z3 are the surface observation points along the strike. The measured subsidence [13] and the theoretical subsidence calculated by (7) are shown in Figure 3. Obviously, there is a large deviation between the theoretical and measured values, especially at the accelerated and decelerated stages. Therefore, the improved Knothe model cannot accurately and reasonably predict deep soil surface subsidence.Figure 3
Comparison between the measured and theoretical subsidence.


## 3. New Time Function Model

### 3.1. Model Solution

Based on the expression of the Knothe model, assume thatφt is an arbitrary function, and φt>0. With the hysteretic characteristics of surface subsidence, the Knothe model can be improved to(8)Tt=qt1−φt−ct−τop,where qt, τop, and c are the same as those of the Kowalski model.Obviously, (8) satisfies the boundary conditions of surface subsidence, namely, t≤τop, Tt=0; t=∞, Tt=1.According to reference [14], the parameter c can be solved by(9)c=−vln0.022H/tanβ+s,where H and v are the average mining depth and advancing rate of the working face, respectively, s is the displacement of the inflection point, and β is the main influence angle.For the 1308 working face, there areH = 780 m, v = 1200 m/a, s = 0, and tan β = 1.6. It can be solved from (9) that c = 4.8 a−1. Based on the measured results of points Z1, Z2, and Z3, the φt can be inverted by (8). With the mean fitting method, φt varies approximately linearly with time, as shown in Figure 4.Figure 4
Fittingφt at different time periods.Based on the abovementioned analysis,φt can be expressed as(10)φt=1+kt−τop,where k is a model coefficient. Based on the dimensional analysis, the dimension of k is the same as the lithology time coefficient c.Substituting (10) into (8),(11)Tt=qt1−1+kt−τop−ct−τop.Equation (11) is the new time function for predicting the deep soil surface subsidence, which can be called the composite function model according to its form.When the opening has reached its critical size, the maximum subsidence is equal to 0.98w0. Assume that tc is the time when the face has reached full mining; then, k can be solved by (11) as follows:(12)k=0.02−ctc−τop−1tc−τop.


### 3.2. Characteristic Analysis

For convenience of expression, assumetm=t−τop≥0, and then, q(t) = 1. Thus, (11) can be rewritten as(13)Ttm=1−1+ktm−ctm.Solving the first derivative of (13), the dimensionless subsidence velocity vtm can be derived as(14)vtm=1+ktm−ctmcln1+ktm+cktm1+ktm.For the 1308 working face,tc = 1.8 a and τop = 0.07 a. It can be calculated from (12) that k = 0.93 a−1. Substituting parameters into (13) and (14), the curves of Ttm and vtm can be obtained as shown in Figures 5 and 6, respectively.Figure 5
TheoreticalTtm of the composite function model.Figure 6
Dimensionlessvtm curve of the composite function model.Figure5 shows that the theoretical results are consistent with the measured data, indicating that the composite function model can accurately and reasonably predict deep soil surface subsidence. Figure 6 shows that the new model can well describe the whole subsidence process. Therefore, the composite function model is superior to the existing time function models in predicting deep soil surface subsidence.


## 3.1. Model Solution

Based on the expression of the Knothe model, assume thatφt is an arbitrary function, and φt>0. With the hysteretic characteristics of surface subsidence, the Knothe model can be improved to(8)Tt=qt1−φt−ct−τop,where qt, τop, and c are the same as those of the Kowalski model.Obviously, (8) satisfies the boundary conditions of surface subsidence, namely, t≤τop, Tt=0; t=∞, Tt=1.According to reference [14], the parameter c can be solved by(9)c=−vln0.022H/tanβ+s,where H and v are the average mining depth and advancing rate of the working face, respectively, s is the displacement of the inflection point, and β is the main influence angle.For the 1308 working face, there areH = 780 m, v = 1200 m/a, s = 0, and tan β = 1.6. It can be solved from (9) that c = 4.8 a−1. Based on the measured results of points Z1, Z2, and Z3, the φt can be inverted by (8). With the mean fitting method, φt varies approximately linearly with time, as shown in Figure 4.Figure 4
Fittingφt at different time periods.Based on the abovementioned analysis,φt can be expressed as(10)φt=1+kt−τop,where k is a model coefficient. Based on the dimensional analysis, the dimension of k is the same as the lithology time coefficient c.Substituting (10) into (8),(11)Tt=qt1−1+kt−τop−ct−τop.Equation (11) is the new time function for predicting the deep soil surface subsidence, which can be called the composite function model according to its form.When the opening has reached its critical size, the maximum subsidence is equal to 0.98w0. Assume that tc is the time when the face has reached full mining; then, k can be solved by (11) as follows:(12)k=0.02−ctc−τop−1tc−τop.


## 3.2. Characteristic Analysis

For convenience of expression, assumetm=t−τop≥0, and then, q(t) = 1. Thus, (11) can be rewritten as(13)Ttm=1−1+ktm−ctm.Solving the first derivative of (13), the dimensionless subsidence velocity vtm can be derived as(14)vtm=1+ktm−ctmcln1+ktm+cktm1+ktm.For the 1308 working face,tc = 1.8 a and τop = 0.07 a. It can be calculated from (12) that k = 0.93 a−1. Substituting parameters into (13) and (14), the curves of Ttm and vtm can be obtained as shown in Figures 5 and 6, respectively.Figure 5
TheoreticalTtm of the composite function model.Figure 6
Dimensionlessvtm curve of the composite function model.Figure5 shows that the theoretical results are consistent with the measured data, indicating that the composite function model can accurately and reasonably predict deep soil surface subsidence. Figure 6 shows that the new model can well describe the whole subsidence process. Therefore, the composite function model is superior to the existing time function models in predicting deep soil surface subsidence.


## 4. Engineering Test Analysis

The Longgu Coal Mine is located in Juye Coalfield in Shandong Province, with an average soil thickness of 650 m. As the first mining of working face, 1301 was mined from September 2009 to April 2011. The maximum surface subsidence occurred at the observation point of K131, with a subsidence value of 3484 mm. Therefore, the point K131 was selected as the test object. According to the observation report [15], the main elements of the 1301 working face are as follows: H = 807 m, v = 1535 m/a, tc = 0.62 a, τop = 0.08 a, s = 50 m, and tanβ = 2.2. It can be solved from (9) and (12) that c = 7.2 a−1 and k = 3.2 a−1. Thus, the theoretical subsidence of K131 calculated from the Knothe model and composite function model can be obtained by (3) and (13), respectively. The theoretical and measured subsidence of K131 are shown in Figure 7.Figure 7
Dimensionless subsidence of K131.As seen from Figure7, the composite function model can well describe the whole subsidence process, while the Knothe function model cannot. Assume that the dimensionless measured subsidence is TM, and then, the mean square error between the measured and theoretical subsidence can be calculated as follows:(15)σ=∑i=1nTMi−Titm2n=∑i=1nΔTi2n,where n is the number of observations and ΔTi is the difference between the measured and theoretical subsidence of nondimension.For the calculation result of the Knothe function model, there isσ=1.461/28=0.228; for the composite function model, there is σ=0.046/28=0.041. The results prove that the composite function model is more accurate to predict deep soil surface subsidence under different mining conditions.


## 5. Discussion

To discuss the application of the composite function model in thin soil strata, the observation point of A27 above the middle of the 29401 working face is taken as the research object, as seen in reference [16]. The main parameters of the 29401 working face are as follows: H = 260 m, v = 730 m/a, s = 28.6 m, and tanβ = 2.1. Then, c = 9.24 a−1 is solved by (9). According to the measured subsidence of A27, the parameter values are approximately selected as follows: tc = 0.8 a and τop = 0.3 a. It can be calculated from (12) that k = 2.66 a−1. Thus, the theoretical subsidence of A27 calculated from the Knothe model and composite function model can be obtained by (3) and (13), respectively. When c = 9.24 a−1 and 16 a−1 and τ=tc−τop/2=0.25a, the theoretical subsidence of A27 calculated from the improved Knothe model can be obtained by (7). The theoretical and measured subsidence of A27 are shown in Figure 8.Figure 8
Dimensionless subsidence of A27.The results show that, under the condition of thin soil strata, the composite function model is superior to the Knothe model, and the improved Knothe model is optimal whenc = 16 a−1. However, when c = 9.24 a−1, the improved Knothe model is significantly smaller in predicting the maximum subsidence, which indicates that the improved Knothe model is constrained in exact parameter values.


## 6. Conclusions

Based on an analysis of existing time function models, the composite function model was proposed. All parameters in the new model can be obtained by specific formulas. The conclusions are as follows:(1)
The composite function model can well describe the whole subsidence process of deep soil surface. Compared with the measured subsidence, the prediction results have high precision.(2)
Model parameters solved by specific formulas are reasonable, which enhances the reliability of the prediction results under different mining conditions.(3)
The composite function model is significantly superior to the Knothe model and easier to solve and apply than the improved Knothe model.(4)
As the innovation and expansion of the mining subsidence theory, the new model has an important guiding significance for preventing subsidence geological disasters and determining the coal mining time under the buildings, the railways, and the water bodies.



---

*Source: 4279401-2020-07-09.xml*