# Estimation of Hydraulic Diffusivity of a Confined Limestone Aquifer at the Xiluodu Dam

Authors: Mingwei Li; Zhifang Zhou; Hailong Huang; Jianxin Liao
Journal: Geofluids
 (2022)

---

## Abstract

Hydraulic diffusivity is a fundamental parameter in studying groundwater flow. An analytical solution for groundwater flow within a finite-length one-dimensional aquifer has been introduced. This analytical solution simplifies the river water level (RWL) by a piecewise linear function, accurately describing the arbitrary variation pattern. An inversion method for hydraulic diffusivity has been provided based on the analytical solution. Then, the method was employed to estimate the hydraulic diffusivity of a limestone aquifer in the Xiluodu reservoir and was verified by the FEFLOW software. The parameter sensitivity increases with the RWL fast increasing and reaches its peak value when the RWL approaches the maximum value. Given the sensitivity analysis, great attention should be paid to the stage when RWL fluctuates drastically for estimating hydraulic diffusivity. The proposed method is flexible to conditions that lack observation data or RWL fluctuating drastically.


---

## Body

## 1. Introduction

Hydraulic diffusivity is a fundamental hydraulic parameter of an aquifer, and it is vital for the transient simulation of groundwater flow. The method of estimating hydraulic diffusivity includes the pumping test, flood wave response analysis, tidal response analysis, and numerical modelling [1, 2]. Analytical models are essential tools for estimating hydraulic diffusivity based on river-aquifer interaction [3]. Pinder et al. [4] proposed a method for calculating the hydraulic diffusivity, in which a series of discrete steps represented the stage hydrograph. Shih [5] derived an inversion model for the inversion of hydraulic diffusivity using a variation of groundwater. Shapiro et al. [6] proposed a method to estimate the hydraulic diffusivity through fluid-injection-induced seismic. The conventional approach to studying hydraulic diffusivity using river-aquifer interaction assumes that the RWL is a periodical variation [1, 7–9]. In addition to the above methods, the recently popular nature-inspired algorithms may be a good choice for estimating hydraulic diffusivity [10–13].There is often irregular fluctuation in the process of periodic variation of RWL. Many slight irregular fluctuations will be ignored when a periodic function fits the RWL. To reflect the fluctuation characteristics, it is necessary to use the arbitrary variation function to fit the RWL. Based on the assumption of an infinite aquifer, some studies have proposed the analytical solution for aquifer water level (AWL) with the arbitrary variation of RWL [14, 15]. Ramp kernels were presented by Singh [14] to accurately calculate responses within an infinite aquifer with an arbitrary RWL. The ramp kernel function requires a fixed time step, which sometimes brings difficulties to practical applications. Most methods estimate the hydraulic diffusivity with the assumption of the infinite-length aquifer [15–18]. However, there is rarely a study that estimates hydraulic diffusivity with the assumption that a finite-length aquifer with arbitrary water fluctuation.This study proposes an inversion method to estimate the hydraulic diffusivity with arbitrary fluctuation of RWL within a finite-length aquifer. We first introduce an analytical solution for the AWL response to RWL variations. Then, the solution was employed to estimate the hydraulic diffusivity of an aquifer buried beneath the Xiluodu reservoir. The accuracy of the proposed method is verified by the numerical solution of the FEFLOW software. Furthermore, the sensitivity of hydraulic diffusivity was discussed to examine the importance of observation data of different stages on parameter estimation.


## 2. Methodology

The confined aquifer extends horizontally, and the change of AWL far away from the river bank is tiny, which can be regarded as a fixed value. The conceptual diagram of the hydraulic process of the confined aquifer is shown in Figure1.Figure 1
Schematic diagram of the hydraulic response of a confined aquifer.Here are some assumptions: (1) The confined aquifer is homogeneous and isotropic, the formation extension is horizontally, and the extension length isl. (2) The boundary condition on one side of the confined aquifer is the constant water level, and the other is consistent with the RWL. (3) The hydraulic gradient of the confined aquifer is minimal under natural conditions. We use a piecewise linear function to represent the RWL. The initial boundary value problem of the hydraulic response of the confined aquifer can be expressed as
(1)∂ux,t∂t=a∂2ux,t∂x2,x>0,t>0,subject to the boundary conditions:
(2)u0,t=βit−βiti−1+∑j=1i−1βitj−tj−1,t>0,ul,t=0,t>0,and the initial condition is
(3)ux,0=0,x>0,where ux,t is the change of AWL and a is the hydraulic diffusivity, which can also be expressed as K/Ss. K is the hydraulic conductivity, and Ss is specific storage. l is the extended length of the aquifer. After the RWL is linearized, ti represents the time demarcation point, and βi represents the slope of time period i and t0=0. The solution of this problem has been given by Zhuang et al. [19, 20]:
(4)ux,t=βit−ti−1+∑j=1i−1βjtj−tj−11−xl−2l2π3a∑j=1iβj∑n=1∞1n3exp−n2π2al2t−tj∅i−j−exp−n2π2al2t−tj−1sinnπxl,where
(5)∅i−j=0j=i1j≠i.To understand the above analytical solution, we take the RWL data of the Jinsha River as an example to analyze the variation process of AWL near the river bank. A total of 54 RWL values from January 10, 1992, to December 20, 1992, were examined. First of all, the RWL was subtracted from the initial value to obtain the fluctuations (see Figure2 and Table 1). Then, it was divided into 42 periods, within which the RWL changes linearly. The RWL at the time segment point between the adjacent periods is equal. As shown in Figure 2, the linearization result can well reflect the main characteristics and ignore the small changes in RWL. Table 1 shows the segmentation point and the slope of the change of the RWL in each period.Figure 2
Linearization of reservoir water level.Table 1
Linear piecewise time node and curve slope.
No.Time nodeslopeNo.Time nodeslopeNo.Time nodeslope10–15157–0.3292330.782270–0.02127161620.40230238–0.1963800.133285171670.1431243–0.6584910.015455181720.59432248–0.03851010.221191770.571332530.3746111–0.028201870.333342580.15271160.31211920.07535263–0.03481210.23622197–0.94436268–0.4349126–0.31223202–0.33837274–0.1561110131–0.072242080.476667382840.028111360.224252130.024392890.37133312141–0.25626218–1.02840294–0.242131470.101667272230.07441304–0.38375141520.136282280.13642355–0.08513Assuming that the confined aquifer's horizontal extension length (l) is 5 km, the hydraulic diffusivity is 70000 m2/d. The hydraulic response duration is calculated at multiple points at different distances from the river bank. The response at points 0, 0.1 l, 0.2 l, 0.5 l, and 0.8 l distance from the river bank are shown in Figure 3. It can be seen that after the same time, the AWL near the river bank varies drastically and reflects the periodic variation characteristics of RWL. With the increased distance from the river bank, the variation amplitude of AWL decreases gradually.Figure 3
Hydraulic response duration curve in a confined aquifer at different points.Although the storativity of the confined aquifer is small, the large horizontal extension length gives the aquifer some water storage capacity, so the hydraulic response of each point within the aquifer has hysteresis. There is a more obvious hysteresis farther away from the bank slope. When the RWL keeps stable periodic fluctuation, the hysteresis of the hydraulic response of the confined aquifer can also be reflected by the phase difference between AWL and RWL. The AWL's phase lag increases with the distance between observed points and the bank slope.Assuming that the confined aquifer's horizontal extension length (l) is 5 km, the hydraulic diffusivity of the aquifer is 10,000, 20,000, 50,000, and 80000 m2/d, respectively. At point 0.1l, the AWL duration curve of the confined aquifer is shown in Figure 4. The hydraulic response has the following characteristics: (1) With the decreased hydraulic diffusivity (a), the variation of the AWL gradually decreases, and hysteresis characteristics become more obvious. (2) When the hydraulic diffusivity is large enough, its influence on the AWL curves is small (see curves 50000 m2/d and 80000 m2/d).Figure 4
Hydraulic response duration curve in a confined aquifer for different hydraulic diffusivity.The hydraulic diffusivity can be estimated by parameter estimation approaches such as the PEST parameter estimation package [21] and the Markov chain Monte Carlo (MCMC) algorithm [22]. All these approaches are through the continuous trial calculation to find the parameters that minimize the difference between the simulation result and the observed data. In this study, the PEST parameter estimation package was used to estimate a. PEST is a “model-independent” program, which constantly provides estimated parameters to the model and obtains the calculation results from the model. Using Equation (4), calculate the hydraulic head in a confined aquifer for any RWL hydrograph. MATLAB software is used to code Equation (4) and then combined with the PEST program. The objective function Φ of PEST, without considering the observation weight, can be expressed as
(6)Φ=∑ui∗−ui2,where ui∗ is the observation value and ui is the calculation value.The parameter means squared error (MSE) is used to estimate the accuracy of the estimation result.MSE=Φ/N, where Φ is the result of the objective function in PEST and N is the number of the observation.


## 3. Field Application

The study area of the Xiluodu hydropower station is located in Yongsheng Basin, on the Jinsha River in Yunnan Province. Before the construction of the hydropower station, many geological surveys were carried out, among which some water level observation wells were set up on both sides of the Jinsha River. Upper Permian Emeishan basalt (P2β) and Lower Permian Yangxin limestone (P1y) are the primarily exposed strata in the dam site (Figure 5). The thickness of the basalt formation in the dam area is 490~520 m, and the thickness of the limestone formation is 512 m. The confined aquifer formed by Lower Permian Yangxin limestone is the main object of this study. The limestone is exposed at the bottom of the river upstream of the dam. In the exposed part of limestone, the exchange between groundwater and surface water is robust, and the degree of karst development is high. Along the downstream direction of the river, the buried depth of limestone gradually increases, and the degree of karst development decreases slowly. Restricted by the upper and lower aquitards, the limestone confined aquifer is only recharged at the basin’s edge and discharged into the Jinsha River [23]. Near the bank slope, the AWL is mainly affected by the RWL for a short time. Four observation wells are located in the slope on both sides of the Jinsha River, X25, X35, X41, and X71. The location of these wells is shown in Figure 6, and the details are shown in Table 2.Figure 5
Generalized geology cross section of the river bank in the research area.Figure 6
Location of study site and monitoring wells.Table 2
Information of observation wells.
WellYearLocationOrifice elevation (m)Well bottom elevation (m)Aquifer thickness (m)Distance with river bank (m)Karst developmentX251992Ground402.07202.97172.6757.10Quartz hole 2 cmX351992Ground420.26119.09256.8785.70Small fractureX411995Ground409.05208.35176.84142.00FractureX711995Adit404.35-14.17399.0168.60Small solution gaps and poresThe water level data for 1992 and 1995 are used to study the hydraulic diffusivity. First, all observation data are subtracted from the initial value to get the water level changes. Then, the RWL data are linearized, as shown in Figures2 and 7. The linearization can accurately reflect the characteristics of water level changes.Figure 7
Linearization of the river water level in the year 1995.The exposed area of the limestone formation near the dam is very far away from the exposed area on the basin edge. The formation on the left and right banks near the dam site is approximately horizontally distributed. It can be assumed that the exposed area in the basin edge has no effect on the short-term groundwater level of the limestone aquifer near the dam site. Hydrological data show that the annual variation of RWL is about 20 meters. The variation of RWL is very small compared to the aquifer length. Considering the hydrogeological conditions and the amplitude of RWL, it is considered that the boundary far from the river bank has no effect on the AWL. In the estimation in this work, we assume that the aquifer length is 5000 m. The impact of aquifer length will be discussed below. The distance between the well and the river bank is set according to the real space (see Table2).The estimation result of hydraulic diffusivity is obtained using the PEST program (Table3). The comparison between the observed and calculated AWL is shown in Figure 8. The calculation results show that the range of hydraulic diffusivity is 1.09 to 6.3×104m2/d except well X41 with diffusivity equal to 2.08×106m2/d. Generally speaking, the hydraulic diffusivity in the upstream aquifer is higher than that of the downstream. The difference of hydraulic diffusivity between the left and right banks is not obvious. The main reason for the difference in diffusivity is the nonuniformity of karst development in the limestone formation (Table 2). We can expect that specific storage is constant for the whole aquifer, so higher diffusivity means higher permeability. And this nicely corresponds to Table 2–Karst development–bigger fractures, and holes means higher permeability (diffusivity). Observation well X41 is located in the upstream area, where karst is developed, and fractures are revealed in the process of drilling. Therefore, the hydraulic diffusivity of well X41 is larger than that of other wells. There are also small solution gaps in well X71, so the hydraulic diffusivity is slightly larger than that of wells X35 and X25. Wells X25 and X71 are located on the left bank and are very close to each other. Comparing these two wells, the ratio of the aquifer thickness is X71/X25=2.31, and the ratio of hydraulic diffusivity is X71/X25=5.78. The slight difference between the two ratios indicates that the aquifer thickness has little effect on hydraulic diffusivity. For limestone, the karst development degree has a great influence on the hydraulic conductivity, but has little influence on the specific storage.Table 3
Calculation results of hydraulic diffusivity.
WellHydraulic diffusivity (proposed method)Hydraulic diffusivity (FePEST)MSEValue (×104 m2/d)MSEValue (×104 m2/d)X250.471.090.521.13X350.512.530.532.81X410.66207.801.20285.46X710.816.300.5011.38Figure 8
Observed and calculated water level fluctuations in wells.
(a)(b)(c)(d)


## 4. Compare with the Numerical Solution

To verify the reliability of the analytical results, the numerical method is also applied by the FEFLOW method [24]. The settings of the 2D FEM model are as follows: the range 5000m×500m; the observed well is settled as the distance of the well shown in Table 2; the total element is 3200; total node is 3609; the simulated periods are 365 days with an initial time-step length of 0.001 d; one boundary condition was set as the observed RWL; and the other boundary condition was set as a constant value. The hydraulic parameters were set as follows: the initial water level of the total model was assumed to be 0 m; the specific storage was assumed 1×10−6m−1. After the finite element model is established by FEFLOW, the model can be imported into the inversion program (FePEST) through the PEST plug-in. The hydraulic conductivity was estimated by the PEST module with the initial value of 0.13 m/d. The hydraulic diffusivity is obtained by the ratio of hydraulic conductivity and specific storage. Figure 8 shows AWL fluctuations calculated using FEFLOW based on the estimation result of FePEST.Table3 shows the estimation results of Equation (4) and the FEFLOW method. The mean square error (MSE) values in the three wells obtained using Equation (4) are lower than those obtained using the FEFLOW method. The MSE in X71 was determined to be 0.81 and 0.50 using Equation (4) and the FEFLOW method, respectively. It indicates that Equation (4) performed better than the FEFLOW method for X71. From a comprehensive point of view of the four wells, the estimation results of the newly proposed method are consistent with those of the numerical method.


## 5. Sensitivity Analysis

To analyze the influence of various factors on the inversion results and determine the sensitivity of RWL to parameters, we conduct a sensitivity analysis. Asa is the parameter we estimate, we apply the normalized sensitivity analysis [25] to analyze the effect of a on the aquifer water level u. The normalized sensitivity of a is defined as
(7)Sa=a∂u∂a.Based on Equation (4), Sa is further derived as
(8)Sa=a∂u∂a=2l2π3a∑j=1iβj∑n=1∞1n3exp−n2π2al2t−tj∅i−j−exp−n2π2al2t−tj−1sinnπxl−2π∑j=1iβj∑n=1∞1n−t−tj∅i−jexp−n2π2al2t−tj∅i−j+t−tj−1exp−n2π2al2t−tj−1sinnπxl.Hydraulic diffusivity measures diffusion speed of pressure disturbances [26]. Figure 9 shows the sensitivity of u to avalue with time at a different place of the aquifer. In the initial stage, the RWL changes gently, and the sensitivity of the entire aquifer is minimal. When the RWL fluctuates slightly, the pressure disturbance in the aquifer is feeble, and the effect of hydraulic diffusivity on the diffusion speed cannot be reflected. In the middle stage, along with the rapid change of RWL, the sensitivity increases rapidly and reaches the peak value. With the RWL’s fast change, there is a great pressure disturbance in the aquifer, and the influence of hydraulic diffusivity is highlighted. In the final stage, with the decrease of the RWL’s change rate, the sensitivity decreases rapidly.Figure 9
Normalized sensitivity of hydraulic headu to hydraulic diffusivity a versus time.Figure9 also reveals that the sensitivity increases at first and then decreases with the distance from the river bank. The low sensitivity near the river bank is that the very short distance is not enough to reflect the difference in the diffusion velocity of water pressure represented by the hydraulic diffusivity. The increased distance from the river bank makes the difference in pressure and diffusion velocity obvious. The sensitivity reaches the peak value when the space of the river bank is about 0.2 l (1000 m). As the influence of the fluctuation of RWL weakens with the increase of distance, the sensitivity decreases when it is very far away from the river bank.If one conducts parameter estimation, it is recommended to focus on the period of RWL’s fast change, according to the sensitivity analysis. In addition, if multiple observation wells are available, an observation well near 0.2l from the bank slope seems to produce more reliable estimation parameters. If the water level is observed in multiple hydrological years, the observation data with rapid changes should be preferred.


## 6. Discussion

To verify the validation of the model proposed, we conduct a comparison with the results in the literature. For sinusoidal RWL, the analytical solution of the AWL was proposed by Ferris [9], which is written as
(9)ux,t=M0exp−xω2asinωt−xω2a,where M0 is the amplitude of RWL fluctuation, ω is radians per time unit, ω=2π/t0, t0 is duration of the flood wave, and a is the hydraulic diffusivity. The river boundary condition is u0,t=sinωt.Singh [14] developed solutions for the aquifer response in an infinite aquifer for arbitrary variation of RWL with time, which is given as
(10)ux,n∆t=∑i=1nHi−Hi−1∆tδx,n−i+1∆t,δx,m∆t=∆t1−F1m−1−m∆t+x22a∆F1m−x∆tπa∆F2m,F1m=erfx2am∆t,F2m=mexp−x24am∆t,where ∆t is the uniform time step size; i, n, and m indices for time steps; H is RWL; δ is the ramp kernel function; a is the hydraulic diffusivity; and ∆F1m=F1m−F1m−1, ∆F2m=F2m−F2m−1.To verify the accuracy of the new method, we consider a synthetic river-aquifer response test. The parameters wereM0=10m; x=50m; ω=0.04π; a=10000m2/d; ∆t=2d; and l=5000m. The AWL is compared to that calculated using the analytical solution, Singh’s solution, and the newly proposed solution (Figure 10). Figure 10 shows that in the initial stage, there is a big gap between these two solutions and the analytical solution. The reason is that both Singh’s and the newly proposed method assume that the water level is horizontal at the initial condition, which is different from Ferris’ assumption. Figure 10 shows that the performance of the newly proposed solution is equivalent to the analytical solution and as good as Singh’s solution.Figure 10
Variation of water level with time.To understand the impact of aquifer length on results, the AWL with different aquifer lengths (l) is calculated. We use MSE to characterize the error between the proposed and analytical methods (Figure 11). It can be found that when the aquifer length is short, there is a large difference between the newly proposed method and the analytical method of Ferris. When the length of the aquifer reaches a certain value (e.g., larger than 1000 m in this condition), the error tends to be stable, and the calculated results are very close to the analytical solution. We define this certain value as the influence length. The above results show that a great difference between finite-length and infinite-length aquifers when the aquifer length is less than the influence length. That is to say, when the aquifer length is less than the influence length, the nonriver boundary conditions have a noticeable effect on the AWL near the river valley. The influence length is related to the amplitude of the river water level. With the increase of the amplitude of the RWL, the influence length will be lengthened.Figure 11
Accuracy of simulation results of different aquifer length.Singh’s method assumes that the time step∆t must have the same length. However, it is difficult in some cases that the observation data are not fully obtained, such as data like 1992, Figure 2. When the RWL fluctuates greatly (Figure 7), Singh’s method needs a shorter time step to ensure the overall simulation accuracy. The newly proposed method not only can automatically determine the time step according to the characteristics of the RWL, but also it is more flexible for the condition when there is a lack of observations. Figure 11 shows that when the aquifer length is less than the influence length, the effect of the nonriver boundary on the AWL is obvious. Therefore, when the extension length of the aquifer is limited, it is more reliable to use the new method to estimate the hydraulic diffusivity.


## 7. Conclusions

This study introduced an analytical solution for groundwater flow in a finite-length aquifer with arbitrary variation boundary conditions of water level. The hydraulic diffusivity of a limestone confined aquifer was estimated by the PEST method. The normalized sensitivity of water level to the hydraulic diffusivity is then conducted. The main conclusions of this study can be drawn as follows:(1)
The arbitrary fluctuation of river water level (RWL) is linearized. The analytical solution of aquifer water level (AWL) within a finite-length aquifer with arbitrary fluctuation boundary conditions was introduced. Based on the proposed method, the hydraulic diffusivity of a confined aquifer was estimated.(2)
The parameter sensitivity increases with the magnitude of change of the RWL, and its peak value is obtained when the RWL changes fast. For estimating hydraulic diffusivity, great attention should be paid to the period of RWL’s fast change. The observation well is recommended to place at the location of 0.2l to improve the accuracy of the parameter estimation.(3)
The calculated hydraulic diffusivity is consistent with those reported by the FEFLOW method. The newly proposed model is usable for cases when there is a lack of observations and the RWL fluctuates greatly capering with the literature model Singh [14].



---

*Source: 8732415-2022-09-05.xml*