Application of the Improved Knothe Time Function Model in the Prediction of Ground Mining Subsidence: A Case Study from Heze City, Shandong Province, China

: Taking into account the inadequacy of the Knothe time function model to predict the dynamic surface subsidence caused by underground mining, a new hypothesis is proposed, and the improved Knothe time function model is established. Theoretical analysis shows the improved model agrees well with surface subsidence dynamic change, velocity change, and acceleration change rules. Combined with ﬁeld measured data, the probability integral method, dual-medium method, and least square method are adopted to determine the time inﬂuence parameter C and the model order n . Based on monitoring data from four monitoring stations in the Guotun coal mine subsidence basin strike main proﬁle from Heze city, Shandong Province, China, the accuracies of the Knothe time function and improved model are compared and analyzed. Results show the improved model can accurately describe the dynamic surface subsidence process and subsidence velocity with mining time. The average relative standard error between the predicted and measured values is only 4.8%—far lower than the Knothe time function model is 23%, verifying the improved model’s accuracy and reliability. show that the surface subsidence initial value is not a ﬀ ected by the model order n , while the time when the surface subsidence reaches the ﬁnal value is a ﬀ ected by the model order n . The greater n is, the less time it takes the surface to reach the ﬁnal subsidence value. As the value of n increases, the maximum subsidence velocity increases, the time of peak and valley acceleration will be shortened, with the absolute peak value and valley acceleration also increasing.


Introduction
The surface subsidence and deformation caused by underground mining is a complex four-dimensional problem varying with a continuous time-space function [1,2]. With progressive mining of the working face, the surfaces' subsidence undergoes three stages of subsidence, namely, slow deformation stage, rapid deformation stage, and basic stable stage. The slow deformation stage only accounts for 10%-15% of the final subsidence and causes minor ground deformation. The rapid deformation stage is the principal subsidence corresponding to about 75% of the final subsidence and causes major ground deformation. The basic stable stage may vary from 5% to 10% of the final subsidence [3]. The deformation that occurs in each of these stages leads to various forms of damage, such as that to buildings, water facilities, and roads. Geological disasters resulting from ground subsidence cause huge economic losses in China every year. Researching the ground subsidence characteristics caused by underground mining and accurately predicting the deformation caused by ground subsidence is one of the most important aspects of protecting buildings and structures above the gob [4][5][6][7].
Surface subsidence is a widespread problem frequently caused by underground mining. In the past several decades, a number of studies have focused on using the geometric theory prediction method, profile function prediction method, mechanical model prediction method, interferometric synthetic aperture radar (InSAR) method, GPS (Global Positioning System) method, and numerical simulation prediction method to estimate the ground subsidence and deformation caused by underground coal mining. Kwinta et al. [8] proposed a one-parameter Knothe's function model based on a growth law and a two-parameter Sroka-Schober's function based on double use of the growth law considering the convergence of extracted excavations and delaying activity of overlain rock. Sheorey et al. [9] proposed a new functional influence method and additional modeling subsidence asymmetry functions for the extraction of edge effects, with the method's accuracy examined by four cases. Lian et al. [10] extended an existing classical dynamic subsidence prediction model, proposing new potential research avenues offered by cellular automata (CA) models. The model verified the deformation measurements collected above a coalmine panel in the Huainan mining area in Anhui Province, China. Tugrul et al. [11] proposed a new subsidence prediction approach using FEM (Finite Element Models) and satellite measurements and gave two case studies as examples for the application of the proposed approach. Hu et al. [12] established a model for calculating the Knothe time function parameters based on a full subsidence angle, verifying the proposed model's reliability and applicability by comparing the predicted and measured values. Chen et al. [13] proposed using intensity tracking technique and the small baseline subset (SBAS) interferometry technique to obtain complete monitoring results. The comparison between experiment monitoring results and GPS measurements showed that the proposed method has the capability to detect and measure the large gradient deformation. Yang et al. [14] proposed a method for completely retrieving three-dimensional (3-D) mining-induced displacements with OT (Offset Tracking)-derived observations of LOS (Line-of-Sight) deformation from a single amplitude pair of SAR (Single Amplitude Pair), and this method had been successfully used in the Daliuta coal mining area. Janez Rošer et al. [15] analyzed surface subsidence above the Velenje Coal Mine's underground mining sites using a modified sigmoid function (logistic function surrogate) and considered the time estimation of the next and final epoch measurement. Luo et al. [16] proposed a range split spectrum interferometry assisted phase unwrapping (R-SSIaPU) method to monitor large gradient mining-induced surface movements, this method agreed with GPS with differences of approximately 4.2 cm in the Tangjiahui mine.
Due to numerous parameters, an extremely complex velocity and acceleration prediction function, and an unreasonable parameter calculation methodology, most of the prediction models are not only different from the actual values but also not conducive to field application. Hence, a dynamic surface subsidence prediction model with high accuracy and few parameters is necessary. In this paper, a new model, based on the Knothe time function model, to predict the dynamic surface subsidence is proposed, and a simple, feasible parameter determination method is provided. Finally, the model accuracy is verified by comparing to the 1301 working face monitoring data in the Guotun coal mine in Shandong Province, China.

The Knothe Time Function Model
In 1952, Knothe used mathematical methods to describe the complex time-space process of surface subsidence caused by underground mining and made the following assumptions: The surface subsidence rate dW(t) dt is proportional to the difference between the final surface subsidence W 0 and the dynamic subsidence W(t) at time t, which may be expressed as [17][18][19]: where c is a time influence parameter related to the overburdened rock mass physical and mechanical properties. This parameter is expressed in the unit 1/year. By separating variables and combining with boundary conditions t = 0, W(t) = 0, the relationship between subsidence and time is obtained as follows: The Knothe time function model is widely used to predict the actual surface subsidence process due to its simple form and few undetermined parameters. However, with expanding research and practice, this model does not accurately describe the dynamic processes of the subsidence [19,20].
The first and second derivatives of Equation (2) are used to obtain the surface subsidence velocity and acceleration as follows: Using Equations (2) and (3), the surface subsidence rule described by the Knothe time function model is shown in Figure 1a. The figure shows that when t = 0, the surface subsidence velocity is at the maximum value cW 0 . With the passage of time, the velocity gradually decreases and finally approaches 0. That is to say, the surface subsidence gradually increases from 0, but the increasing range is shortened. The surface subsidence acceleration curve shows that from the beginning of surface subsidence to final stability, the whole subsidence process acceleration is always less than 0. In other words, surface subsidence velocity is always decreasing.
Using Equations (2) and (3), the surface subsidence rule described by the Knothe time function model is shown in Figure 1a. The figure shows that when = 0, the surface subsidence velocity is at the maximum value . With the passage of time, the velocity gradually decreases and finally approaches 0. That is to say, the surface subsidence gradually increases from 0, but the increasing range is shortened. The surface subsidence acceleration curve shows that from the beginning of surface subsidence to final stability, the whole subsidence process acceleration is always less than 0. In other words, surface subsidence velocity is always decreasing.
Numerous research and field monitoring data show that the surface subsidence curve caused by coal mining is approximately an "S"-type curve [21,22]. When the distance between set-up entry and the working face reaches the surface starting distance, the surface starts sinking. The subsidence velocity is 0 at that point. Then, the subsidence slowly increases, but subsidence velocity rapidly increases. As the working face continuously advances, the surface subsidence increases significantly, and the subsidence velocity gradually increases to the maximum value. When the working face has advanced a sufficient distance away from the entry set-up, the surface subsidence shape curve remains basically unchanged, and the subsidence continues to increase, but the velocity increase is very slow and eventually tends toward 0. Figure 1b shows the actual measured subsidence, subsidence velocity, and subsidence acceleration curve at the z7 monitoring point on the 35,101 working face in the Sandaogou mine [23] in Shanxi Province, China. Comparing Figure 1a,b, the Knothe time function model is clearly quite different from the actual situation in describing surface subsidence, subsidence velocity, and subsidence acceleration. The Knothe time function model is suitable only for describing the basic stable stage and cannot accurately predict the whole dynamic surface subsidence process [24].   Numerous research and field monitoring data show that the surface subsidence curve caused by coal mining is approximately an "S"-type curve [21,22]. When the distance between set-up entry and the working face reaches the surface starting distance, the surface starts sinking. The subsidence velocity is 0 at that point. Then, the subsidence slowly increases, but subsidence velocity rapidly increases. As the working face continuously advances, the surface subsidence increases significantly, and the subsidence velocity gradually increases to the maximum value. When the working face has advanced a sufficient distance away from the entry set-up, the surface subsidence shape curve remains basically unchanged, and the subsidence continues to increase, but the velocity increase is very slow and eventually tends toward 0. Figure 1b shows the actual measured subsidence, subsidence velocity, and subsidence acceleration curve at the z7 monitoring point on the 35,101 working face in the Sandaogou mine [23] in Shanxi Province, China. Comparing Figure 1a,b, the Knothe time function model is clearly quite different from the actual situation in describing surface subsidence, subsidence velocity, and subsidence acceleration. The Knothe time function model is suitable only for describing the basic stable stage and cannot accurately predict the whole dynamic surface subsidence process [24].

The Establishment of an Improved Knothe Time Function Model
By analyzing the Knothe time function hypothesis and Equation (3), Equation (4) can be obtained: Equation (4) shows that the surface subsidence velocity to acceleration ratio and the (n-1) acceleration derivative time to n-th derivative ratio described by the Knothe time function model are time-independent, and are all constant values − 1 c at any time. This is obviously inconsistent with the actual dynamic surface subsidence velocity and acceleration relationship, leading to a significant difference between the predicted surface subsidence and the field monitoring values. The Knothe time function model premise assumption is special, not general, so an improved Knothe time function model is required and a new assumption is made: the n-th surface subsidence derivative d n W(t) d n t is proportional to the n-th power of the difference between the final surface subsidence W 0 and the dynamic subsidence W(t) at time t, a general differential Equation (5) can be obtained: where C is a time influence parameter of the improved Knothe time function model, with this parameter expressed in the unit 1/(year n ). The solution of Equation (5) is as follows: where A is the integral constant. Considering the initial condition t = 0, W (t) = 0, the integral constant A = W 0 , then Equation (6) can be simplified as: Equation (7) shows that the dynamic surface subsidence function W(t) has a power-law relationship with time. When 0 ≤ t < ∞, W(t) ∈ (0, W 0 ), the surface dynamic subsidence velocity can be derived from Equation (7): The dynamic subsidence acceleration obtained from the Equation (8) derivative is as follows: When n = 1, upon comparing Equations (2) and (3), and Equations (7)-(9), we can find that the Knothe time function model is actually a special case of the improved Knothe time function model. If a(t) = 0, Equation (9) can be converted into: According to Equations (8)- (10), in the interval 0, n−1 C 1 n , the dynamic subsidence acceleration a(t) is always greater than 0, the dynamic surface subsidence velocity gradually increases, and the dynamic surface subsidence gradually increases from 0. In the interval n−1 C 1 n , ∞ , the surface dynamic subsidence acceleration a(t) is always less than 0, the surface dynamic subsidence velocity V(t) gradually decreases and finally approaches 0, and the dynamic surface subsidence increase gradually slows and finally approaches W 0 . According to the function extreme value principle, when t = n−1 C 1 n , the maximum dynamic subsidence velocity value V(t) max is shown in Equation (11): The corresponding surface subsidence is shown in Equation (12): The surface subsidence curves of different models order n, represented by the improved Knothe time function model, the subsidence velocity curve, and the subsidence acceleration curve are shown in Figures 2-4, respectively. Figures 2-4 demonstrate that the improved Knothe time function model is consistent with the actual surface subsidence, and can be used to describe the dynamic surface subsidence process. Figures 2-4 show that the surface subsidence initial value is not affected by the model order n, while the time when the surface subsidence reaches the final value is affected by the model order n. The greater n is, the less time it takes the surface to reach the final subsidence value. As the value of n increases, the maximum subsidence velocity increases, the time of peak and valley acceleration will be shortened, with the absolute peak value and valley acceleration also increasing.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 14 According to Equations (8)- (10), in the interval 0， , the dynamic subsidence acceleration ( ) a t is always greater than 0, the dynamic surface subsidence velocity gradually increases, and the dynamic surface subsidence gradually increases from 0. In the interval ，∞ , the surface dynamic subsidence acceleration ( ) is always less than 0, the surface dynamic subsidence velocity ( ) gradually decreases and finally approaches 0, and the dynamic surface subsidence increase gradually slows and finally approaches . According to the function extreme value principle, when = , the maximum dynamic subsidence velocity value ( ) is shown in Equation (11): The corresponding surface subsidence is shown in Equation (12): The surface subsidence curves of different models order n, represented by the       The greater n is, the less time it takes the surface to reach the final subsidence value. As the value of n increases, the maximum subsidence velocity increases, the time of peak and valley acceleration will be shortened, with the absolute peak value and valley acceleration also increasing.

Determination of Model Parameters
When we use the improved Knothe time function to predict progressive subsidence, the key process is to determine the model parameters, which are closely related to the mechanical characteristics of the overburden strata, the coal mining speed, and the roof management method [25][26][27].

Determination of Model Order n
Chang and Wang [28] showed that when the surface point subsidence velocity reaches the maximum value, the surface subsidence is about half the final subsidence value. According to Equation (12), we obtain the following formula: The model order = 3.26.

Determination of Model Parameters
When we use the improved Knothe time function to predict progressive subsidence, the key process is to determine the model parameters, which are closely related to the mechanical characteristics of the overburden strata, the coal mining speed, and the roof management method [25][26][27].

Determination of Model Order n
Chang and Wang [28] showed that when the surface point subsidence velocity reaches the maximum value, the surface subsidence is about half the final subsidence value. According to Equation (12), we obtain the following formula: The model order n = 3.26.

Determination of Time Influence Parameter C
According to Equation (7), the surface subsidence curve with different time influence parameter C is shown in Figure 5.

Determination of Time Influence Parameter C
According to Equation (7), the surface subsidence curve with different time influence parameter C is shown in Figure 5.  Comparing Figures 2 and 5, it can be seen that the impact of time influence parameter C on the surface subsidence curve is similar to the model order n, both of which have no effect on the initial and final surface subsidence values, but only affect the time to reach the final subsidence value. The larger C is, the less time it takes for the surface subsidence to reach a stable state. When the model n is determined, the key to the improved time function model accuracy is the determination of C.
According to the probability integral method, when the mining range reaches the critical value of full mining, the maximum subsidence is approximately equal to 0.98W 0 [29,30]. If we assume that the mining velocity is v and the critical gob dimension is L f , the critical time to achieve full mining can be expressed as L f v . According to Equation (7), the time influence parameter C calculation formula is as follows: Because the mining velocity is easily determined, the key to finding time influence parameter C is to determine the critical gob dimension when full mining is realized. Presently, three methods are used to calculate the critical gob dimension [31]: (1) Bedrock decision method: Determines final surface subsidence based only on the bedrock with no loose layer consideration. Obviously, the L f determined by this method is only suitable for a shallow loose layer, not applicable for mining areas with thick loose layers, such as those found in central and eastern China. This method's calculation model is shown in Figure 6a and L f = 2H j tanϕ j relating to the geometric relationship. (2) Comprehensive impact parameter method: Regards the loose layer and bedrock as a composite medium, so the bedrock and loose layers mining influence angle are equal. This method's calculation model is seen in Figure 6b, with L f = 2H 0 tanϕ relating to the geometric relationship. This method ignores the mechanical properties difference between the loose layer and the bedrock, so errors exist. Generally speaking, the loose layer mining influence angle is larger than that of the bedrock.
(3) Dual-medium method: The coal mining influence on the bedrock and loose layer deformation is considered comprehensively. The critical gob dimension obtained is close to the actual value.
The calculation model is shown in Figure 6c. According to the geometric relationship, L f is determined as follows: where H 0 is the average seam depth, ϕ is the full subsidence angle, H j is the bedrock thickness, ϕ j is the bedrock full subsidence angle, H s is the loose layer thickness, and ϕ s is the loose layer full subsidence angle.
3) Dual-medium method: The coal mining influence on the bedrock and loose layer deformation is considered comprehensively. The critical gob dimension obtained is close to the actual value. The calculation model is shown in Figure 6c. According to the geometric relationship, is determined as follows: where is the average seam depth, is the full subsidence angle, is the bedrock thickness, is the bedrock full subsidence angle, is the loose layer thickness, and is the loose layer full subsidence angle. Substituting Equation (15) into Equation (14), the calculation formula of C is as follows: Substituting Equation (15) into Equation (14), the calculation formula of C is as follows: A significant amount of monitoring data shows that it is not precise and that the corresponding surface subsidence is just half of W 0 when the subsidence velocity reaches the maximum value. Therefore, there is an error in the value of model order n. The time influence parameter C obtained according to Equation (16) will also produce a certain error. However, in the absence of significant field surface subsidence monitoring data, the above parameter determination method can be used as a reference.
When there is sufficient field monitoring data, the least square method can be used to determine parameters n and C; the multiple group subsidence data is selected from the monitoring points, Equation (7) is used to fit the subsidence data, multiple groups of n and C are obtained, and then the average value is taken as the final parameter. This method is easy to operate and has high accuracy in obtaining parameters.

Model Validation
To verify the improved Knothe time function model accuracy and reliability established in this paper, based on the surface movement observation station field monitoring data above the 1301 working face in Guotun coal mine from Heze city, Shandong Province, China, the monitoring data and the model predicted results were compared and analyzed. The surface movement observation stations' layout is shown in Figure 7, and the spacing between each observation station is 30 m. Due to a long observation time span and a large amount of monitoring data, clearer monitoring and prediction data comparison and analysis were achieved by randomly selecting only four stations, z38, z41, z44, and z47, as the research study objects for improved and traditional Knothe time function model accuracy analysis in surface subsidence prediction.
The working face 1301 elevation is from −748.9 to −817 m; the coal seam thickness is 1 to 3.6 m, with the average being 2.9 m, and the coal seam dip angle is 12 • ; the roof is managed by full caving methodology. Due to mining area limitations, the observation point along the strike direction sinks to a flat bottom, with no flat bottom phenomenon along the inclined surface, that is to say, the working face strike is full mining, and the working face incline is not. Table 1 shows the monitoring data of z38, z41, z44, and z47 for each period. The traditional and improved Knothe time function models were, respectively, used to predict the surface subsidence and subsidence velocity. The two models' prediction comparison results are shown in Figures 8 and 9. and the model predicted results were compared and analyzed. The surface movement observation stations' layout is shown in Figure 7, and the spacing between each observation station is 30 m. Due to a long observation time span and a large amount of monitoring data, clearer monitoring and prediction data comparison and analysis were achieved by randomly selecting only four stations, z38, z41, z44, and z47, as the research study objects for improved and traditional Knothe time function model accuracy analysis in surface subsidence prediction. The working face 1301 elevation is from −748.9 to −817 m; the coal seam thickness is 1 to 3.6 m, with the average being 2.9 m, and the coal seam dip angle is 12°; the roof is managed by full caving methodology. Due to mining area limitations, the observation point along the strike direction sinks  As can be seen from the monitoring data in Table 1 and Figure 8, the duration of the rapid deformation stage of z37, z41, z44, and z47 monitoring sites was 310 days, 303 days, 276 days, and 198 days, respectively. Moreover, the duration of the rapid deformation stage of z44 and z47 monitoring sites was shorter than that of the slow deformation stage and basic stable stage. However, the data of the four monitoring sites indicated that the dynamic surface subsidence curve was approximately "S"-type. Therefore, the proposed model in the manuscript can be used to predict the vertical surface deformation. In addition, due to the lack of surface horizontal movement monitoring data and the focus of the manuscript is the vertical surface deformation. The horizontal deformation law of the surface will be discussed in the following research. Figures 8 and 9 demonstrate that the improved Knothe time function model prediction accuracy for surface subsidence is significantly better than the Knothe time function model, and can comprehensively reflect the whole dynamic surface subsidence velocity process from beginning acceleration to attenuation. The prediction results were consistent with the monitoring data. The prediction value of the Knothe time function model was quite different from the monitoring data. The improved Knothe time function model parameters obtained by the least square method are shown in Table 2.
Standard deviation and relative standard deviation can measure the deviation between monitored and predicted values [20,31]. To further verify the improved Knothe time function model accuracy, the standard deviation and relative standard deviation of the monitoring value and the predicted value of four monitoring points were analyzed without considering the monitoring error. The calculation formula is shown in Equation (17), and the calculation results are shown in Table 3. The improved Knothe time function model's maximum relative standard deviation was only 5.7%, and the average relative standard deviation was only 4.8%. In contrast, the Knothe time function model maximum relative standard deviation was as high as 24.8%, and the average relative standard deviation was 23%, much higher than the improved Knothe time function model, which confirms the improved Knothe time function model accuracy and reliability established in this paper.
where m is standard deviation; f is relative standard deviation; W m is monitored value; W p is predicted value; k is monitoring data volume.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 14 to a flat bottom, with no flat bottom phenomenon along the inclined surface, that is to say, the working face strike is full mining, and the working face incline is not. Table 1 shows the monitoring data of z38, z41, z44, and z47 for each period. The traditional and improved Knothe time function models were, respectively, used to predict the surface subsidence and subsidence velocity. The two models' prediction comparison results are shown in Figures 8 and  9.  As can be seen from the monitoring data in Table 1 and Figure 8, the duration of the rapid deformation stage of z37, z41, z44, and z47 monitoring sites was 310 days, 303 days, 276 days, and 198 days, respectively. Moreover, the duration of the rapid deformation stage of z44 and z47 the data of the four monitoring sites indicated that the dynamic surface subsidence curve was approximately "S"-type. Therefore, the proposed model in the manuscript can be used to predict the vertical surface deformation. In addition, due to the lack of surface horizontal movement monitoring data and the focus of the manuscript is the vertical surface deformation. The horizontal deformation law of the surface will be discussed in the following research.  Figures 8 and 9 demonstrate that the improved Knothe time function model prediction accuracy for surface subsidence is significantly better than the Knothe time function model, and can comprehensively reflect the whole dynamic surface subsidence velocity process from beginning acceleration to attenuation. The prediction results were consistent with the monitoring data. The prediction value of the Knothe time function model was quite different from the monitoring data. The improved Knothe time function model parameters obtained by the least square method are shown in Table 2.

Conclusions
(1) Taking into account the inadequacy of the Knothe time function model to predict the dynamic surface subsidence, an improved Knothe time function model was established. The new model can accurately predict the dynamic surface subsidence caused by underground mining.
(2) The time influence parameter C and the model order n do not affect the initial and final surface subsidence values but only affect the time to reach the final subsidence value. The larger C and n are, the less time it takes to reach the final subsidence value. When the surface subsidence monitoring data is lacking, the time influence parameter C and the model order n can be determined according to the probability integral and the dual-medium methods. When a significant amount of monitoring data exists, the time influence parameter C and the model order n can be determined according to the least-squares methods.
(3) Four surface movement observation station monitoring data sets on the mine subsidence basin's main profile in the Guotun coal mine were randomly selected. The monitoring data indicated that the dynamic surface subsidence curve is approximately "S"-type. Therefore, the proposed model in the manuscript can be used to predict the vertical surface deformation. By comparing the curve of improved Knothe time function model, Knothe time function model, and monitoring data, the results showed that the improved Knothe time function model surface subsidence value agrees well with the monitoring value, having better accuracy than the Knothe time function model. The average relative standard deviation between predicted and monitoring values was only 4.8%, much lower than the 23% of the Knothe time function model. The results verify the accuracy and reliability of the improved Knothe time function model.