Effects of soil parameter variation and grid deformation on the numerical simulation of aquitard consolidation

Changes in soil parameters such as saturated unit weight (body force), hydraulic conductivity, and elastic modulus must be considered for simulating aquitard consolidation. In addition to changes in soil parameters, grid intervals change the subsequent soil deformation in such simulations. This study incorporated and examined the effects of both the grid deformation and soil parameter variation on the numerical simulation of aquitard consolidation caused by hydraulic head decline in a multilayered aquifer system. According to the nondimensional analysis, the deformation number can be regarded as a dimensionless parameter for assessing whether the grid deformation and soil parameter variation affect this simulation. The results revealed that when the deformation number exceeds 0.01 (i.e., when the aquitard is thick or soft), both the grid deformation and soil parameter variation must be simultaneously considered to increase accuracy. Aquitard consolidation was always overestimated when only the soil parameter variation was considered. However, when only the grid deformation was considered, the aquitard consolidation was either overestimated or underestimated depending on whether the influence of the change in the body force was more substantial than that of the elastic modulus. Article history: Received 5 July 2017 Revised 28 September 2017 Accepted 1 November 2017


IntroductIon
The rigorous theory of soil consolidation originated from the effective stress principle proposed by Terzaghi (1925).He hypothesized that the total stress of saturated soil remains constant, soil deformation occurs only slightly in the vertical dimension, and fluid is incompressible and satisfies Darcy's law.With these hypotheses, he developed the widely used one-dimensional consolidation equation.Jacob (1950) hypothesized that fluid is compressible, and combined the conventional two and three-dimensional groundwater flow equation (Theis 1938) with Terzaghi's one-dimensional consolidation theory.Thus, the Terzaghi-Jacob two-and three-dimensional consolidation groundwater flow equation was formulated, and the physical significance of the storage coefficient was clarified.Following the development of the mass conservation relationships for soil, water, and air, and the static equilibrium for solid-phase soil, Verruijt (1969) used the concept of multiphase flow to systematically estab-lish a three-dimensional poroelastic consolidation model for the consolidation of unsaturated soil.Under the hypotheses that the total stress in soil remains constant and that soil is saturated and consolidated only in the vertical direction, the groundwater flow equation proposed by Jacob (1950) is a special case of the three-dimensional poroelastic consolidation equation proposed by Verruijt (1969).Other researchers such as Lewis and Schrefler (1978), Safai and Pinder (1980), Bear and Corapcioglu (1981), Lewis et al. (1991), Tarn and Lu (1991), Yeh et al. (1996), Gambolati et al. (2000), and Ye et al. (2016) have investigated the problems of soil deformation by empolying the poroelastic consolidation model.
When soil is consolidated, the resulting deformation changes its porosity and causes variation in the soil parameters, for example, saturated unit weight (body force), hydraulic conductivity, and elastic modulus.On the basis of poroelastic theory, Mei (1985) and Fallou et al. (1992) have investigated the body force effect of the porosity changes on soil consolidation due to surface loading.Tsai et al. (2006) and Tseng et al. (2008) have examined the body force effect Terr. Atmos. Ocean. Sci., Vol. 29, No. 4, 455-464, August 2018 on aquitard consolidation caused by hydraulic head decline in multilayered aquifer systems.Wang and Hsu (2009) analyzed the effect of changes in hydraulic conductivity and elastic modulus on soil consolidation due to surface loading.Tsai and Jang (2014) further examined the combined influences of variations in body force, hydraulic conductivity, and elastic modulus on aquitard consolidation.
When aquitard consolidation is simulated numerically, the soil deformation causes changes not only in the body force, hydraulic conductivity, and elastic modulus, but also in the grid intervals.Therefore, all such changes must be simultaneously considered in the numerical simulation of aquitard consolidation.The goal of this study was to examine the influence of the changes in grid intervals and soil parameters on the numerical simulation of aquitard consolidation caused by hydraulic head decline in a multilayered aquifer system.

GovErnInG EquAtIons
Neglecting the compressibility of the fluid and soil solids and combing the mass conservation of the fluid and soil yield the groundwater flow equation of a saturated porous medium satisfying Darcy's law (Verruijt 1969) as follows: where K represents the hydraulic conductivity tensor, u is the soil displacement vector, and t represents time.h denotes the hydraulic head.
To calculate the soil deformation caused by the changes in the hydraulic head during the consolidation process, the incremental static equilibrium equation of a saturated porous medium using the elastic constitutive relation (Verruijt 1969;Tsai and Jang 2014) can be expressed as follows: where G and m are Lamé constants, n 0 represents the initial porosity.w t and s t denote the fluid density and soil density, respectively.h e = hh 0 is the change in the hydraulic head, and h 0 is the initial hydraulic head.g is the gravitational acceleration.The rightmost term on the left-hand side of Eq. (2) represents the body force effect (Tsai et al. 2006).
Sandy soil exhibits higher permeability and lower compressibility than clay.The concept of quasi-three-dimensional flow can be used to simplify the groundwater flow in a multilayered aquifer system.The groundwater flow in the aquifer is horizontal whereas that in the aquitard is vertical.The soil displacement in the aquitard is assumed to occur only vertically.In the multilayered aquifer system shown in Fig. 1, according to Eqs. ( 1) and ( 2), the groundwater flow equation and incremental static equilibrium equation in the aquitard can be respectively simplified as follows: where K and u z represent the hydraulic conductivity and soil displacement in the vertical direction, respectively, and G 2 s a m = + ^h represents the elastic modulus of the soil.By using the Kozeny-Carman equation (Bear 1972), Wang and Hsu (2009) proposed the following relational equation between hydraulic conductivity and soil deformation, where K 0 represents the initial hydraulic conductivity and u z z 2 2 f = denotes the soil deformation rate per unit volume.Tsai and Jang (2014) derived the following relationship between the elastic modulus and soil deformation, where s 0 a represents the initial elastic modulus.According to Eqs. (3) and ( 4), an aquitard consolidation equation can be established by considering the body force effect.The changes in the hydraulic conductivity and elastic modulus due to soil deformation can be further incorporated on the basis of Eqs. ( 5) and (6).To solve the Fig. 1.Schematic illustration of multilayered aquifer system.aforementioned aquitard consolidation equation, the initial and boundary conditions are required.The initial conditions of the hydraulic head and soil displacement are respectively expressed as follows: , () Where h 0 (z) represents the initial hydraulic head.The hydraulic head of the interface between the aquitard and an aquifer must be continuous.The upper and lower boundary conditions of the hydraulic head are respectively expressed as follows: , () where h u (t) and h l (t) represent the hydraulic heads of the upper and lower aquifers, respectively, with respect to time.
The lower boundary condition for the soil displacement is expressed as follows:

dIscrEtIzInG thE GovErnInG EquAtIons
The groundwater flow and incremental static equilibrium equations of the aquitard are expressed as Eqs.( 3) and (4), respectively.Using the implicit central difference method, the two aforementioned equations can be discretized at the time n + 1 and at the grid point i as follows (Fig. 2): where Δt and ( ) D + represent the time and grid intervals, respectively; i = 1, 2, …, M, where 1 and M are the grid numbers of the upper and lower boundaries, respectively; and n = 0, 1, …, NT, where 0 and NT represent the initial and final simulated time, respectively.Equations ( 13) and ( 14) can be further expressed as follows: f + , and elastic modulus ( ) a + are respectively expressed as follows: In addition, the initial conditions of the hydraulic head and soil displacement are respectively discretized as follows: where ΔZ represents the initial grid interval.The upper and lower boundary conditions of the hydraulic head are respectively discretized as follows: The lower and upper boundary conditions of the soil displacement are respectively discretized as follows:

cAlculAtInG GrId dEformAtIon And solutIon ProcEdurE
The effect of grid deformation on the numerical simulation of aquitard consolidation was considered in this study.When the aquitard is consolidated, the resulting soil displacement alters the grid intervals.The relationship between the grid intervals and soil displacement is expressed as follows: The soil displacement and hydraulic head can be numerically solved based on an iterative coupling method using the discretized governing equations shown in Eqs. ( 15) -( 16) together with the soil parameter variations, boundary conditions, initial conditions, and grid deformation shown in Eqs. ( 17) -( 27).The convergent condition in each time step is the absolute relative error between the two iterative values less than 0.001 at each grid point, respectively for the soil displacement and hydraulic head.

EffEcts of GrId dEformAtIon And soIl PArAmEtEr vArIAtIon
Before the effects of the grid deformation and soil parameter variation were examined, the grid independent test was conducted using a multilayered aquifer system with an aquitard between two aquifers as shown in Fig. 1.The initial steady state was assumed and the changes in the hydraulic heads of the upper and lower aquifers were -4 and -3 m, respectively.The hydraulic head declines of the aquifers linearly increased to the maximum values in 6 months and then remained constant.The following values were used in the test: initial thickness B 0 = 15 m, initial elastic modulus s 0 a = 2 × 10 6 N m -2 , initial hydraulic conductivity K 0 = 1 × 10 -3 m day -1 , initial porosity n 0 = 0.3, and specific weight G s = 2.65.The simulation results of the steady-state aquitard consolidation using the various time and grid intervals for three methods are listed in Table 1.Methods 1 and 2 considered only the effects of soil parameter variation and grid deformation, respectively.Method 3 accounted for the two effects combined.Table 1 indicates that the simulation results are identical when the time and grid interval are less than 10 days and 0.05 m, respectively.Therefore, to ensure high accuracy, the time interval Δt = 5 days and the grid interval ΔZ = 0.01 m were used in the following examinations.
Five scenarios involving an aquitard between two aquifers were used to examine the effects of the changes in the grid intervals and soil parameters on the simulation results of steady-state aquitard consolidation.Table 2 lists the initial thickness, maximum hydraulic head decline, initial hydraulic conductivity, and initial elastic modulus for each scenario.The specific weight and initial porosity of the aquitard in all scenarios were 2.65 and 0.3, respectively.The initial steady state was assumed in the examination.The hydraulic head declines of the upper and lower aquifers in each of the scenarios equally and linearly increased to the maximum values in 6 months and then remained constant.As shown in Table 2, Scenario 2 differed from Scenario 1 in that its initial elastic modulus was only one-thirtieth that of Scenario 1.In other words, the aquitard in Scenario 2 was 30 times as soft as that of Scenario 1.The maximum hydraulic head declines in Scenarios 3 and 4 were respectively twice and one-fifth that of Scenario 2. The aquitard in Scenario 5 was initially twice as thick as that of Scenario 4.
Table 3 presents the simulation results of the steadystate aquitard consolidation using the three methods in Scenarios 1 -5.In Scenario 1, all three methods produced identical simulation results.Figures 3 and 4 display the simulation results of the hydraulic head variations and soil displacements in the vertical direction with respect to time in Scenario 1.In the other scenarios, the simulation results of the three methods differed not only when the elastic modulus was small but also when the thickness was increasing.The simulation results revealed that the steady-state aquitard consolidation was always overestimated when only the effect of soil parameter variation was considered.However, when only the effect of grid deformation was incorporated, the steady-state aquitard consolidation could be either underestimated or overestimated.
As shown in Table 3, in Scenarios 2, 3, 4, and 5, the simulation results of the steady-state aquitard consolidations using Method 1 were always greater than those using Method 3. The two methods differed according to whether the effect of grid deformation was considered.The simulated steady-state aquitard consolidation was smaller because the grid intervals decreased following soil consolidation.In other words, the simulated steady-state aquitard consolidation decreased when the effect of grid deformation was considered.In Scenarios 2 and 3, the simulation results of the steady-state aquitard consolidations using Method 2 were greater than those using Method 3; the opposite was true for Scenarios 4 and 5.The two methods differed depending on whether the effects of the variations in soil parameters such as body force and elastic modulus were considered.When the effect of the change in the body force was included, the simulated steady-state aquitard consolidation increased because the saturated unit weight of the soil increased.However, when the effect of the change in the elastic modulus was considered, the simulated steady-state aquitard consolidation decreased because the consolidation enlarged the elastic modulus.Accordingly, the effect of the change in the body force on aquitard consolidation was the reverse of that of the change in the elastic modulus.Therefore, when the effect of the change in the body force was smaller than that of the change in the elastic modulus, Method 2 produced a higher steady-state aquitard consolidation than did Method 3, as revealed in Scenarios 2 and 3. Conversely, Method 2 had a lower steady-state aquitard consolidation than did Method 3 when the effect of the change in the body force was larger than that of the change in the elastic modulus, as shown in Scenarios 4 and 5.
combined effect of both changes was still larger than that of the grid deformation.In addition, in Scenarios 4 and 5, because the effect of the change in the body force was larger than that of the change in the elastic modulus, the combined effect of both changes was also greater than that of the grid deformation.However, Table 3 shows that in Scenario 3, Method 1 produced a higher steady-state aquitard consolidation than did Method 1 not only because the effect of the change in the body force was less than that of the change in the elastic modulus, but also because the combined effect was lower than that of the grid deformation.
Percent relative errors were used to further demonstrate the differences in the simulation results using the three methods as shown in Table 4.A percent relative error is defined as follows: where u Method 3 denotes the simulation obtained using Method 3, which includes the combined effects of soil parameter variation and grid deformation.u Methods 1 or 2 represents the simulation results using either Methods 1 or 2, which includes only the effect of the soil parameter variation or grid deformation.
As shown in Table 4, the percent relative errors involving only the effect of soil parameter variation in Scenarios 2 -5 were all positive.In Scenarios 2 and 3, the percent relative errors involving only the effect of the grid deformation were positive, but those in Scenarios 4 and 5 were negative.As revealed in the comparison of Scenario 1 with the other scenarios, when the elastic modulus was reduced (i.e., the aquitard was softened), the effects of the soil parameter variation and grid deformation on the steady-state aquitard consolidation intensified.In addition, the comparison of Scenarios 4 and 5 revealed that the thicker the aquitard was, the more substantial the effects of the soil parameter variation and grid deformation on steady-state aquitard consolidation became.

nondimensionalized Governing Equations
Let B 0 be the initial aquitard thickness representing the characteristic length, s 0 a denote the characteristic elastic modulus, and K 0 represent the characteristic hydraulic conductivity.The groundwater flow equation for the aquitard can be nondimensionalized as follows: where u u B The nondimensionalized upper and lower boundary conditions of the hydraulic head are respectively expressed as follows: , ( ) The nondimensionalized lower and upper boundary conditions of the soil displacement are respectively written as follows: scenario 1 scenario 2 scenario 3 scenario 4 scenario 5 Method 2 0% 0.96% 5.68% -3.05% -7.00% Method 1 0% 2.53% 4.79% 0.53% 0.57% In addition, the nondimensionalized grid intervals are expressed as follows: where , and

nondimensional Analysis results
A nondimensional steady-state analysis was first performed on a multilayered aquifer system involving an aquitard between two aquifers.The initial steady state was assumed in the multilayered aquifer system.The hydraulic heads of the upper and lower aquifers were consistent and remained unchanged after linearly reaching their maximum declines.The nondimensionalized time interval t * D = 0.001 and the nondimensionalized grid interval Z * D = 0.001 were used based on the grid independent test.Table 5 illustrates the percent relative errors of the nondimensionalized steadystate aquitard consolidation with various initial porosities, DNs, and nondimensionalized maximum hydraulic head declines.
As shown in Table 5, the percent relative errors did not differ significantly with the various initial porosities.When the effect of only the soil parameter variation was considered, all the percent relative errors were positive and increased following the increase in the nondimensionalized maximum hydraulic head decline.However, the percent relative errors subsequently changed from negative to positive when only the effect of grid deformation was considered.When the DN was higher than 0.01, the percent relative errors became significant.In other words, when the DN is higher than 0.01 in the aquitard consolidation simulation the effects of the grid deformation and soil parameter variation must be simultaneously considered.These results from the nondimensional steady-state analysis were consistent with the simulation results from Scenarios 1 -5.The DNs were 0.00245 for Scenario 1, 0.0735 for Scenarios 2 -4, and 0.147 for Scenario 5. Accordingly, the combined effect of the grid deformation and soil parameter variation must be considered in Scenarios 2 -5 to avoid the overestimation or underestimation of aquitard consolidation.
A nondimensional transient analysis was then conducted on the same aquifer system subjected to varying hydraulic head declines as follows: Figure 5 illustrates the nondimensionalized consolidations with respect to time using the aforementioned three methods with the initial porosity of 0.3 and various DNs.
When the DN was higher than 0.01, the simulation results from the three methods differed significantly and the differences increased alongside the increase in the DN.Therefore, the aquitard consolidation simulation must simultaneously consider the effects of both grid deformation and soil parameter variation when the DN is higher than 0.01.

conclusIon
This study investigated the effects of both grid deformation and soil parameter variation on the numerical simulation of the aquitard consolidation caused by hydraulic head decline in a multilayered aquifer system.According to the simulations, the steady-state aquitard consolidation involving only the effect of the soil parameter variation was always higher than that involving the combined effects of soil parameter variation and grid deformation.However, the steadystate aquitard consolidation considering only the effect of grid deformation was either higher or lower than that considering both effects depending on whether the effect of the soil parameter variation was dominated by the change in the body force or elastic modulus.This result indicates that the aquitard consolidation may be overestimated or underestimated when the combined effects of the grid deformation and soil parameter variation are not simultaneously considered.The steady-state and transient results of the nondimensional analysis indicates that the differences between the aquitard consolidation simulation involving the effect of only the grid deformation or soil parameter variation and that involving both effects reached significance when the DN exceeded 0.01.Therefore, in numerical simulations of aquitard consolidation, when the DN is higher than 0.01, the combined effects of the grid deformation and soil parameter variation must be considered.

Fig. 3 .
Fig. 3. Simulation results of hydraulic head variations with respect to time in Scenario 1.

Fig. 4 .
Fig. 4. Simulation results of soil displacement with respect to time in Scenario 1.
parameter and called the deformation number (DN) that increases when the thickness increases or the elastic modulus decreases.The initial conditions of the hydraulic head and soil displacement are respectively nondimensionalized as follows: ,