Probabilistic analysis of land subsidence due to pumping by Biot poroelasticity and random field theory

Land subsidence is a global problem in urban areas. The main cause of land subsidence is the pumping of subsurface water. It is of great significance to study the subsurface settlement and water flow of the lands due to pumping. In this study, the probabilistic analysis of land subsidence due to pumping is performed by Biot’s poroelasticity and random field theory based on a case study. The results show that the change of deformation of the aquifer is far less significant than the hydraulic head over the years. When considering the spatial variability of soil strength, the land subsidence suffers from great uncertainty when the correlation length is large. Nevertheless, the spatial variability of soil strength on the uncertainty of hydraulic head can be ignored. When considering the spatial variability of soil hydraulic conductivity, the uncertainty of the hydraulic head is mainly located near the bedrock and increases markedly along with the rise of the correlation length. Time is another important factor to increase the uncertainty of the hydraulic head. However, its contribution to the uncertainty of displacement is insignificant.


Introduction
Land subsidence is the gradual or rapid sinking of the ground surface due to the deformation of subsurface earth materials, which is a global problem in urban areas [18,43]. The main cause of land subsidence is the pumping of subsurface water [8,14,31]. It is of great significance to study the subsurface settlement and water flow of the lands due to pumping.
The land subsidence is often simulated or evaluated based on the soil consolidation theory, which is a process of volumetric changes of soil due to water pressure. Early methods to model soil consolidation are based on Terzaghi's theory. It assumes that the settlement and flow of water are vertical. Ignoring the horizontal deformation does not allow for a complete analysis of problems of consolidation. If the horizontal deformation needs to be considered, this one-dimensional theory of consolidation may not be valid. In recent decades, the more rigorous Biot's poroelasticity considering horizontal and vertical components of elastic deformation has been widely used for the problems of land subsidence. Bear and Corapcioglu [3] developed a mathematical model for regional subsidence due to pumping from an aquifer based on Biot's theory on coupled three-dimensional consolidation. Chiou and Chi [11] studied the settlement induced by surface loading and land subsidence due to pumping for saturated layered soils. Xu et al. [42] presented the prediction approaches on land subsidence employed in China and found that Biot's consolidation can simulate the field data better. Ferronato et al. [16] proposed a coupled Biot model based on a three-field formulation to predict the land subsidence in the Chaobai River alluvial fan, China.
However, most of the studies related to land subsidence did not consider the uncertainty of geo-properties. It is well recognized that the subsurface geo-properties such as seepage and strength parameters are remarkably variable and heterogeneously suffering from great uncertainty. To understand the uncertainty of soil consolidation, probabilistic analysis by Monte Carlo simulation is always adopted regarding the different engineering geological backgrounds. The parameters in Biot's formulations are modeled as random variables to account for the uncertainty of subsurface geo-properties or further modeled as random fields to consider spatial variability. For example, Houmadi et al. [23] used a collocation-based stochastic response surface method for the probabilistic analysis of a consolidation problem of a single clayey layer, and the deterministic model is based on a Biot consolidation analysis using the finite difference code FLAC 3D. Cheng et al. [10] integrated random field simulation of soil spatial variability with numerical modeling of coupled flow and deformation to investigate consolidation in spatially random unsaturated soil. Zhang et al. [49] proposed a probabilistic method to calibrate coupled hydro-mechanical slope stability model with the integration of multiple types of field data. Houmadi et al. [24] analyzed the impact on surface settlement due to a uniform surcharge loading on the ground surface with a twodimensional spatially varying Young's modulus by the subset simulation method. Savvides and Papadrakakis [30] presented a stochastic analysis to study the consolidation phenomenon of clayey interaction. In summary, based on the models of Biot's consolidation, the uncertainty of many geotechnical issues including land reclamation, embankments, tunnels, and excavation are evaluated by several researchers. However, the probabilistic analysis for the problem of land subsidence is seldom involved. Therefore, in this study, the probabilistic analysis of land subsidence due to pumping is performed by Biot's poroelasticity and random field theory. First, based on Leake and Hsieh [26], the numerical model of an aquifer underlain by a bedrock step and pumping is established. Second, to consider soil spatial variability, two key parameters (i.e., Young's modulus and hydraulic conductivity) in Biot's equations are viewed as heterogeneous properties and generated by random field theory. Finally, the influence of correlation length and time on the uncertainty of pumping responses (i.e., displacement and hydraulic head) are investigated.

Biot's poroelasticity
In this study, the built-in module in COMSOL Multiphysics [13] is adopted to simulate land subsidence. Based on Biot's poroelastic theory [4,5], the constitutive relations for the poroelastic behavior are: where σ is the total stress; ":" stands for the double-dot tensor product; c denotes the elasticity matrix of solid; ε is the strain tensor; p is the fluid pore pressure; I is the identity matrix; α b is the Biot-Willis coefficient representing the coupling between the stress and the pore pressure. The value of α b is less than unity, indicating the extent to which the pore pressure contributes on elastic deformation.
The form of force balance equation is: where ρ represents the average density of solid and fluid; ρ f and ρ s are the density of the fluid and solid, respectively; ϕ is the porosity; g represents the acceleration of gravity. Note that Eq. (1) is the linear theory of elasticity, implying that the general theory proposed by Biot is the linear poroelasticity. Biot's equations can be extended to nonlinear poroelasticity, such as elastoplastic materials, by changing the form of Eq. (1) [2].
Based on the mass conservation equation, with the increase of the rate of expansion of the pore space, the volume available for the fluid also increases and thereby gives rise to liquid sink [22]: where t is time; k is the hydraulic conductivity; ε v is the volumetric strain, ε v = ε x + ε y +ε z , which is the trace of ε; S b is the storage coefficient of Biot's poroelasticity, which is related to the compressibility of the fluid and solid phases. When both the solid and the fluid are assumed compressible, it can be calculated from basic material properties as [7]: where K f is the fluid bulk modulus, which is the inverse of the fluid compressibility χ f , and K s is the solid bulk modulus and K s ¼ E 3ð1−2νÞ for elastic materials. E and ν are Young's modulus and Poisson's ratio, respectively.
For saturated soil, some studies assumed that the water and soil are incompressible. Therefore, the values of S b and α b can be 0 and 1, respectively, and ρ is equal to the density of soil [6,[23][24][25]34]. While some other studies, such as oil reservoir simulation, considered the contributions of S b and α b [20,47]. Since the poroelasticity of Biot's consolidathangion is a built-in module in COMSOL Multiphysics, the solution of the above equations is very convenient. As a result, the compressible nature of soil and water is taken into consideration in this study.

Numerical model of an aquifer
The numerical model of land subsidence is referenced from Leake and Hsieh [26]. There is an aquifer system overlying an impermeable bedrock in a basin. The height of the aquifer is 420 m, and the length exceeds 4000 m. The bedrock is a fault and acts as a step near a mountain front. The aquifer system includes a middle compressible confining unit, which is 20 m below the ground surface ( Fig. 1).
In this study, the predefined mesh grid of COMSOL is adopted. The finer element size (Fig. 1b) is chosen for simulation. The maximum element size is 117 m. Note that the finite element mesh is usually finer than the random field grid to capture the information on the spatial variability. However, the overly fine mesh will lead to high computational costs. There is a trade-off between the accuracy of the solution and computational efficiency. In this study, the maximum element size is larger than some examined correlation lengths because the area to be simulated is very large. To overcome this problem, the midpoint discretization method [32,37] is employed to determine grid points of the random field. Shen et al. [33] illustrated that this method is sufficient to obtain accurate statistics of model responses. Please refer to Shen et al. [33] for the discussion of finite element meshes and discretization error.
For the deterministic model, the parameters of an aquifer, semi-confined layer, and water are summarized in Table 1. The hydraulic and physical properties are set as the alluvial basin in the southwestern USA [21]. The values of porosity for aquifer ϕ a and semi-confined layer ϕ i are 0.25 and 0.025, respectively. The hydraulic conductivity k a of the aquifer is 25 m/day whereas k i = 0.01 m/day for the semi-confined layer. Young's modulus is assumed to be different. E a = 800 MPa for aquifer and E i = 80 MPa for semi-confined layer. Except for the above parameters, the Poisson's ratio and density of soil are the same for the aquifer and the semi-confined layer. The Poisson's ratio ν and ρ are assumed to be 0.25 and 2750 kg/m 3 , respectively. The constants for the water of compressibility χ f and density ρ f are 4 × 10 −10 1/Pa and 1000 kg/m 3 , respectively.
The boundary conditions of the aquifer model are shown in Fig. 1. The hydraulic head in JA-AB-BC is specified as zero-constant during the entire period of simulation to assume that no consolidation occurs in this part. IH is fixed with a head that linearly declines by 60 m over 10 years. Other boundaries are no-flow. For the mechanical

Random field
It is well known that the soil properties of an aquifer are variant but correlated in space due to the geological processes. Site investigation can only obtain limited samples of soil parameters. From the point of view of probability, the statistical characteristics of soil parameters can be obtained from limited samples with randomness. Therefore, random field theory is used to characterize the spatial variability of soil properties. Soil parameters such as Young's modulus and hydraulic conductivity are positive and fit well with log-normal distributions [1,29,48]. Therefore, the natural logarithm of a certain soil parameter follows a normal distribution. Its mean value μ ln and the standard deviation σ ln are calculated as follows: where μ and σ are the mean value and the standard deviation of soil parameters, respectively.
In random field theory, the covariance function is proposed to illustrate the spatial correlation of a certain soil parameter. It is a function related to coordinates x = [(x 1 , z 1 ), (x 2 , z 2 )] in the domain. The horizontal and vertical correlation lengths (l x and l z ) are thresholds to determine the relevance of a soil parameter of two positions in the domain. In this study, an empirical covariance function C(x) is used to simulate the spatial variability of soil parameters [40,41,44]:   To generate random fields, the covariance function C(x) is decomposed by the Karhunen-Loève expansion method as previous studies [45,46]. More details of this method can be found in Ghanem and Spanos [19].
The land subsidence based on Biot's consolidation is a coupled hydro-mechanical problem. Therefore, two parameters, E a and k a , of strength and hydraulic conductivity for the aquifer are modeled by random field to consider their spatial variability. Correspondingly, two cases are used to illustrate the effects of spatial variability of soil strength and hydraulic conductivity on the uncertainty of model responses. It is recognized that the hydraulic properties of soil suffer from great uncertainty. According to previous studies, the CoV of the saturated coefficient of hydraulic conductivity can be ranged from 50 to 450% [9,48]. Relatively, the CoVs of soil strength parameters are small, around 5~50% [12,28]. Therefore, the CoVs of k a and E a are assumed to be 80% and 30% in this study, respectively. The first case considers the spatial variability of soil strength. The mean and coefficient of variation (COV) of E a are 800 MPa and 0.3, respectively. The same idea applies to k a for the second case, where it no longer goes into details. Please refer to Table 2 accordingly.
The selection of the correlation lengths for the parametric study is mainly based on the following facts: (1) For natural soil parameters, the vertical correlation length varies from less than 1 m to more than 20 m [15,17,36]. The horizontal correlation length is generally much larger than the vertical length due to the stratification of natural deposits. (2) For the practice of probabilistic study, many studies set the correlation lengths as a ratio of the model size for uncertainty or reliability analysis [35,50]. It is suggested that the correlation length of the soil parameters can be taken as 0.02~2 times the model size. 3. In geotechnical engineering, the site-scale models are generally adopted, and the model size is commonly less than 100 m, while the model size in this study is in large basin-scale. The large-scale models, such as watershed-scale, are considered as references. It is reported that the correlation length can exceed 650 m [38,39]. Therefore, in this study, l x and l z vary from 200 to approximately 800 m and 40 to approximately 160 m, respectively. Typical realization of lognormal random fields with different correlation lengths is shown in Fig. 2.
The uncertainty of the model responses can be determined by running the model repeatedly with different random soil parameters to arrive at an estimate of the standard deviation of the model responses, i.e., the so-called Monte Carlo simulation. A sensitivity analysis is conducted to determine the number of random fields for Monte Carlo simulation. Figure 3 presents the effect of the number of random fields on the mean values of subsidence at surface nodes. There is almost no fluctuation of the estimation of mean values when the number of random fields is less than 500. Therefore, a total number of 500 random fields is generated to assess the uncertainty of the model responses, which is also consistent with the previous study suggested by Peng et al. [27].

Results and discussions
Deterministic results Figure 4 shows the deterministic results of displacement over the years. The displacement at the upper boundary indicates the surface subsidence. The surface subsidence exceeds 2 m, and it is gradually grown over the years. With the increase of depth, the displacement is decreased and less sensitive to time. Figure 5 shows the deterministic results of hydraulic heads over the years. The hydraulic head in the whole domain is reduced rapidly as a result of pumping. The hydraulic head around the bedrock is reduced from − 4 to − 40 m for 10 years, which nearly drops 4 m per year due to pumping. Comparably, the change of deformation of the aquifer is far less significant than the hydraulic head over these 10 years.

Case 1: Spatial variability of soil strength
The effect of the correlation length of E a on the standard deviation of displacement (σ s ) is illustrated in Fig. 6. With the increase in correlation length of E a , σ s increases dramatically. The maximum value of σ s for surface settlement is around 0.6 m with the largest correlation length. The effect of the correlation length of E a on the standard deviation of displacement is significant. The land subsidence due to pumping suffers from great uncertainty when the correlation length of soil strength properties is large.  Figure 7 shows the standard deviation of the hydraulic head (σ h ) considering the spatial variability of E a . Although the σ h rises with the increase of the correlation length of E, the values of σ h are very small compared to displacement. Even when l x = 800 m and l z = 160 m, the maximum σ h is only 0.004 m. It is indicated that the spatial variability of E a has a slight influence on the uncertainty hydraulic head for land subsidence due to pumping.

Case 2: Spatial variability of soil hydraulic conductivity
The effects of the correlation length of k a on the uncertainty of displacement are displayed in Fig. 8. The spatial variability of soil hydraulic conductivity has a minor influence on the uncertainty of displacement. The maximum σ s is only 0.03 m in year 10. It implies that the uncertainty of displacement is insignificant when dealing with the spatial variability of soil hydraulic conductivity. Figure 9 presents the effect of the correlation length of k a on the uncertainty of the hydraulic head. When l x = 200 m and l z = 40 m changes to l x = 800 m and l z = 160 m, respectively, the maximum value of σ h increases from 2 to 6 m, which is nearly tripled. The correlation length of k a is strongly influential to σ h . The σ h increases markedly along with the rise in the correlation length of k a . In addition, like Fig. 7, the σ h around the bedrock is comparatively large, which illustrates that the uncertainty of the hydraulic head is mainly located near the bedrock. Effects of time Case 1: Spatial variability of soil strength Figure 10 shows the uncertainty of displacement over the years considering the spatial variability of E a . The σ s is steadily increased with years. In the tenth year, the maximum σ s is approximated to 0.30 m near the surface. It illustrates that the uncertainty of land subsidence rises gradually over the years due to pumping. The effects of spatial variability of E a on σ h over the years are shown in Fig. 11. There is no difference among them, indicating the σ h is constant with time. Although the hydraulic head is gradually increased, its uncertainty is invariable over the years if only considering the spatial variability of E a . Besides, the values of σ h are all around the order of a millimeter, indicating the trivial effect on the uncertainty of the hydraulic head. To conclude, the spatial variability of E a on the uncertainty of hydraulic head for land subsidence due to pumping can be ignored.

Case 2: Spatial variability of soil hydraulic conductivity
The uncertainty of displacement over the years considering spatial variability of k a is shown in Fig. 12. Obvious changes in the σ h appeared, but σ h is only approximated to 0.01 m even in the 10-year settlement for land subsidence. The contribution of spatial variability of soil hydraulic conductivity to the uncertainty of displacement is unimportant. In Fig. 13, the uncertainty of hydraulic head over the years considering the spatial variability of k a is shown. In year 1, the maximum of σ h is approximately 0.2 m, and it is increased to 2 m in year 10. Therefore, besides the correlation length of hydraulic conductivity, time is another important factor to increase the uncertainty of hydraulic head for land subsidence due to pumping.

Effects of boundary conditions
The effect of boundary conditions on the uncertainty of land subsidence is further investigated. The results of two different boundary conditions are shown in Fig. 14. Figure 14a shows the hydraulic head and land subsidence at year 1 with hydraulic head boundary condition. Figure 14b shows the corresponding result with flux boundary condition in the final steady state. It can be seen that the two different boundary conditions produce the same results.
The effects of boundary conditions on the uncertainty of hydraulic head considering spatial variability of k a are shown in Fig. 15. When choosing head boundary condition, a large uncertainty appeared around the bedrock (Fig. 15(a)). However, in Fig. 15(b), the uncertainty around the flux boundary condition is large. The standard deviation of the hydraulic head exceeds 0.25 m. Therefore, flow boundary condition has an obvious impact on the uncertainty of the hydraulic head.

Conclusions
In this study, the probabilistic analysis of land subsidence due to pumping is performed by Biot's poroelasticity and random field theory based on a case study. First, the numerical model of an aquifer underlain by a bedrock step and pumping is established. Second, to consider soil spatial variability, two key parameters in Biot's equations controlling deformation and hydraulic head are viewed as heterogeneous properties and generated by random field theory. Finally, the influences of correlation length and  1. The total surface settlement exceeds 2 m for 10 years of land subsidence due to pumping. The hydraulic head around the bedrock nearly drops 4 m per year due to pumping. In general, the change of deformation of the aquifer is far less significant than the hydraulic head over these 10 years. 2. When considering the spatial variability of soil strength, it suffers from great uncertainty when the correlation length is large. The uncertainty of displacement gradually rises over the years. Nevertheless, the spatial variability of Young's module on the uncertainty of hydraulic head can be ignored. 3. When considering the spatial variability of soil hydraulic conductivity, the uncertainty of the hydraulic head is mainly located near the bedrock and increases markedly along with the rise of the correlation length. Time is another important factor to increase the uncertainty of the hydraulic head. However, its contribution to the uncertainty of displacement is insignificant.