Research Article | Open Access

# A One-Dimensional Fractional-Derivative Viscoelastic Model for the Aquitard Consolidation of an Aquifer System

**Academic Editor:**Constantinos Loupasakis

#### Abstract

Land subsidence resulting from the withdrawal of groundwater is one of the problems of major importance to geosciences and geomechanics. In order to analyze the phenomena of land subsidence, this study attempts to constitute a more reliable prediction model to simulate the consolidation process of the aquifer system under the drawdown of groundwater. A fractional-derivative viscoelastic model is introduced to characterize the rheological properties of the aquitard, and the equation governing the consolidation process is established on the basis of the one-dimensional consolidation theory. The semianalytical solutions for the pore-water pressure, the settlement of aquitard, and the degree of consolidation are deduced using both finite sine transform and Laplace transform. Then, the comparisons of two special cases are performed to validate the correctness of the proposed model. The variations in the degree of consolidation and settlement of the aquitard are simulated to analyze the consolidation behavior of the aquifer system, and the influences of model parameters and the drawdown speed of groundwater are investigated.

#### 1. Introduction

Land subsidence including both gentle downward warping and the sudden sinking of parts of the ground surface is one kind of widespread geological disaster that occurs in many countries and regions all over the world [1–7]. In China, it has been found that land subsidence mainly occurred in 17 provinces in the eastern and central regions [5, 6]. Commonly, land subsidence can result from natural processes (such as mineral dissolution), human activities, or a combination of both. Perhaps, one of the best-known causes is the withdrawal of groundwater, which can generate land subsidence by changing the effective stress in the soil of the aquifer system [8]. This phenomenon is particularly found in an aquifer system with an aquitard. In order to simulate the process of land subsidence, some aquifer system models and consolidation solutions have been developed after the seminal work on the consolidation problems conducted by Terzaghi [9] and Biot [10]. In the 1940s, Jacob [11] studied the problem of an elastic artesian aquifer system with linear leakage. After this work, Hantush and Jacob [12] studied the problem of groundwater flow in an infinite leaky aquifer system, and then Hantush [13, 14] extended this study to the leaky aquifer system under the conditions of a constant hydraulic head and constant flow. Selvadurai and Kim [4] studied the withdrawal process of groundwater from saturated porous medium with a deformable skeleton. Bear and Corapcioglu [15] presented a mathematical model for regional land subsidence caused by pumping from aquifer systems. Ren and Su [16] analyzed the land subsidence and rebound induced by pumping and flooding combining observed values from two different locations in Shanghai. Luo et al. [17] derived theoretical solutions for analyzing the settlement process resulting from decompression in an aquifer. Wang et al. [6] studied the land subsidence process resulting from hydraulic head variation in a multilayer aquifer system. Zeng and Wang [18] deduced an analytical solution for studying the consolidation behavior of an aquifer system under unsteady groundwater variation. Xie et al. [19] established a one-dimensional consolidation model for an aquitard in a leakage system and an analytical solution for pore water pressure and consolidation degree were given. Li et al. [20] derived an analytical solution to characterize the drawdown variation in an aquifer system taking into account the large strain behavior. Hu et al. [21] studied the consolidation behavior in a layered aquifer system taking into consideration the structured characteristics of the aquitard. In all of these studies, the aquitard soil was simulated using a linear elastic model for developing simple analytical solutions. However, the soft soil generally shows remarkable rheological properties during the consolidation process [22, 23]. Therefore, other viscoelastic models may be more suitable for describing the consolidation characteristics of the aquitard soil.

In order to capture the rheological properties of the soft soil, a number of viscoelastic models have been introduced to extend the consolidation theories since the 1960s, but most of the existing studies have been undertaken to simulate the consolidation process subjected to external loading for a one-dimensional consolidation system or that for a vertical drain consolidation system [22, 24–27]. For the aquifer system, Wu and Miao [28] developed the viscoelastic theory of flow in an aquifer system using the standard Merchant model. Using this model, Liu et al. [5] simulated the consolidation process of the viscoelastic aquitard induced by the withdrawal of groundwater. Tsai [29] proposed a viscoelastic-plastic model for the aquitard consolidation caused by hydraulic head variations in the aquifer system on the basis of the standard Voigt model. Zhang et al. [30, 31] adopted several models including linear elastic, viscoelastic, elastic-plastic, and viscoelastic-plastic to characterize the deformation behavior of the aquifer system under complex changing patterns of the groundwater level. All of the constitutive models adopted in these studies are standard ones that were defined by derivatives of integer order. As a branch of mathematics, fractional calculus deals with the generalization of integrals and derivatives to all real (and even complex) orders [32, 33]. It has been shown that fractional calculus is an effective mathematical instrument for describing the memory phenomena and rheological properties of various materials. Some early studies indicated that the rheological properties of materials were simulated more accurately by a fractional-derivative viscoelastic model than by an integer-order model (i.e., a standard model) [32–37]. In order to better capture the mechanical properties of geomaterials, Yin et al. [38] adopted fractional calculus to develop a fractional-derivative model and its reasonability was demonstrated by the triaxial test results. Hou et al. [39] presented a fractional-derivative model for frozen soil considering the strengthening and weakening effects. Zhang et al. [40] proposed a fractional-derivative Merchant model to study the time-dependent behavior of a simply supported plate on viscoelastic foundation. Using this model, Zhu et al. [41] predicted the deformation characteristics of Nansha clay subjected to one-dimensional compression. Zhang et al. [33] established a four-element fractional-derivative model to study the creep deformation of the surrounding rock around circular tunnels. For the consolidation problem of soft clay, Wang et al. [42] extended the one-dimensional consolidation theory using a two-element fractional-derivative model. Huang and Li [43] investigated the consolidation process of the vertical drain system installed in soft clay using a fractional-derivative Merchant model. For the aquifer system, however, there are few studies on the consolidation process by considering a fractional-derivative model.

In this paper, the idea of the fractional-derivative model is introduced to develop a general model to present the consolidation process of an aquifer system caused by the withdrawal of groundwater. Incorporating the fractional-derivative Merchant model, the equation governing the consolidation process of the aquifer system is established on the basis of the one-dimensional consolidation theory [9]. The semianalytical solutions for the pore-water pressure, the settlement of an aquitard, and the degree of consolidation are obtained using both the finite sine transform and Laplace transform. Then, the comparisons of two special cases are performed to validate the correctness of the proposed model. The variations in the degree of consolidation and settlement of an aquitard are simulated to analyze the consolidation behavior of the aquifer system, and the influences of model parameters are investigated.

#### 2. Mathematical Model

Consider an aquifer system consisting of an unconfined aquifer (sandy stratum), a confined aquifer (sandy stratum), and a viscoelastic aquitard (clay stratum), as presented in Figure 1. In general, the sandy strata are much more permeable than the clay stratum. If the permeability of the neighboring sandy strata exceeds that of the clay stratum over two orders of magnitude, the flow in the clay stratum can be reasonably considered to be vertical [44]. In this study, the flow in the aquitard is assumed to occur only in a vertical direction. In addition, the aquitard is supposed to deform only along the vertical direction. The weight of the soil is ignored, and the total stress in the aquitard is assumed to be unchanged during the whole consolidation process. The thickness of the aquitard is given as , and the vertical coefficient of permeability is .

In adopting all the hypotheses of the one-dimensional consolidation theory by Terzaghi [9], the equation governing the consolidation process of the aquitard is given by where is the vertical strain of the aquitard; is the pore-water pressure in the aquitard, including the excess pore-water pressure and the static pore-water pressure; and are depth and time; and is the water unit weight, respectively.

Applying the continuity of the hydraulic head, the initial pore-water pressure in the aquitard can be expressed as where and are the initial hydraulic heads in the aquifers above and below the aquitard, respectively. Similarly, at the top and bottom boundaries of the aquitard, the pore-water pressures can be written as where and are the drawdown of the hydraulic heads with time in the aquifers above and below the aquitard, respectively.

Based on the principle of effective stress [9], the vertical effective stress in the aquitard is where represents the total vertical stress in the aquitard, and it is assumed to be constant during the whole consolidation process. At the initial time, the effective vertical stress in the aquitard equals to zero. Thus, we have

The fractional-derivative Merchant model, as shown in Figure 1, is adopted to simulate the rheological characteristics of the aquitard. This model is composed of a fractional-derivative Voigt body in a series with an independent Hookean spring of modulus . The fractional-derivative Voigt body consists of a fractional-derivative element in parallel connection with a Hookean spring of modulus . The fractional-derivative element is defined using the idea of a spring-pot model [34–36]. It has a property wherein its constitutive equation has continuity from the ideal solid state to the ideal fluid state. Based on this model, the vertical strain of the aquitard can be expressed in a convolution form as [40, 43] where stands for the modulus ratio of the Hookean springs; is defined as the creep time of the fractional-derivative Voigt body with the viscosity coefficient ; is the fractional order of the fractional-derivative element, and it is commonly considered to be a dimensional coefficient concerning the memory of the real materials; is called the convolution of two functions; represents the Laplace transform of a real function and its inverse operator is denoted as ; and is a complex variable, respectively. Moreover, it has been found [43] that where is the Dirac delta function and is the time. It is obvious that the substitution of equation (7) into equation (6) gives the constitutive relationship of the linear elastic model if , while it produces that of the standard Merchant model if . That is, both the linear elastic model and standard Merchant model are two special cases of the fractional-derivative Merchant model.

#### 3. Solution Derivation

Substituting equation (6) into equation (1) with consideration of equations (4) and (5), the equation governing the consolidation process of the aquitard can be obtained as where is the coefficient of vertical consolidation and is the two-parameter function of the Mittag-Leffler type. In order to derive the solution of equation (8), both the finite sine transform and Laplace transform (see Appendix) [45, 46] are performed. The transformed variable is defined as where stands for the finite sine transform with respect to .

Applying the transform of to both sides of equation (8) with the definite conditions of equations (2) and (3), we can obtain the following linear algebraic equation: where , , and .

Solving for in equation (10) results in where

Implementing the inverse transform of on both sides of equation (11) produces where where and can be any real number greater than the real part of all the singularities of . For the general case of , equation (14) cannot be worked out analytically; thus, the widely used Crump’s method (see Appendix) [42, 47] is adopted to approximate the result of . In the time domain, the result of can be computed using the following series: where is a real parameter satisfying the condition of .

Substituting equations (5) and (13) into equation (4) produces

The average degree of consolidation is usually used as one of the criteria for assessing the effectiveness of soil consolidation and improvement [48]. In general, there are two common methods widely used to define the average degree of consolidation [49, 50]. One is based on the effective stress (i.e., pore-water pressure) in the soil stratum, and the other is based on the settlement attained by the soil stratum. For a simple case with the assumptions of the linear elastic model and constant permeability, the degrees of consolidation evaluated by these two methods are equal to each other; while for the case with the presence of viscosity, the equality does not hold. Herein, the average degree of consolidation is evaluated for the aquitard by comparing the average effective stress and the total stress as follows: where is the average degree of consolidation of the aquitard.

In general, the bottom boundary of the aquitard connected to the nearly rigid sandy stratum yields zero deformation [29]; thus, the settlement of the aquitard can be obtained by integrating equation (6) along the vertical direction. The result is where is the settlement of the aquitard, and it can be numerically calculated from where

#### 4. Simplification and Discussion

##### 4.1. Simplification for Two Special Cases

In this part, two special cases of the proposed aquifer system are discussed to validate the correctness of the above-derived semianalytical solution. As stated, the fractional-derivative Merchant model can be simplified into the linear elastic model and the standard Merchant model for the cases of and , respectively. When it comes to the aquifer system with the linear elastic aquitard (i.e., ), equation (14) can be obtained as

By substituting equation (21) into equation (13), the pore-water pressure can be expressed as

Equation (22) is the same as the solution that was put forward by Zeng and Wang [18]. Further, for the case of , , , and , this equation can be reduced to where and . Equation (23) is the same as the analytical solution proposed by Xie et al. [19].

When it comes to the aquifer system with the viscoelastic aquitard that obeys the standard Merchant model (i.e., ), equation (14) can be obtained as where

Substituting equation (24) into equation (13) gives the expression of pore-water pressure as

For the case of , , , and , equation (26) can be further simplified into

Equation (27) is the same as the analytical solution proposed by Liu et al. [5].

For verifying the accuracy of the above-derived semianalytical solution, two case studies are analyzed by using the semianalytical solution and the available analytical solutions. One is for the aquifer system with the linear elastic aquitard, and the other is for the aquifer system with the viscoelastic aquitard that obeys the standard Merchant model. For these two case studies, the results of the pore-water pressure are calculated using the semianalytical solution and the analytical solutions of Xie et al. [19] and Liu et al. [5]. The calculated results are shown in Figure 2. The parameters of groundwater level and aquitard soil used in Figures 2(a) and 2(b) are the same as References [5, 19], respectively.

**(a)**

**(b)**

As can be seen from Figure 2, the isochrones of the pore-water pressure given by the semianalytical solution are identical to the analytical results presented in the study of Xie et al. [19] for the aquifer system with the linear elastic aquitard and match well with those presented in the study of Liu et al. [5] for the aquifer system with the viscoelastic aquitard that obeys the standard Merchant model. These comparisons show that the newly proposed semianalytical solution concerning the fractional-derivative Merchant model is consistent with the existing consolidation theories for the aquifer system and is workable for the aquifer system with the linear elastic aquitard and the aquifer system with the viscoelastic aquitard obeying the standard Merchant model. That is, the aquifer system models respectively put forward by Liu et al. [5] and Xie et al. [19] are two simplified cases of the fractional-derivative model presented herein. Moreover, the fractional-derivative model is more general to simulate the consolidation characteristics of the aquifer system during the withdrawal of groundwater.

##### 4.2. Parametric Study and Discussion

A typical example is presented to illustrate the consolidation characteristics of the proposed aquifer system. Actually, the hydraulic head (i.e., groundwater level) in the aquifer system does not decline instantaneously after withdrawing groundwater [5]. In this discussion, the drawdown of the hydraulic head in the aquifer below the aquitard is considered to obey an exponential law. This is shown as follows: where is a dimensionless parameter reflecting the drawdown speed of the hydraulic head in the aquifer below the aquitard. If its value approaches infinite, the drawdown of the hydraulic head happens instantaneously. Equation (28) is the same as that adopted by Liu et al. [5]. For the sake of simplicity, the drawdown of the hydraulic head above the aquitard is considered to be unchanged throughout the time of the aquitard consolidation (i.e., ). The parameters of groundwater level and aquitard soil are supposed to be as follows [5]: m, m, m, m, m/s, and MPa. A set of parametric studies is also conducted to investigate the influences of the fractional-derivative model parameters (i.e., , , and ) and the drawdown speed of the groundwater level (i.e., ). Using the derived semianalytical solutions, the settlement of the aquitard and the degree of consolidation have been calculated, and the results are illustrated in Figures 3–6. In these figures, the logarithmic scale is chosen as that of the vertical time factor in the horizontal coordinate.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

As can be seen in Figures 3–6, the settlement of the aquitard and the degree of consolidation increase gradually with the passing of time after the groundwater level declines. In Figure 3, the fractional order (i.e., ) varies from 0 to 1. It can be clearly observed that with different fractional orders, both the settlement of the aquitard and the degree of consolidation have significant differences in the consolidation process. As shown in Figure 3(a), the settlement curves having different fractional orders intersect with each other. A larger fractional order leads to a smaller settlement at the initial time, while it results in a greater settlement at a later stage. This phenomenon is identical to that of the soil deformation behavior of the soft soil subjected to external loading [41–40] and that of the vertical drain consolidation system [43]. As shown in Figure 3(b) with different fractional orders, the calculated results of the degree of consolidation also appear as a set of curves intersecting with each other. Different from the developing feature of the settlement, an increase in the fractional order initially increases the degree of consolidation, while it subsequently decreases the degree of consolidation. In other words, for the aquifer system featured as having a greater fractional order, the dissipation of the pore-water pressure in the aquitard is faster at the initial time, while it becomes slower at a later stage after withdrawing groundwater. This phenomenon occurs because of the mechanical properties of the fractional-derivative element whose constitutive equation has continuity from the ideal solid state to the ideal fluid state when the fractional order changes from zero to one. For the case of a small fractional order, the fractional-derivative element performs more like the ideal solid (i.e., the spring); thus, its mechanical properties will be mainly controlled by the stiffness and the deformation. For the case of a large fractional order, the fractional-derivative element works more like the ideal fluid (i.e., the dashpot); thus, its mechanical properties will be mainly dominated by the viscosity and the deformation rate coefficient. As shown in Figure 3(a), the deformation attained by the aquitard is very small at the initial stage of the consolidation process; thus, the vertical effective stress produced by the deformation is insignificant and can be ignored. Therefore, the vertical effective stress in the aquitard is mainly controlled by the deformation rate at this stage. As a result, an increase in the fractional order increases the vertical effective stress and accelerates the dissipation of the pore-water pressure further. With the progress of the consolidation, the deformation attained by the aquitard gradually increases and even becomes very large at a later stage, as shown in Figure 3(a); thus, the influences of the deformation upon the vertical effective stress turns out to be significant. At this stage, the deformation mainly controls the variation of the vertical effective stress in the aquitard. As a result, a decrease in the fractional order increases the vertical effective stress and accelerates the dissipation of the pore-water pressure further. This phenomenon is the same as that found in the consolidation of the vertical drain system [43].

Figure 4 illustrates the results of the settlement of the aquitard and the degree of consolidation under different values of the modulus ratio (i.e., ). The value of the modulus ratio changes from 0.1 to 10. In this discussion, both the consolidation time and the settlement of the aquitard are given as dimensionless variables that are associated with the modulus of the independent spring (i.e., ); thus, the effect of is eliminated herein. As a result, the modulus ratio represents the effect of the modulus of the spring (i.e.,) in the fractional-derivative Voigt body. For a given value of , a greater modulus ratio gives a smaller , which means the aquitard is softer. As shown in Figure 4, an increase in the modulus ratio delivers a larger settlement of the aquitard, while it results in a smaller degree of consolidation. In other words, after withdrawing groundwater, the softer aquitard deforms more rapidly while it undergoes a slower increase in the vertical effective stress (i.e., a slower dissipation of the pore-water pressure). At last, the curves in Figure 4(b) tend to be the same regardless of the values of the modulus ratio.

Figure 5 presents the results of the settlement of the aquitard and the degree of consolidation under different values of the creep time factor (i.e., ). The value of the creep time factor varies from 0.01 to 100. Herein, the creep time factor reflects the effect of the viscosity coefficient (i.e., ) upon the consolidation process of the aquifer system. Figure 5(a) indicates that a greater creep time factor (i.e., a larger viscosity coefficient) brings about a slower settling of the aquitard after the decline of the groundwater level. This phenomenon happens because of the viscoelastic characteristics of the aquitard stratum. In general, a larger viscosity coefficient delivers a slower deforming/settling process. With different values of the creep time factor, the differences in the settlement curves achieve a maximum at a given time. However, the differences reduce gradually with the increases in time. As shown in Figure 5(b), the curves with different values of the creep time factor intersect with each other with the increase in time. In other words, for the viscoelastic aquitard, a larger viscosity coefficient initially accelerates the dissipation of the pore-water pressure and further induces a larger degree of consolidation; but at a later stage, a larger viscosity coefficient decelerates the dissipation of the pore-water pressure and further results in a smaller degree of consolidation. This phenomenon is consistent with that found in the consolidation of the viscoelastic soil subjected to external loading [19], and it happens as a result of the mechanical properties of the fractional-derivative Merchant model. As presented in Figure 1, this model is made of an independent spring in a series with a fractional-derivative Voigt body. The stress acting on the independent spring is equal to the total of the stress components undertaken by the fractional-derivative Voigt body, so it represents the vertical effective stress in the aquitard. After withdrawing groundwater, the deformation attained by the aquitard is very small at the initial time, as shown in Figure 5(a), so that its influence on the vertical effective stress is insignificant and can even be neglected. As a result, the vertical effective stress mainly comes from the fractional-derivative element. For this element, a larger viscosity coefficient generally leads to greater vertical effective stress and then brings about a more rapid dissipation of the pore-water pressure throughout the aquitard. With the consolidation progressing, the deformation in the aquitard gradually turns out to be remarkable at a later stage, as shown in Figure 5(a); thus, the effects of the deformation upon the vertical effective stress becomes significant. Under this scenario, the vertical effective stress is mainly controlled by the spring in the Voigt body. Figure 5(a) shows that the greater is the creep time factor (i.e., the viscosity coefficient), the smaller is the deformation of the aquitard. As a result, the greater creep time factor results in a smaller vertical effective stress and then decelerates the dissipation of the pore-water pressure throughout the aquitard. However, the impact of the creep time factor (i.e., the viscosity coefficient) disappears gradually with the increase in time.

Figure 6 presents the influences of the drawdown speed of groundwater upon the settlement of the aquitard and the degree of consolidation. The dimensionless parameter varies from 0.5 to 10. It can be observed that the drawdown speed of groundwater has great influences on the settlement of the aquitard and the degree of consolidation. An increase in the drawdown speed of groundwater accelerates the processes of the pore-water pressure dissipating and the aquitard settling. In other words, the longer the time spent for the same drawdown of groundwater, the later the land subsidence occurs. In addition, these figures indicate that both the degree of consolidation and the aquitard settlement tend to be the same regardless of the drawdown speed of groundwater at last. That is, the drawdown speed of groundwater does not influence the degree of consolidation and the final subsidence of the aquitard.

#### 5. Conclusions

Land subsidence induced by the withdrawal of groundwater is a severe problem in many countries and regions all over the world. In order to simulate the process of land subsidence resulting from the withdrawal of groundwater, this study introduces the fractional-derivative Merchant model to establish a more general model for the aquitard consolidation of an aquifer system. With the use of finite sine transform and Laplace transform, the semianalytical solutions for the pore-water pressure, the settlement of the aquitard, and the degree of consolidation are developed. The comparisons of two special cases are presented to validate the correctness of the developed solution. The consolidation characteristics of the proposed aquifer system are simulated. The results indicated that the proposed aquifer system can be degenerated into an aquifer system with a linear elastic aquitard and an aquifer system with a viscoelastic aquitard obeying the standard Merchant model. The derived semianalytical solution is identical to the available analytical solutions, and it is more common to investigate the consolidation behavior of the aquifer system. The settlement of the aquitard and the degree of consolidation resulting from the drawdown of groundwater increase gradually with the increase in time. The parameters of the fractional-derivative Merchant model as well as the drawdown speed of groundwater have a great effect on the consolidation behavior of the aquitard. After the drawdown of groundwater, an increase in the fractional order initially decreases the aquitard settlement and increases the degree of consolidation, while it increases the aquitard settlement and reduces the degree of consolidation at a later stage. The greater modulus ratio means that the softer aquitard leads to a slower dissipation of the pore-water pressure as well as a more rapid settling. For the larger viscosity coefficient, the settlement of the aquitard progresses more slowly after a decline of the groundwater level, but the dissipation of the pore-water pressure is initially faster while becoming slower at a later stage. The longer time spent for the same drawdown of groundwater, the later the land subsidence occurs. The conclusions of this study may be helpful in understanding the process of land subsidence resulting from the withdrawal of groundwater. A coupled or three-dimensional aquifer system incorporating a fractional-derivative model will be developed to simulate the process of land subsidence induced by the withdrawal of groundwater in a future work.

#### Appendix

#### The Finite Sine Transform and Laplace Transform

The finite sine transform of a real function is defined by [45, 46] where . The corresponding inverse transform is given by

One of the elemental properties of the finite sine transform is

Let the function be real valued and piecewise continuous on and of exponential order (i.e., ), then the Laplace transform of is defined by [47]

The inverse Laplace transform is given by the well-known Bromwich integral: where and can be any real number greater than . Using the idea of the fast Fourier transform, the approximation of equation (A.5) can be expressed as the summation of the Fourier series. Crump’s method uses the following summation formula [42, 51]: where is a real parameter satisfying the condition of . In general, the convergence of the series on the right side in equation (A.6) is slow. For accelerating the convergence, Crump adopted the following epsilon algorithm [47]: with

Equation (A.7) gives a sequence of successive approximations that often better approximates the sum of the series on the right side in equation (A.6) than the untransformed one.

Suppose that the inverse Laplace transform is desired over a range of for which the largest is and the relative error is to be smaller than ; then, the two parameters in equation (A.6) can be selected as follows [42, 47]: and , in which is simply chosen to be a number slightly larger than the real part of the singularity of that has the largest real part.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was partially supported by the National Key Research and Development Plan of China (2016YFC0800200), the National Natural Science Foundation of China (51508180), the China Postdoctoral Science Foundation (2015M570677), and the Natural Science Foundation of Hunan Province (2019JJ50055), China.

#### References

- D. L. Galloway and T. J. Burbey, “Review. Regional land subsidence accompanying groundwater extraction,”
*Hydrogeology Journal*, vol. 19, no. 8, pp. 1459–1486, 2011. View at: Publisher Site | Google Scholar - H. Li and G. Guo, “Surface subsidence control mechanism and effect evaluation of gangue-backfilling mining: a case study in China,”
*Geofluids*, vol. 2018, Article ID 2785739, 9 pages, 2018. View at: Publisher Site | Google Scholar - P. Ezquerro, C. Guardiola-Albert, G. Herrera, J. A. Fernández-Merodo, M. Béjar-Pizarro, and R. Bonì, “Groundwater and subsidence modeling combining geological and multi-satellite SAR data over the Alto Guadalentín Aquifer (SE Spain),”
*Geofluids*, vol. 2017, Article ID 1359325, 17 pages, 2017. View at: Publisher Site | Google Scholar - A. P. S. Selvadurai and J. Kim, “Ground subsidence due to uniform fluid extraction over a circular region within an aquifer,”
*Advances in Water Resources*, vol. 78, pp. 50–59, 2015. View at: Publisher Site | Google Scholar - J. C. Liu, G. G. Lei, and G. X. Mei, “One-dimensional consolidation of visco-elastic aquitard due to withdrawal of deep-groundwater,”
*Journal of Central South University*, vol. 19, no. 1, pp. 282–286, 2012. View at: Publisher Site | Google Scholar - Y. Wang, M. S. Zhang, F. S. Hu, Y. Dong, and K. Yu, “A coupled one-dimensional numerical simulation of the land subsidence process in a multilayer aquifer system due to hydraulic head variation in the pumped layer,”
*Geofluids*, vol. 2018, Article ID 4083439, 12 pages, 2018. View at: Publisher Site | Google Scholar - A. Ortega-Guerrero, J. A. Cherry, and D. L. Rudolph, “Large-scale aquitard consolidation near Mexico City,”
*Ground Water*, vol. 31, no. 5, pp. 708–718, 1993. View at: Publisher Site | Google Scholar - B. Zapata-Norberto, E. Morales-Casique, and G. S. Herrera, “Nonlinear consolidation in randomly heterogeneous highly compressible aquitards,”
*Hydrogeology Journal*, vol. 26, no. 3, pp. 755–769, 2018. View at: Publisher Site | Google Scholar - K. Terzaghi,
*Theoretical Soil Mechanics*, John Wiley, New York, 1943. View at: Publisher Site - M. A. Biot, “General theory of three-dimensional consolidation,”
*Journal of Applied Physics*, vol. 12, no. 2, pp. 155–164, 1941. View at: Publisher Site | Google Scholar - C. E. Jacob, “Radial flow in a leaky artesian aquifer,”
*Transactions American Geophysical Union*, vol. 27, no. 2, pp. 198–205, 1946. View at: Publisher Site | Google Scholar - M. S. Hantush and C. E. Jacob, “Non-steady radial flow in an infinite leaky aquifer,”
*Transactions, American Geophysical Union*, vol. 36, no. 1, pp. 95–100, 1955. View at: Publisher Site | Google Scholar - M. S. Hantush, “Non-steady flow to flowing wells in leaky aquifers,”
*Journal of Geophysical Research*, vol. 64, no. 8, pp. 1043–1052, 1959. View at: Publisher Site | Google Scholar - M. S. Hantush, “Modification of the theory of leaky aquifers,”
*Journal of Geophysical Research*, vol. 65, no. 11, pp. 3713–3725, 1960. View at: Publisher Site | Google Scholar - J. Bear and M. Y. Corapcioglu, “Mathematical model for regional land subsidence due to pumping: 1. Integrated aquifer subsidence equations based on vertical displacement only,”
*Water Resources Research*, vol. 17, no. 4, pp. 937–946, 1981. View at: Publisher Site | Google Scholar - F. Y. Ren and H. Y. Su, “One-dimensional consolidation of aquitard due to water extracting and filling,”
*Shanghai Geology*, vol. 2, pp. 44–55, 1982. View at: Google Scholar - G. Luo, H. Pan, H. Cao, and X. L. Yin, “Analysis of settlements caused by decompression of confined water,”
*Rock and Soil Mechanics*, vol. 25, no. 2, pp. 196–200, 2004. View at: Google Scholar - J. Zeng and X. D. Wang, “Analytical solution to one-dimensional consolidation of aquitards for unsteady flow in confined aquifers,”
*Journal of Nanjing Technology*, vol. 34, no. 2, pp. 85–90, 2012. View at: Google Scholar - K. H. Xie, L. W. Tao, Y. L. Wang, and C. X. Li, “One-dimensional consolidation solution and analysis for aquitard in leakage system,”
*Journal of Shenyang University of Technology*, vol. 34, no. 5, pp. 581–585, 2012. View at: Google Scholar - Z. Li, Z. Zhou, M. Li, B. Zhang, and B. Dai, “Delayed drainage of a largely deformed aquitard due to abrupt water head decline in adjacent aquifer,”
*Geofluids*, vol. 2018, Article ID 2326491, 12 pages, 2018. View at: Publisher Site | Google Scholar - A. F. Hu, C. Q. Xia, H. Wu, K. H. Xie, and L. H. Yan, “A study on one-dimensional consolidation of layered structured aquitard soils in a leakage system,”
*Marine Georesources & Geotechnology*, vol. 35, no. 3, pp. 318–329, 2016. View at: Publisher Site | Google Scholar - K.-H. Xie, X.-Y. Xie, and X.-B. Li, “Analytical theory for one-dimensional consolidation of clayey soils exhibiting rheological characteristics under time-dependent loading,”
*International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 32, no. 14, pp. 1833–1855, 2008. View at: Publisher Site | Google Scholar - J. H. Qian, S. Y. Wang, and Z. P. Guo, “Rheology and soft clay engineering,” in
*STP923-EB Marine Geotechnology and Nearshore/Offshore Structures*, R. C. Chaney and H. Y. Fang, Eds., pp. 112–121, ASTM International, West Conshohocken, PA, USA, 1986. View at: Publisher Site | Google Scholar - Y. Q. Cai, C. J. Xu, and H. M. Yuan, “One-dimensional consolidation of layered and visco-elastic soils under arbitrary loading,”
*Applied Mathematics and Mechanics*, vol. 22, no. 3, pp. 353–360, 2001. View at: Publisher Site | Google Scholar - J. C. Liu, G. H. Lei, and X. D. Wang, “One-dimensional consolidation of visco-elastic marine clay under depth-varying and time-dependent load,”
*Marine Georesources & Geotechnology*, vol. 33, no. 4, pp. 337–347, 2015. View at: Publisher Site | Google Scholar - X. W. Liu, K. H. Xie, Q. Y. Pan, and G. X. Zeng, “Viscoelastic consolidation theories of soils with vertical drain well,”
*Chinese Civil Engineering Journal*, vol. 31, no. 1, pp. 10–19, 1998. View at: Google Scholar - T. W. Hsu and H. J. Liu, “Consolidation for radial drainage under time-dependent loading,”
*Journal of Geotechnical and Geoenvironmental Engineering*, vol. 139, no. 12, pp. 2096–2103, 2013. View at: Publisher Site | Google Scholar - J. F. Miao and L. G. Wu, “Mathematical modelling of land subsidence due to pumping of a multi-aquifer system with viscoelastic properties,” in
*Land Subsidence Proceedings of the Fourth International Symposium on Land Subsidence*, vol. 200, pp. 593–602, Houston, Texas, May 1991, IAHS Publ. View at: Google Scholar - T. L. Tsai, “A coupled one-dimensional viscoelastic-plastic model for aquitard consolidation caused by hydraulic head variations in aquifers,”
*Hydrological Processes*, vol. 29, no. 22, pp. 4779–4793, 2015. View at: Publisher Site | Google Scholar - Y. Zhang, Y. Xue, J. Wu, and Z. Wang, “Compaction of aquifer units under complex patterns of changing groundwater level,”
*Environmental Earth Sciences*, vol. 73, no. 4, pp. 1537–1544, 2015. View at: Publisher Site | Google Scholar - Y. Zhang, Y. Xue, J. Wu, H. Wang, and J. He, “Mechanical modeling of aquifer sands under long-term groundwater withdrawal,”
*Engineering Geology*, vol. 125, pp. 74–80, 2012. View at: Publisher Site | Google Scholar - I. Podlubny,
*Fractional Differential Equations*, Academic Press, California, 1999. - J. Z. Zhang, X. P. Zhou, and P. Yin, “Visco-plastic deformation analysis of rock tunnels based on fractional derivatives,”
*Tunnelling and Underground Space Technology*, vol. 85, pp. 209–219, 2019. View at: Publisher Site | Google Scholar - W. Smit and H. de Vries, “Rheological models containing fractional derivatives,”
*Rheologica Acta*, vol. 9, no. 4, pp. 525–534, 1970. View at: Publisher Site | Google Scholar - W. Zhao, S. Yang, G. Wen, and X. Ren, “Fractional-order visco-plastic constitutive model for uniaxial ratcheting behaviors,”
*Applied Mathematics and Mechanics*, vol. 40, no. 1, pp. 49–62, 2019. View at: Publisher Site | Google Scholar - R. C. Koeller, “Applications of fractional calculus to the theory of viscoelasticity,”
*Journal of Applied Mechanics*, vol. 51, no. 2, p. 299, 1984. View at: Publisher Site | Google Scholar - A. Gemant, “Applications of fractional calculus to the theory of viscoelasticity,”
*Journal of Applied Physics*, vol. 7, no. 8, pp. 311–317, 1936. View at: Google Scholar - D. Yin, H. Wu, C. Cheng, and Y. Chen, “Fractional order constitutive model of geomaterials under the condition of triaxial test,”
*International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 37, no. 8, pp. 961–972, 2013. View at: Publisher Site | Google Scholar - F. Hou, Q. Li, E. Liu et al., “A fractional creep constitutive model for frozen soil in consideration of the strengthening and weakening effects,”
*Advances in Materials Science and Engineering*, vol. 2016, Article ID 5740292, 12 pages, 2016. View at: Publisher Site | Google Scholar - C. Zhang, H. Zhu, B. Shi, and L. Liu, “Theoretical investigation of interaction between a rectangular plate and fractional viscoelastic foundation,”
*Journal of Rock Mechanics and Geotechnical Engineering*, vol. 6, no. 4, pp. 373–379, 2014. View at: Publisher Site | Google Scholar - H.-H. Zhu, C.-C. Zhang, G.-X. Mei, B. Shi, and L. Gao, “Prediction of one-dimensional compression behavior of Nansha clay using fractional derivatives,”
*Marine Georesources & Geotechnology*, vol. 35, no. 5, pp. 688–697, 2016. View at: Publisher Site | Google Scholar - L. Wang, D. Sun, P. Li, and Y. Xie, “Semi-analytical solution for one-dimensional consolidation of fractional derivative viscoelastic saturated soils,”
*Computers and Geotechnics*, vol. 83, pp. 30–39, 2017. View at: Publisher Site | Google Scholar - M. H. Huang and J. C. Li, “Consolidation of viscoelastic soil by vertical drains incorporating fractional-derivative model and time-dependent loading,”
*International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 43, no. 1, pp. 239–256, 2018. View at: Publisher Site | Google Scholar - S. P. Neuman and P. A. Witherspoon, “Theory of flow in a confined two aquifer system,”
*Water Resources Research*, vol. 5, no. 4, pp. 803–816, 1969. View at: Publisher Site | Google Scholar - L. C. Andrews and B. K. Shivamoggi,
*Integral Transforms for Engineers*, SPIE, Bellingham, 1999. View at: Publisher Site - L. Debnath and D. Bhatta,
*Integral Transforms and Their Applications, Second Edition*, Chapman and Hall/CRC, Boca Raton, 2006. View at: Publisher Site - K. S. Crump, “Numerical inversion of Laplace transforms using a Fourier series approximation,”
*Journal of the Association for Computing Machinery*, vol. 23, no. 1, pp. 89–96, 1976. View at: Publisher Site | Google Scholar - P. Ni, K. Xu, G. Mei, and Y. Zhao, “Effect of vacuum removal on consolidation settlement under a combined vacuum and surcharge preloading,”
*Geotextiles and Geomembranes*, vol. 47, no. 1, pp. 12–22, 2019. View at: Publisher Site | Google Scholar - B. C. Hawlader, B. Muhunthan, and G. Imai, “Viscosity effects on one-dimensional consolidation of clay,”
*International Journal of Geomechanics*, vol. 3, no. 1, pp. 99–110, 2003. View at: Publisher Site | Google Scholar - P. Ni, G. Mei, and Y. Zhao, “Surcharge preloading consolidation of reclaimed land with distributed sand caps,”
*Marine Georesources & Geotechnology*, pp. 1–12, 2018. View at: Publisher Site | Google Scholar - F. Durbin, “Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method,”
*The Computer Journal*, vol. 17, no. 4, pp. 371–376, 1974. View at: Publisher Site | Google Scholar

#### Copyright

Copyright © 2019 Ming-hua Huang and Dun Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.