- Email: [email protected]
Depth-dependent relation between hydraulic conductivity and electrical resistivity in geologic formations
Journal of Hydrology 578 (2019) 124081
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Research papers
Depth-dependent relation between hydraulic conductivity and electrical resistivity in geologic formations ⁎
T
⁎
Kaixuan Lia, Tongchao Nana, , Xiankui Zenga, Lihe Yinb, Jichun Wua, , Jun Zhangb a b
Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China Xi’an Center of Geological Survey, China Geological Survey, China
A R T I C LE I N FO
A B S T R A C T
This manuscript was handled by Huaming Guo, Editor-in-Chief
Predictability of groundwater flow and transport models greatly relies on adequate characterization of aquifer heterogeneity such as hydraulic conductivity (K), porosity (ϕ) etc. The traditional model combining Archie formula and Kozeny-Carman equation (thus called AKC model in this study) has been proven in a large number of studies to be helpful in evaluation of K through electrical resistivity (Re). However, AKC model does not take into account the widely-reported dependence of K on depth, thus often fails in deep formations (> 100 m), as was observed in our study. In order to consider the impact of depth on the K-Re relation, a new model (called LAKC model) was introduced, which combined Louis decay model with AKC model. The effectiveness of LAKC model was verified against K measurements from stratified pumping tests and applied in a real aquifer in northern Ordos Basin, China. The results revealed that LAKC model can characterize the depth-dependent K-Re relation properly and that calibrated LAKC model is capable of providing adequate K estimates based on detailed Re measurements in a similar context of lithology.
Keywords: Aquifer characterization Electrical resistivity Northern Ordos Basin Depth dependence Stratified pumping test
1. Introduction Regional groundwater flow is critical to understand various geologic processes, such as fluid flow in sedimentary basins (Jiang et al., 2009; Wang et al., 2011), surface-ground water interaction (Winter, 1978), heat and solute transport (Cardenas, 2007; Niederau et al., 2017) etc. Aquifer heterogeneity is one of the key factors of properly modeling groundwater flow and solute transport. Unfortunately, due to technical and financial limitations, it is extremely impractical to reveal the spatial variability of real aquifers using traditional approaches like pumping tests and drilling (Binley et al., 2015; Nan and Wu, 2017; Yeh and Liu, 2000). It is still challenging for hydrogeologists to collect sufficient and reliable information to infer aquifer heterogeneity such as detailed distributions of hydraulic parameters (e.g. hydraulic conductivity, K and porosity, ϕ). The predictability and reliability of groundwater models (especially those of solute migration) can thus be significantly undermined in highly heterogeneous media (Mizukami et al., 2017; Nan et al., 2019; Yoon et al., 2013; Zheng et al., 2011). In recent decades the emergence of and advances in hydrogeophysics have made it possible to collect information of aquifer heterogeneity with growing amounts and accuracy, and even to estimate hydraulic parameters quantitatively (Dhakate and Singh, 2005; Frohlich, 1994; Niwas and de Lima, 2003; Rezaei et al., 2016; Zakari ⁎
et al., 2015). The geophysical approaches provide fast and effective technical supports for aquifer exploration and evaluation, and have become alternatives to traditional hydrogeological approaches to obtain hydrogeological parameters (Binley et al., 2015; Singh, 2005). The similarity between Ohm's law and Darcy’s law as well as linkage between geoelectric and hydraulic properties of formations made it possible to estimate hydrogeological parameters through geoelectric methods (Louis et al., 2004). There have been many studies on the relationship between electrical resistivity (Re) and hydrogeological parameters in sandy media, among which the Archie formula had been accepted by many researchers (Archie, 1942; Dick et al., 2018; Robinson et al., 2012; Singha and Gorelick, 2006; Whitman and Yeboah-Forson, 2015; Yeh et al., 2002; Zha et al., 2017). In a large number of studies Re obtained by geophysical methods were converted into hydraulic conductivity (K) through a model combining Archie equation and Kozeny-Carman equation (termed “AKC model” for brevity), and the performance of AKC model was acceptable in shallow formations (Kazakis et al., 2016; Niwas and Celik, 2012; Sikandar and Christen, 2012). Though in near-surface environments (< 1km) heterogeneity (e.g. lithofacies) is a controlling factor of permeability (or equivalently, K) (Gleeson et al., 2011; Ranjram et al., 2014), permeability and K were also found dependent on locale depth. In fact, a number of studies
Corresponding authors. E-mail addresses: [email protected] (T. Nan), [email protected] (J. Wu).
https://doi.org/10.1016/j.jhydrol.2019.124081 Received 13 May 2019; Received in revised form 23 August 2019; Accepted 26 August 2019 Available online 27 August 2019 0022-1694/ © 2019 Elsevier B.V. All rights reserved.
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
fine sandstone and pelite. Since no continuous aquifuge layer exists above J2a, the formations within 800 m deep can be treated as a unified aquifer of large thickness (over 800 m) (Yin et al., 2011). The total depths of boreholes ZKH2, ZKH3 and ZKH5 are 958 m, 703 m and 753 m, respectively. Three boreholes expose Quaternary Holocene aeolian sand layer (Q4eol) at depth 0 m, 0–10 m and 0–1.6 m, respectively. Huanhe formation (K1h) is found in the three boreholes at 0–474 m, 10–330 m and 1.6–574 m, respectively. Luohe formation (K1l) are 474–880 m, 330–701 m and 574–753 m, respectively. Anding formation (J2a) is exposed by ZKH2 and ZKH3 at 880 m and 701 m (see Fig. 2). In this study, XSKC-92 logging instrument was used for well logging. By using this instrument, the borehole temperature and measured resistivity (R0) values of depth 5–700 m were collected from ZKH5, ZKH2 and ZKH3 boreholes, with electrode space equal to 5 cm. The range of temperature in three boreholes at all depths are 19–30 °C, and the average is about 24 °C. Considering the influence of ground temperature variation on R0 within 700 m depth, the temperature correction for R0 is carried out by Arps' formula (Arps, 1953):
demonstrated that permeability generally decays with depth (Jiang et al., 2009; Louis, 1974; Manning and Ingebritsen, 1999; Ranjram et al., 2014; Snow, 1968; Wang et al., 2009). Such decay behaviors can be remarkable in thick aquifers and significantly change the conditions of regional groundwater flow systems (Jiang et al., 2010; Saar and Magna, 2004; Zhao, 1998). But quantitative descriptions of depth-related decay is still limited (Sanford, 2017). There were studies on depth’s impact on permeability through empirical or semi-empirical methods (see Jiang et al., 2011; Kuang and Jiao, 2014; and references therein). To describe the general relationship between permeability and depth a few models were proposed, among which two major types were widely used, i.e. the power-law model (Bredehoeft et al., 1992; Ingebritsen and Manning, 2002, 2010; Manning and Ingebritsen, 1999; Stober and Bucher, 2007) and the exponential model (Jiang et al., 2010; Louis, 1974; Saar and Magna, 2004; Williams and Narasimhan, 1989). These models were aimed to a general rule of depth-permeability relation so that intrinsic heterogeneity of medium (such as lithofacies) was ignored and variations in observed permeability data in a large range of depth were fully attributed to depth (Kuang and Jiao, 2014). In thick and highly heterogeneous aquifers, it is better to take into account the impacts of both formation heterogeneity and depth. In this paper, we develop a depth-dependent K-Re model which combines a depth-permeability model of Louis (1974) exponential type and widely used AKC model (hence the proposed model is called LAKC model). Rather than analyzing the impact of gravity stress (depth) on each AKC parameter or intermediate separately, a lumped exponential function is used to fully account for gravity impact to the final K estimates. Our goal is to consider the impacts of aquifer heterogeneity and depth simultaneously so that K in thick regional aquifers can be estimated by Re measurements collected by geophysical methods in a parsimonious manner. To analyze the applicability of LAKC model, we investigate the performance of LAKC model in a cross-validation experiment based on K and Re measurements collected in three boreholes of depth > 700 m in a thick aquifer. That is, LAKC model is calibrated via PEST code using measurements from any two boreholes and validated using measurements from the third borehole. To demonstrate the necessity of including depth dependence, AKC model is also applied in a similar procedure and compared to LAKC model. After validation of LAKC model, we apply it in estimating detailed K values from 5 cm-spacing Re data throughout the boreholes. This paper is structured as follows. Section 2 explains AKC model, LAKC model and calibration procedures. Section 3 presents the information of the study area and available data. In Section 4, we analyze the predictability of LAKC model in a cross-validation experiment compared to AKC model and apply LAKC model throughout depth to obtain detailed K distributions and discuss depth-dependency of effective porosity. Section 5 provides a summary and the main conclusions.
Re =
T0 + 21.5 R0 Te + 21.5
(1)
where T0 (°C) is the borehole temperature, R0 (Ω·m) is the measured resistivity and Re (Ω·m) is the electrical resistivity at a base-line temperature Te (equal to average temperature in three boreholes). Stratified pumping tests based on Solexperts Heavy-Duty Double Packer System (HDDP) were conducted in 44 borehole segments of length 20 m or 50 m. HDDP is a set of drilling water sample acquisition system consisting of two balloon embolisms (packers). In this method, two packers were inflated with high-pressure nitrogen, and in close contact with the well wall to produce water-stopping effect. The nontarget sections at both ends of the target zone (segment) were isolated through Packer air bag, and then the target zone was pumped by submersible pump. Meanwhile, the groundwater level of the target zone (segment) was monitored and recorded in real time by a probe. Groundwater flows were approximated as radial flow in horizontal direction in data interpretation. The interpreted (measured) 44 K values represent the mean horizontal K on 44 segments and were shown in Fig. 3. The interpreted K varies over three orders of magnitude, showing high heterogeneity as well as a plausible decay trend with depth. Total dissolved solids (TDS) was measured in laboratory based on water samples collected in stratified pumping tests. Rw profiles (at average temperature 24 °C) were calculated based on interpolated TDS measurements (e.g. Fig. 5(d) below) and reciprocal relation of TDS-Rw. The relevant data can be found in Nan (2019). 3. Materials and methodology 3.1. AKC model: a classic K-Re relation based on Archie and KozenyCarman formulae
2. Study area and available data Archie formula reflects the relationship between (electrical) resistivity of granular medium and medium porosity, which reads as:
The hydrogeological exploration boreholes ZKH5, ZKH2 and ZKH3 in this study are located in the middle and low desert plain of Wushenqi area, the lake-concentrated area in northern Ordos Basin, China (Fig. 1). The strata in this area are almost horizontal, and all the strata are terrigenous clastic deposits. The hydrogeological profile of the study area is shown in Fig. 2. The aquifer under study mainly consists of Lower Cretaceous Luohe formation (K1l) and Lower Cretaceous Huanhe formation (K1h). The lithofacies are semi-cemented siltstone, fine-/medium-/coarse-sandstone, gritstone and sandwiched silty mudstone, with well-developed pores and low shale content. At the top is discontinuously-distributed Quaternary Holocene aeolian sand layer (Q4eol, thickness < 10 m). The underlying formation below K1l is Jurassic Anding formation (J2a, depth > 800 m in most areas), consisting of low-permeability dense
R e = αRw ϕ−mS
−n
(2)
where Re is volume resistivity of rock matrix; Rw is resistivity of fluid (water in aquifers); ϕ is the porosity of the medium; n is saturation index; m is cementation index, which is related to lithology, pore structure and diagenesis; α is the value related to the aqueous medium; S is the saturation of the medium, and for the saturated aquifer, S = 1. In general, the values of α and m are obtained by testing core samples in laboratory and then applied to formations of similar lithology (Kazakis et al., 2016). If there is a lack of drilling samples in the study area, empirical values can be given according to lithofacies and documented ranges in literature (Aristodemou and Thomas-Betts, 2000). As mentioned above, Re and Rw data in this study are all corrected to base-line 2
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 1. Sketch map of the study area and borehole locations.
thus ignored due to small temperature variations (< 11 °C) in three boreholes; g is the gravitational acceleration 9.81(m/s2); de is the particle size (mm), set to be average particle size of the corresponding lithological section according to the bore histogram of three boreholes. The range of de used in this study is listed in Table 1. The classical AKC model is obtained by combining Eqs. (2) and (3) to eliminate porosity ϕ, that is,
temperature 24 °C. Note that in formations of low clay content, temperature has similar impacts on Re and Rw and hence the ratio Re/Rw can be considered temperature-independent (Waxman and Thomas, 1974). Though temperature correction is not really required to apply Eq. (2) in this area, these corrections are still carried out in both Re and Rw data, simply in order to filter out temperature-induced trend in data visualization (e.g. Re and TDS profiles in Fig. 5 below). In porous media ϕ is closely related to K. The pore characteristics (such as number, size, shape, direction, degree of connectivity and spatial changes) determines the soil storage, retention, release and transmission performance of water (Niu et al., 2015). For two similar sandy media, one with a higher ϕ usually has a higher K (Sikandar and Christen, 2012). Kozeny-Carman formula can be expressed as: (Carman, 1937; Kozeny, 1927)
KK−C =
ρw gde2 ϕ3 180μ w (1 − ϕ)2
KAKC =
ρw gde2 (αR e / Rw )1/ m · 180μ w [(Rw / αR e )1/ m − 1]2
(4)
KAKC represents K value calculated by AKC model. Traditionally, AKC model is used to convert rock resistivity collected by geophysical method into K. This method of estimating K by ϕ has been found effective in granular media with low clay content (Flint and Selker, 2003; Urumović and Sr, 2016).
(3)
3.2. LAKC model: a new K-Re relation based on Louis negative exponential model
where KK-C is calculated hydraulic conductivity (m/s); ρw is water density, equal to 1000 kg/m3; μw is dynamic viscosity of water, equal to 0.0014 kg/m·s in this study (Fetter, 1994). Temperature dependence of water density and viscosity are very small (Al-Shemmeri, 2012) and
Although AKC model is useful in estimating K in shallow aquifers (Kazakis et al., 2016; Niwas and Celik, 2012; Sikandar and Christen,
Fig. 2. Hydrogeological profile of the study area. Red sections indicate the tested segments in stratified pumping tests. Dash lines show formation division. 3
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 3. Depth cover and interpreted K values at 44 test segments in stratified pumping test at ZKH5, ZKH2 and ZKH3.
K decays with depth. Researchers tried to find one or several universal attenuation rate(s) for Earth’s crust (Kuang and Jiao, 2014; and references therein), but strictly speaking, A actually depends on compressibility of the medium and should have different values for different lithofacies. Fortunately, the variability of rock compressibility is much weaker than K (Domenico and Mifflin, 1965) so that one can approximately use a single A value within a stratum without introducing remarkable errors. In this study, two individual values of A are separately calibrated for Huanhe formation (K1h) and Luohe formation (K1l).
2012), it does not consider K’s decay with depth and generally fails in thick aquifers where differences in self-weight stress are significant. Here we employ a negative exponential model which was proposed based on statistical data by Louis (1974) and then widely accepted in applications (Anderman and Hill, 2003; Zhang and Franklin, 1993). Louis’ negative exponential model is as follows:
K (h) ≈ K 0·exp(−Ah)
(5)
where K0 is K at surface; K(h) is K at depth h; A is an attenuation coefficient. The larger A is, the faster K decreases with depth. Since AKC model works properly in near-surface sandy aquifers, AKC-based K estimate is taken as K0 in Eq. (5) here, i.e., K0 = KAKC. By combining Eqs. (4) and (5), one immediately obtains a depth-dependent model for K- Re relation i.e. Eq. (6) or (7), which we refer to as LAKC model:
K (h) ≈ KAKC·exp(−Ah)
K (h) =
3.3. Model calibration, cross-validation and error evaluation To investigate the performance of AKC and LAKC models, we calibrate the models against K and Re data from two boreholes and validate the models using data from the third borehole. Cross-validations are carried out for all three distinct calibration-validation combinations of boreholes, i.e. C1: (ZKH2 + ZKH5 | ZKH3), C2: (ZKH2 + ZKH3|ZKH5) and C3: (ZKH3 + ZKH5|ZKH2). A parameter-estimation software, PEST (Doherty, 2008), is used for automatic model calibration. PEST is independent of calculation models and has been widely used in groundwater and soil water transport models (see e.g. Fang et al., 2010;
(6)
ρw gde2 (αR e / Rw )1/ m · ·exp(−Ah) 180μ w [(Rw / αR e )1/ m − 1]2
(7)
While AKC has two parameters α and m, LAKC has one more parameter A to be determined. As mentioned above, A reflects how fast Table 1 The range of particle size (de) used in this study. Lithology
Gritstone
coarse sandstone
Medium sandstone
Fine sandstone
Siltstone
Clay
Particle diameter (de) (mm)
>1
0.5–1
0.25–0.5
0.05–0.25
0.002–0.05
< 0.002
4
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
parameter ranges does not improve the fitting much (not shown). The decay of K with depth has to be considered in these boreholes for better fitting.
Zyvoloski et al., 2003). The ranges and initial guesses of parameters are given according to empirical values in literature (Cai et al., 2017). That is, the ranges of parameters α and m are set to 0.48–1.5 and 1.1–3.0, and initial values are 1.0 and 1.2, respectively, in both models. The initial value and range of A is set to 0.016 m−1 and 0.003–0.03 m−1. Note that A has two independent values in Huanhe (K1h) and Luohe formations (K1l). The root mean square error (RMSE) is employed to evaluate the model performances. RMSE can be calculated by the following Eq. (8):
RMSE =
1 n
4.2. LAKC model Fig. 6 reports the cross-validation results of LAKC model. It can be seen that the estimation results seem reasonable and close to each other for different schemes of data usage. The estimation result suggests that calibrated LAKC model is consistent with all the data. With different value of α and m for three data usage, the RMSE for all data in C1, C2 and C3 data usage are 0.073, 0.124 and 0.079, respectively. Fig. 7 shows profiles of fitted and measured K values in three boreholes for C1 with parameter estimates α = 0.56, m = 1.60, A1 = 0.0155 m−1 and A2 = 0.0070 m−1. The agreement between measured and calculated K profiles confirms the validity of calibrated LAKC model. Similar agreements can be found for the other two groups of estimates thus omitted for brevity. To directly compare the goodness of fits for AKC and LAKC models, scatter plots of measured and calculated K values in two models for data usage C1 are shown in Fig. 8. Coefficient of determination R2 and RMSEs are also included. Fig. 8(a) and 8(b) report the results in calibration and validation stages, respectively. In calibration, RMSEs for AKC and LACK models are 1.410 and 0.068 respectively, and R2 relative to y = x line are −1.038 and 0.944. In validation, RMSEs for AKC and LAKC are 1.030 and 0.087, respectively and the corresponding R2 values are 0.902 and 0.999. Both high RMSE and low R2 indicate that AKC model is not proper to fit the data. A negative R2 in AKC model means that the trend of line y = x is inconsistent with trend in data. In contrast, very low RMSEs and high R2 show that fitted K values by LAKC model are close to measured K values and that LAKC model performs far better than AKC model in both calibration and validation. AIC and BIC values of AKC and LAKC models (number of parameters, p = 2 and 4, respectively) are calculated using Eqs. (9) and (10). AIC of AKC and LAKC are 28.2 and −221.8, respectively; BIC values of them are 31.8 and −214.6, respectively. Clearly, LAKC model has smaller values than AKC model. LAKC model performs much better than AKC model, not because of parameter number. The good performance of LAKC may result from two major reasons, i.e. (1) the ideal geological environment at the site and (2) the significance of gravity impact. The formations mainly consist of sandstones, siltstone and gritstone, and have very low content of clay. Thus, it provides a good situation in which Archie and Kozeny-Carman formulae properly apply. In fact, AKC model works fine in fitting measurements in shallow segments (depth < 100 m). When the entire profile is considered, gravity impact on the profile is primary and overwhelming among those factors which AKC doesn’t take into account. As discussed above, the attenuation coefficient A in LAKC model reflects the variation degree of K and relates to compressibility of media. In this study, the attenuation coefficient A1 in Huanhe formation is larger than A2 in Luohe formation. It implies the decrease of A with increase in depth and compaction degree of strata, which is consistent
n
∑ (Km (i) − K c (i))2 i=1
(8)
where Km is the measured K at the i-th test segment and Kc is the arithmetic mean of calculated K within the i-th test segment according to AKC or LAKC model, and n is the total number of test segments. RMSE reflects the overall deviation of the model from reference (data). The smaller RMSE is, the better the model fits the data. The Akaike information criterion(AIC) (Akaike, 1973) and Bayesian information criterion (BIC) (Schwarz, 1978) were also used to contrast the effects of LAKC model and AKC model in parameter estimation. The AIC and BIC can not only measure the deviation between measured data and model predicted data, but also reflect the instability caused by the different number of model parameters (Enemark et al., 2019). The smaller value of AIC and BIC is, the better the parameter estimation model performs. They are calculated as follows:
AIC = nln(ESS/n )+ 2p BIC = nln(ESS/ n) + p ln n
(9) (10)
ESS is the residual sum of squares between Km and the Kc calculated by AKC and LAKC model, n is the total number of test segments, and p is the number of model parameters. 4. Results and discussions 4.1. AKC model Fig. 4 reports calibration results in cross-validation experiment. In all calibrations, both parameters α and m of AKC model consistently stop at lower bounds of given ranges (i.e. α = 0.48, m = 1.1). Since estimated α and m are the same for C1, C2 and C3, the calibrated AKC models are identical and the resultant RMSEs are all equal to 1.317. It implies that no optimal fitting is found for AKC model within given parameter ranges. Fig. 5 reports fitted K profiles based on AKC model and estimated parameters (α = 0.48, m = 1.1) compared to measured K values. Re profiles are also included for comparison. The fitted K values tend to be much (1–2 orders of magnitude) larger than measured K values, especially at deep segments. The remarkable differences between measured and fitted K values confirm that AKC model fails in fitting the data for given parameter ranges without depth attenuation effect. Note that the given ranges of two parameters include all documented possible values for AKC model so far. Even further extension of
Fig. 4. The cross-validation results of AKC model. Final estimates are shown by solid lines inside boxes while top and bottom of boxes show 95% confidence intervals (CIs). Red line indicates initial values (IV); dash lines are upper and lower bounds (UB and LB) of parameter ranges. 5
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 5. (a–c) K values calculated by AKC model compared to measured K values in three boreholes for data usage C1. Profiles of resistivity Re are also included for comparison. (d) TDS data collected in ZKH3.
Fig. 6. The cross-validation results of LAKC model. Final estimates are shown by solid lines inside boxes while top and bottom of boxes show 95% confidence intervals (CIs). Red line indicates initial values (IV); dash lines are upper and lower bounds (UB and LB) of parameter ranges.
with the finding of Jiang et al. (2011). Parameters α and m are usually determined through core electricity experiments. In this study, the parameters were obtained in parameter estimation and then verified against the empirical values in similar lithotypes in literature. The results show that the estimated values of α and m in the study are close to those in similar lithotypes in literature (Carothers, 1968). In addition, although the Cretaceous Huanhe
Formation and Luohe Formation are two different sedimentary strata, they are lithologically similar. Therefore, the same set of parameters are used in both Huanhe Formation and Luohe Formation in this paper. 4.3. K profiles of 5-cm spacing based on Re data Based on Re data of 5-cm spacing in well log of three boreholes, 6
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 7. K values calculated by LAKC model compared to measured K values in three boreholes for data usage C1.
Fig. 8. Calculated K versus measured K values for data usage C1.
siltstone (e.g. at 103–134 m and 484–524 m in ZKH5, 83–87 m and 468–474 m in ZKH2, 306–313 m in ZKH3) or a thin layer of clay of thickness 0.5 m or less. Besides, the same type of rocks at different depths tend to have different K. For example, the averaged K of medium sandstone in ZKH2 at depth of 90–112 m, 333–376 m and 448–468 m are 0.327, 0.077 and 0.022 m/d, respectively. It suggests once again that apart from lithology, the influence of depth on K must be considered at large depth. A comparison of measured and the log K calculated by previous permeability-depth models, such as the power law model (Manning and Ingebritsen, 1999), exponential model (Saar and Magna, 2004) and Kuang-Jiao model (Kuang and Jiao, 2014), are made here (Fig. 10). The formula K = ρw gk/μw was used to convert permeability k (m2) to hydraulic conductivity K (m/d), and the ρw, g, μw take the same value as
estimated K profiles of 5-cm spacing from AKC model and LAKC models are shown in Fig. 9. For AKC model, α = 0.48, m = 1.1 according to parameter estimation above; for LAKC model, α = 0.56, m = 1.60, A1 = 0.0155 m−1 and A2 = 0.0070 m−1. It is found that despite of scale difference, K profiles from LAKC model (blue line) in three boreholes are all compatible with measured K values at segments (red line) in trend but reveal much more details of heterogeneity. On the other hand, K profiles from AKC model (black line) and measured K profiles are quite divergent, especially the trend in ZKH5 and ZKH2. According to Fig. 9 and lithology records in boreholes, strata with relatively high K generally corresponds to coarse sandstone and medium coarse sandstone, such as at (385–485 m) in ZKH5, (113 m–320 m) in ZKH2, and (120–252 m), (362–513 m) in ZKH3. And depth with low K often corresponds to fine sandstone, argillaceous 7
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 9. Hydraulic conductivity curve of borehole ZKH5/2/3 for data usage C1. Red line: K measured in pumping test; blue line: K calculated by LAKC model; black line: K calculated by AKC model.
contribution of pore water to the electrical resistivity, while KozenyCarman formula invokes an effective porosity ϕe for seepage. ϕe and ϕ are close to each other in shallow sandy aquifers and can be eliminated as intermediate variables in AKC model. But in aquifers of large depth, ϕe and ϕ may be remarkably different and one doesn’t have AKC model (Eq. (4)). If we assume that Kozeny-Carman formula is still valid, and denote
the parameters in AKC and LAKC model. As shown in Fig. 10, these existing permeability-depth models can reflect the attenuation characteristics of K with depth, but the accuracy is far from enough for the study of aquifer heterogeneity. Furthermore, it is difficult to describe the outlier values of K (data points in red circles in Fig. 10), which actually impact groundwater migration significantly. It can be seen that general permeability-depth models are not good enough in aquifer characterization and that local heterogeneity has to be included. The RMSE and R2 for different models are also shown in Fig. 10. The method of calculating average log K in different models is the same as that used in subsections above. The RMSE of exponential model is the lowest in three permeability-depth models, and the RMSE of power law model is the highest because of its spurious large values near the surface, which is an inherent defect of the model (Kuang and Jiao, 2014). The exponential and Kuang-Jiao models seem better at fitting the data except four “outliers” in red circles (ZKH3), yielding moderate R2 values. AKC model does not consider depth and has a result as bad as those of depth-only models or worse. Negative values of R2 which appear in fitting of AKC model and power law model, reflects that the trend of log K calculated by these two models is not consistent with that of measured K. LAKC model has the lowest RMSE and highest R2 among all the models for all data combinations. It shows that LAKC model considering both depth and lithology has higher accuracy than those models considering only one single factor.
the effective porosity as ϕe, then K (h) = LAKC model, one gets:
K (h) =
ρ w gde2
·
ϕe3
180μ w (1 − ϕe)2
ϕe3 ρw gde2 ρ gde2 ϕ3 · · exp(−Ah) = w 180μ w (1 − ϕe)2 180μ w (1 − ϕ)2
. Upon adopting
(11)
Eq. (11) allows one to calculate ϕe and ϕ. Fig. 11 reports the results of ϕe and ϕe/ϕ with depth in ZKH2. The jumps in both subplots at formation interface (depth 475 m) are related to our selection of independent A parameter for Luohe and Huanhe formations. Fig. 11 shows that despite discontinuity at formation interface, both ϕe and ϕe/ ϕ have a decreasing trend. ϕe/ϕ becomes far below 50% for large depth, which is consistent with our analysis of the possible reason of AKC failure above. 5. Conclusions Aquifer heterogeneity has significant influences on groundwater flow field, solute and heat transport. The traditional AKC model which relates K to resistivity Re does not consider the impact of depth on K-Re relation. According to stratified pumping tests at 44 segments in three boreholes, K of shallow strata (< 100 m) is much (1–2 orders of magnitude) larger than that in deep strata. Based on the traditional AKC
4.4. Depth dependency of effective porosity ϕe and ratio ϕe/ϕ A possible explanation of AKC model’s failure is the large difference between total porosity ϕ and effective porosity ϕe in thick aquifers. Archie formula requires a total porosity ϕ to fully account the 8
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 10. Comparison of measured log K and log K calculated by models proposed by predecessors (In these models, k is the permeability (m2) and z is the depth (km)). (“All data” represents the all 44 K values measured in ZKH5, ZKH2, ZKH3.)
Fig. 11. The effective porosity ϕe and ϕe/ϕ in ZKH2. 9
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
model and the newly established LAKC model, K within an aquifer of 700 m deep was estimated using Re data from well logs. In the attempt of fitting data from stratified pumping tests using AKC model, it was found that AKC model was unable to fit the data to an acceptable extent, especially for the decay trend in K. Though depth probably also impacts rock resistivity, standard AKC model seems inappropriate in thick aquifers and a depth-based modification to AKC model is necessary. As a depth-modified version of AKC model, LAKC model was able to match all the data very well and presented the decay trend in K with increasing depth. Results from cross-validation suggested that LAKC model calibrated by data from any two boreholes provides reliable K estimates for the third boreholes. According to lithology record and pumping tests in three boreholes, K values of similar lithofacies at different depths are quite different, showing a decay trend with depth. In order to infer K in thick aquifers based on borehole resistivity, one should take into account not only lithology-induced resistivity variation but also depth-led K decay simultaneously. This study demonstrated not only the potential of LAKC model in K estimation based on resistivity, but also the necessity to consider K decay with depth in quantitative description of heterogeneity in thick aquifers. It is also possible to reflect depth impact in AKC parameters (e.g. cementation index m) instead of introducing Louis decay function, which will be explored in our following study. In the future, more geophysical methods (such as the magnetotelluric, transient electromagnetic method, etc.) can be combined to transform resistivity data to K, which will be a promising way to obtain the spatial distribution of K with a larger scale and higher dimension.
Cardenas, M.B., 2007. Potential contribution of topography-driven regional groundwater flow to fractal stream chemistry: residence time distribution analysis of Tóth flow. Geophys. Res. Lett. 34. https://doi.org/10.1029/2006GL029126. Carman, P.C., 1937. Fluid flow through granular beds. Chem. Eng. Res. Des. 75, S32–S48. https://doi.org/10.1016/S0263-8762(97)80003-2. Carothers, J.E., 1968. A statistical study of the formation factor relation. Log Anal. 9, 13–20. Dick, J., Tetzlaff, D., Bradford, J., Soulsby, C., 2018. Using repeat electrical resistivity surveys to assess heterogeneity in soil moisture dynamics under contrasting vegetation types. J. Hydrol. 559, 684–697. Dhakate, R., Singh, V.S., 2005. Estimation of hydraulic parameters from surface geophysical methods, Kaliapani Ultramafic Complex, Orissa, India. J. Environ. Hydrol. 13, 1–11. Doherty, J., 2008. PEST, model independent parameter estimation user manual. Watermark Computing, Corinda, Australia. Domenico, P.A., Mifflin, M.D., 1965. Water from low-permeability sediments and land subsidence. Water Resour. Res. 1, 563–576. Enemark, T., Peeters, L.J.M., Mallants, D., Batelaan, O., 2019. Hydrogeological conceptual model building and testing: a review. J. Hydrol. 569. https://doi.org/10. 1016/j.jhydrol.2018.12.007. Fang, Q.X., Green, T.R., Ma, L., Erskine, R.H., Malone, R.W., Ahuja, L.R., 2010. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions. Soil Sci. Soc. Am. J. 74, 1897–1913. https://doi.org/10.2136/sssaj2009.0380. Fetter, C.W., 1994. Applied Hydrogeology, third ed. Prentice-Hall Inc., New Jersey. Flint, L.E., Selker, J.S., 2003. Use of porosity to estimate hydraulic properties of volcanic tuffs. Adv. Water Resour. 26, 561–571. Frohlich, R.K., 1994. The electric-hydraulic relationship. A geophysical model. Trends Hydrogeol. 1, 347–358. Gleeson, T., Smith, L., Moosdorf, N., Hartmann, J., Dürr, H.H., Manning, A.H., Van Beek, L.P.H., Jellinek, A.M., 2011. Mapping permeability over the surface of the Earth. Geophys. Res. Lett. 38, 93–104. Ingebritsen, S.E., Manning, C.E., 2010. Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism. Geofluids 10, 193–205. Ingebritsen, S.E., Manning, C.E., 2002. Diffuse fluid flux through orogenic belts: implications for the world ocean. PNAS 99, 9113–9116. Jiang, X.W., Wan, L., Wang, X.S., Ge, S.M., Liu, J., 2009. Effect of exponential decay in hydraulic conductivity with depth on regional groundwater flow. Geophys. Res. Lett. 36, 88–113. Jiang, X.W., Wang, X.S., Li, W., Ge, S., 2011. An analytical study on stagnant points in nested flow systems in basins with depth-decaying hydraulic conductivity. Water Resour. Res. 47, 128–139. Jiang, X.W., Wang, X.S., Wan, L., 2010. Semi-empirical equations for the systematic decrease in permeability with depth in porous and fractured media. Hydrogeol. J. 18, 839–850. Kazakis, N., Vargemezis, G., Voudouris, K.S., 2016. Estimation of hydraulic parameters in a complex porous aquifer system using geoelectrical methods. Sci. Total Environ. 550, 742–750. Kozeny, J., 1927. Uber Kapillare Leitung des Wasser im Boden, Wien Akad. Wiss. Sitzungsber. Akad.Wiss. Wien, Math. Naturw Klasse, Abt.II A. Kuang, X., Jiao, J.J., 2014. An integrated permeability-depth model for Earth“s crust. Geophys. Res. Lett. 41, 7539–7545. Louis, C., 1974. Rock Hydraulics BT – Rock Mechanics. Springer Vienna, Vienna, pp. 299–387. https://doi.org/10.1007/978-3-7091-4109-0_16. Louis, I.F., Karantonis, G.A., Voulgaris, N., Louis, F., 2004. The contribution of geophysical methods in the determination of aquifer parameters: the case of Mornos River delta, Greece. Res. J. Chem. Environ 8, 41–49. Manning, C.E., Ingebritsen, S.E., 1999. Permeability of the continental crust: implications of thermal data and metamorphic systems. Rev. Geophys. 37, 127–150. Mizukami, N., Clark, M., Newman, A.J., Wood, A.W., Samaniego, L., 2017. Toward seamless large domain parameter estimation for hydrologic models. Water Resour. Res. 53. Nan, T., 2019, Hydraulic conductivity and resistivity measurements in boreholes in Ordos, China, OSF, https://doi.org/10.17605/OSF.IO/78CQW. Nan, T., Wu, J., 2017. Application of ensemble H-infinity filter in aquifer characterization and comparison to ensemble Kalman filter. Water Sci. Eng. 10 (1), 25–35. https://doi. org/10.1016/j.wse.2017.03.009. Nan, T., Xue, L., Wu, J., 2019. Efficient identification of preferential flow path in heterogeneous media based on stream function. J. Hydrol. 577, 123961. Niederau, J., Ebigbo, A., Marquart, G., Arnold, J., Clauser, C., 2017. On the impact of spatially heterogenous permeability on free convection in the Perth Basin, Australia. Geothermics 66, 119–133. Niu, Q., Fratta, D., Wang, Y.H., 2015. The use of electrical conductivity measurements in the prediction of hydraulic conductivity of unsaturated soils. J. Hydrol. 522, 475–487. Niwas, S., Celik, M., 2012. Equation estimation of porosity and hydraulic conductivity of Ruhrtal aquifer in Germany using near surface geophysics. J. Appl. Geophys. 84, 77–85. Niwas, S., de Lima, O.A., 2003. Aquifer parameter estimation from surface resistivity data. Groundwater 41, 94–99. Ranjram, M., Gleeson, T., Luijendijk, E., 2014. Is the permeability of crystalline rock in the shallow crust related to depth, lithology or tectonic setting? Geofluids 15, 106–119. Rezaei, M., Saey, T., Seuntjens, P., Joris, I., Boënne, W., Van Meirvenne, M., Cornelis, W., 2016. Predicting saturated hydraulic conductivity in a sandy grassland using proximally sensed apparent electrical conductivity. J. Appl. Geophys. 126, 35–41. Robinson, J., Slater, L., Schaefer, L., 2012. Evidence for spatial variability in hydraulic
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 41730856, 41602250, U1503282) and the project “Hydrogeological investigation at 1:50 000 scale in the lakeconcentrated areas of the northern Ordos Basin” of China Geological Survey (Grant No. DD20160293). Data in all figures will be available upon publication of this paper. We declare that we have no conflict of interest. References Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Second International Symposium on Information Theory. Budapest, pp. 261–304. Al-Shemmeri, T., 2012. Engineering Fluid Mechanics. Ventus Publishing ApS, Holland. Anderman, B.E.R., Hill, M.C., 2003. MODFLOW-2000, the U.S. Geological Survey modular ground-water model-three additions to the hydrogeologic-unit flow (HUF) package: alternative storage for the uppermost active cells (SYTP parameter type), flows in hydrogeologic units, and the hydraulic-con. U.S. Geological Survey OpenFile Report 03-347. Archie, G., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 146, 54–62. Aristodemou, E., Thomas-Betts, A., 2000. DC resistivity and induced polarisation investigations at a waste disposal site and its environments. J. Appl. Geophys. 44, 275–302. Arps, J.J., 1953. The effect of temperature on the density and electrical resistivity of sodium chloride solutions. J. Petrol. Technol. 5, 17–20. https://doi.org/10.2118/ 953327-G. Binley, A., Hubbard, S.S., Huisman, J.A., Revil, A., Robinson, D.A., Singha, K., Slater, L.D., 2015. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res. 51, 3837–3866. Bredehoeft, J.D., Belitz, K., Sharp-Hansen, S., 1992. The hydrodynamics of the Big Horn Basin: a study of the role of faults. Am. Assoc. Pet. Geol. Bull. 76, 530–546. Cai, J., Wei, W., Hu, X., Wood, D.A., 2017. Electrical conductivity models in saturated porous media: a review. Earth Sci. Rev. 171, 419–433. https://doi.org/10.1016/j. earscirev.2017.06.013.
10
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
San Andreas fault: a testing of hypotheses. Earth Planet. Sci. Lett. 92, 131–143. Winter, T.C., 1978. Numerical simulation of steady state three-dimensional ground water flow near lakes. Water Resour. Res. 14, 245–254. Whitman, D., Yeboah-Forson, A., 2015. Electrical resistivity and porosity structure of the upper Biscayne Aquifer in Miami-Dade County, Florida. J. Hydrol. 531, 781–791. Yeh, T.-C.J., Liu, S., 2000. Hydraulic tomography : Development of a new aquifer test method 36, 2095–2105. Yeh, T.-C.J., Liu, S., Glass, R.J., Baker, K., Brainard, J.R., Alumbaugh, D., Labrecque, D., 2002. A geostatistically based inverse model for electrical resistivity surveys and its applications to vadose zone hydrology. Water Resour. Res. 38, 1–13. https://doi.org/ 10.1029/2001WR001204. Yin, L., Hou, G., Dou, Y., Tao, Z., Li, Y., 2011. Hydrogeochemical and isotopic study of groundwater in the Habor Lake Basin of the Ordos plateau, NW China. Environ. Earth Sci. 64. https://doi.org/10.1007/s12665-009-0383-z. Yoon, H., Hart, D.B., Mckenna, S.A., 2013. Parameter estimation and predictive uncertainty in stochastic inverse modeling of groundwater flow : comparing null-space Monte Carlo and multiple starting point methods. Water Resour. Res. 49, 536–553. https://doi.org/10.1002/wrcr.20064. Zakari, A., Robert, N., Philippe, N., Jamal, A., 2015. Aquifers productivity in the PanAfrican context. J. Earth Syst. Sci. 124, 527–539. Zha, Y., Yeh, T.J., Illman, W.A., Onoe, H., Mok, C.M.W., Wen, J., Huang, S., Wang, W., 2017. Incorporating geologic information into hydraulic tomography: a general framework based on geostatistical approach. Water Resour. Res. 53, 2850–2876. https://doi.org/10.1002/2016WR019185. Zhang, L., Franklin, J.A., 1993. Prediction of water flow into rock tunnels: an analytical solution assuming an hydraulic conductivity gradient. Int. J. Rock Mech. Mining Sci. Geomech. Abstracts 30, 37–46. Zhao, J., 1998. Rock mass hydraulic conductivity of the Bukit Timah granite, Singapore. Eng. Geol. 50, 211–216. Zheng, C.M., Bianchi, M., Gorelick, S.M., 2011. Lessons learned from 25 years of research at the MADE site. Ground Water 49, 649–662. Zyvoloski, G., Kwicklis, E., Eddebbarh, A.A., Arnold, B., Faunt, C., Robinson, B.A., 2003. The site-scale saturated zone flow model for Yucca Mountain: calibration of different conceptual models and their impact on flow paths. J. Contam. Hydrol. 62, 731–750.
redistribution within an oak-pine forest from resistivity imaging. J. Hydrol. 430, 69–79. Saar, A., Magna, M.O., 2004. Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints. J. Geophys. Res. 109, B04204. https://doi.org/10.1029/2003JB002855. Sanford, W.E., 2017. Estimating regional-scale permeability–depth relations in a fractured-rock terrain using groundwater-flow model calibration. Hydrogeol. J. 25, 405–419. Schwarz, G., 1978. Estimating the Dimension of a Model. Ann. Stat. 6 (2). https://doi. org/10.1214/aos/1176344136. Sikandar, P., Christen, E.W., 2012. Geoelectrical sounding for the estimation of hydraulic conductivity of alluvial aquifers. Water Resour. Manage. 26, 1201–1215. Singh, K.P., 2005. Nonlinear estimation of aquifer parameters from surficial resistivity measurements. Hydrol. Earth Syst. Sci. Discuss. 2, 917–938. Singha, K., Gorelick, S.M., 2006. Effects of spatially variable resolution on field-scale estimates of tracer concentration from electrical inversions using Archie’s Law. Geophysics 71, G83–G91. Snow, D.T., 1968. Hydraulic character of fractured metamorphic rocks of the front range and implications of the rocky mountain arsenal well. Colorado School of Mines Quarterly 63, 167–199. Stober, I., Bucher, K., 2007. Hydraulic properties of the crystalline basement. Hydrogeol. J. 15, 213–224. Urumović, K., Sr, K.U., 2016. The referential grain size and effective porosity in the Kozeny-Carman model. Hydrol. Earth Syst. Sci. 20, 1669–1680. Wang, X.S., Jiang, X.W., Wan, L., Ge, S.M., Li, H.L., 2011. A new analytical solution of topography-driven flow in a drainage basin with depth-dependent anisotropy of permeability. Water Resour. Res. 47, 3101–3106. Wang, X.S., Jiang, X.W., Wan, L., Song, G., Xia, Q., 2009. Evaluation of depth-dependent porosity and bulk modulus of a shear using permeability–depth trends. Int. J. Rock Mech. Min. Sci. 46, 1175–1181. Waxman, M., Thomas, E., 1974. Electrical conductivities in shaly sands: i. Relation between hydrocarbon saturation and resistivity index; ii. The temperature coefficient of electrical conductivity. Spe Journal 12, 213–225. Williams, C.F., Narasimhan, T.N., 1989. Hydrogeologic constraints on heat flow along the
11
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Research papers
Depth-dependent relation between hydraulic conductivity and electrical resistivity in geologic formations ⁎
T
⁎
Kaixuan Lia, Tongchao Nana, , Xiankui Zenga, Lihe Yinb, Jichun Wua, , Jun Zhangb a b
Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China Xi’an Center of Geological Survey, China Geological Survey, China
A R T I C LE I N FO
A B S T R A C T
This manuscript was handled by Huaming Guo, Editor-in-Chief
Predictability of groundwater flow and transport models greatly relies on adequate characterization of aquifer heterogeneity such as hydraulic conductivity (K), porosity (ϕ) etc. The traditional model combining Archie formula and Kozeny-Carman equation (thus called AKC model in this study) has been proven in a large number of studies to be helpful in evaluation of K through electrical resistivity (Re). However, AKC model does not take into account the widely-reported dependence of K on depth, thus often fails in deep formations (> 100 m), as was observed in our study. In order to consider the impact of depth on the K-Re relation, a new model (called LAKC model) was introduced, which combined Louis decay model with AKC model. The effectiveness of LAKC model was verified against K measurements from stratified pumping tests and applied in a real aquifer in northern Ordos Basin, China. The results revealed that LAKC model can characterize the depth-dependent K-Re relation properly and that calibrated LAKC model is capable of providing adequate K estimates based on detailed Re measurements in a similar context of lithology.
Keywords: Aquifer characterization Electrical resistivity Northern Ordos Basin Depth dependence Stratified pumping test
1. Introduction Regional groundwater flow is critical to understand various geologic processes, such as fluid flow in sedimentary basins (Jiang et al., 2009; Wang et al., 2011), surface-ground water interaction (Winter, 1978), heat and solute transport (Cardenas, 2007; Niederau et al., 2017) etc. Aquifer heterogeneity is one of the key factors of properly modeling groundwater flow and solute transport. Unfortunately, due to technical and financial limitations, it is extremely impractical to reveal the spatial variability of real aquifers using traditional approaches like pumping tests and drilling (Binley et al., 2015; Nan and Wu, 2017; Yeh and Liu, 2000). It is still challenging for hydrogeologists to collect sufficient and reliable information to infer aquifer heterogeneity such as detailed distributions of hydraulic parameters (e.g. hydraulic conductivity, K and porosity, ϕ). The predictability and reliability of groundwater models (especially those of solute migration) can thus be significantly undermined in highly heterogeneous media (Mizukami et al., 2017; Nan et al., 2019; Yoon et al., 2013; Zheng et al., 2011). In recent decades the emergence of and advances in hydrogeophysics have made it possible to collect information of aquifer heterogeneity with growing amounts and accuracy, and even to estimate hydraulic parameters quantitatively (Dhakate and Singh, 2005; Frohlich, 1994; Niwas and de Lima, 2003; Rezaei et al., 2016; Zakari ⁎
et al., 2015). The geophysical approaches provide fast and effective technical supports for aquifer exploration and evaluation, and have become alternatives to traditional hydrogeological approaches to obtain hydrogeological parameters (Binley et al., 2015; Singh, 2005). The similarity between Ohm's law and Darcy’s law as well as linkage between geoelectric and hydraulic properties of formations made it possible to estimate hydrogeological parameters through geoelectric methods (Louis et al., 2004). There have been many studies on the relationship between electrical resistivity (Re) and hydrogeological parameters in sandy media, among which the Archie formula had been accepted by many researchers (Archie, 1942; Dick et al., 2018; Robinson et al., 2012; Singha and Gorelick, 2006; Whitman and Yeboah-Forson, 2015; Yeh et al., 2002; Zha et al., 2017). In a large number of studies Re obtained by geophysical methods were converted into hydraulic conductivity (K) through a model combining Archie equation and Kozeny-Carman equation (termed “AKC model” for brevity), and the performance of AKC model was acceptable in shallow formations (Kazakis et al., 2016; Niwas and Celik, 2012; Sikandar and Christen, 2012). Though in near-surface environments (< 1km) heterogeneity (e.g. lithofacies) is a controlling factor of permeability (or equivalently, K) (Gleeson et al., 2011; Ranjram et al., 2014), permeability and K were also found dependent on locale depth. In fact, a number of studies
Corresponding authors. E-mail addresses: [email protected] (T. Nan), [email protected] (J. Wu).
https://doi.org/10.1016/j.jhydrol.2019.124081 Received 13 May 2019; Received in revised form 23 August 2019; Accepted 26 August 2019 Available online 27 August 2019 0022-1694/ © 2019 Elsevier B.V. All rights reserved.
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
fine sandstone and pelite. Since no continuous aquifuge layer exists above J2a, the formations within 800 m deep can be treated as a unified aquifer of large thickness (over 800 m) (Yin et al., 2011). The total depths of boreholes ZKH2, ZKH3 and ZKH5 are 958 m, 703 m and 753 m, respectively. Three boreholes expose Quaternary Holocene aeolian sand layer (Q4eol) at depth 0 m, 0–10 m and 0–1.6 m, respectively. Huanhe formation (K1h) is found in the three boreholes at 0–474 m, 10–330 m and 1.6–574 m, respectively. Luohe formation (K1l) are 474–880 m, 330–701 m and 574–753 m, respectively. Anding formation (J2a) is exposed by ZKH2 and ZKH3 at 880 m and 701 m (see Fig. 2). In this study, XSKC-92 logging instrument was used for well logging. By using this instrument, the borehole temperature and measured resistivity (R0) values of depth 5–700 m were collected from ZKH5, ZKH2 and ZKH3 boreholes, with electrode space equal to 5 cm. The range of temperature in three boreholes at all depths are 19–30 °C, and the average is about 24 °C. Considering the influence of ground temperature variation on R0 within 700 m depth, the temperature correction for R0 is carried out by Arps' formula (Arps, 1953):
demonstrated that permeability generally decays with depth (Jiang et al., 2009; Louis, 1974; Manning and Ingebritsen, 1999; Ranjram et al., 2014; Snow, 1968; Wang et al., 2009). Such decay behaviors can be remarkable in thick aquifers and significantly change the conditions of regional groundwater flow systems (Jiang et al., 2010; Saar and Magna, 2004; Zhao, 1998). But quantitative descriptions of depth-related decay is still limited (Sanford, 2017). There were studies on depth’s impact on permeability through empirical or semi-empirical methods (see Jiang et al., 2011; Kuang and Jiao, 2014; and references therein). To describe the general relationship between permeability and depth a few models were proposed, among which two major types were widely used, i.e. the power-law model (Bredehoeft et al., 1992; Ingebritsen and Manning, 2002, 2010; Manning and Ingebritsen, 1999; Stober and Bucher, 2007) and the exponential model (Jiang et al., 2010; Louis, 1974; Saar and Magna, 2004; Williams and Narasimhan, 1989). These models were aimed to a general rule of depth-permeability relation so that intrinsic heterogeneity of medium (such as lithofacies) was ignored and variations in observed permeability data in a large range of depth were fully attributed to depth (Kuang and Jiao, 2014). In thick and highly heterogeneous aquifers, it is better to take into account the impacts of both formation heterogeneity and depth. In this paper, we develop a depth-dependent K-Re model which combines a depth-permeability model of Louis (1974) exponential type and widely used AKC model (hence the proposed model is called LAKC model). Rather than analyzing the impact of gravity stress (depth) on each AKC parameter or intermediate separately, a lumped exponential function is used to fully account for gravity impact to the final K estimates. Our goal is to consider the impacts of aquifer heterogeneity and depth simultaneously so that K in thick regional aquifers can be estimated by Re measurements collected by geophysical methods in a parsimonious manner. To analyze the applicability of LAKC model, we investigate the performance of LAKC model in a cross-validation experiment based on K and Re measurements collected in three boreholes of depth > 700 m in a thick aquifer. That is, LAKC model is calibrated via PEST code using measurements from any two boreholes and validated using measurements from the third borehole. To demonstrate the necessity of including depth dependence, AKC model is also applied in a similar procedure and compared to LAKC model. After validation of LAKC model, we apply it in estimating detailed K values from 5 cm-spacing Re data throughout the boreholes. This paper is structured as follows. Section 2 explains AKC model, LAKC model and calibration procedures. Section 3 presents the information of the study area and available data. In Section 4, we analyze the predictability of LAKC model in a cross-validation experiment compared to AKC model and apply LAKC model throughout depth to obtain detailed K distributions and discuss depth-dependency of effective porosity. Section 5 provides a summary and the main conclusions.
Re =
T0 + 21.5 R0 Te + 21.5
(1)
where T0 (°C) is the borehole temperature, R0 (Ω·m) is the measured resistivity and Re (Ω·m) is the electrical resistivity at a base-line temperature Te (equal to average temperature in three boreholes). Stratified pumping tests based on Solexperts Heavy-Duty Double Packer System (HDDP) were conducted in 44 borehole segments of length 20 m or 50 m. HDDP is a set of drilling water sample acquisition system consisting of two balloon embolisms (packers). In this method, two packers were inflated with high-pressure nitrogen, and in close contact with the well wall to produce water-stopping effect. The nontarget sections at both ends of the target zone (segment) were isolated through Packer air bag, and then the target zone was pumped by submersible pump. Meanwhile, the groundwater level of the target zone (segment) was monitored and recorded in real time by a probe. Groundwater flows were approximated as radial flow in horizontal direction in data interpretation. The interpreted (measured) 44 K values represent the mean horizontal K on 44 segments and were shown in Fig. 3. The interpreted K varies over three orders of magnitude, showing high heterogeneity as well as a plausible decay trend with depth. Total dissolved solids (TDS) was measured in laboratory based on water samples collected in stratified pumping tests. Rw profiles (at average temperature 24 °C) were calculated based on interpolated TDS measurements (e.g. Fig. 5(d) below) and reciprocal relation of TDS-Rw. The relevant data can be found in Nan (2019). 3. Materials and methodology 3.1. AKC model: a classic K-Re relation based on Archie and KozenyCarman formulae
2. Study area and available data Archie formula reflects the relationship between (electrical) resistivity of granular medium and medium porosity, which reads as:
The hydrogeological exploration boreholes ZKH5, ZKH2 and ZKH3 in this study are located in the middle and low desert plain of Wushenqi area, the lake-concentrated area in northern Ordos Basin, China (Fig. 1). The strata in this area are almost horizontal, and all the strata are terrigenous clastic deposits. The hydrogeological profile of the study area is shown in Fig. 2. The aquifer under study mainly consists of Lower Cretaceous Luohe formation (K1l) and Lower Cretaceous Huanhe formation (K1h). The lithofacies are semi-cemented siltstone, fine-/medium-/coarse-sandstone, gritstone and sandwiched silty mudstone, with well-developed pores and low shale content. At the top is discontinuously-distributed Quaternary Holocene aeolian sand layer (Q4eol, thickness < 10 m). The underlying formation below K1l is Jurassic Anding formation (J2a, depth > 800 m in most areas), consisting of low-permeability dense
R e = αRw ϕ−mS
−n
(2)
where Re is volume resistivity of rock matrix; Rw is resistivity of fluid (water in aquifers); ϕ is the porosity of the medium; n is saturation index; m is cementation index, which is related to lithology, pore structure and diagenesis; α is the value related to the aqueous medium; S is the saturation of the medium, and for the saturated aquifer, S = 1. In general, the values of α and m are obtained by testing core samples in laboratory and then applied to formations of similar lithology (Kazakis et al., 2016). If there is a lack of drilling samples in the study area, empirical values can be given according to lithofacies and documented ranges in literature (Aristodemou and Thomas-Betts, 2000). As mentioned above, Re and Rw data in this study are all corrected to base-line 2
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 1. Sketch map of the study area and borehole locations.
thus ignored due to small temperature variations (< 11 °C) in three boreholes; g is the gravitational acceleration 9.81(m/s2); de is the particle size (mm), set to be average particle size of the corresponding lithological section according to the bore histogram of three boreholes. The range of de used in this study is listed in Table 1. The classical AKC model is obtained by combining Eqs. (2) and (3) to eliminate porosity ϕ, that is,
temperature 24 °C. Note that in formations of low clay content, temperature has similar impacts on Re and Rw and hence the ratio Re/Rw can be considered temperature-independent (Waxman and Thomas, 1974). Though temperature correction is not really required to apply Eq. (2) in this area, these corrections are still carried out in both Re and Rw data, simply in order to filter out temperature-induced trend in data visualization (e.g. Re and TDS profiles in Fig. 5 below). In porous media ϕ is closely related to K. The pore characteristics (such as number, size, shape, direction, degree of connectivity and spatial changes) determines the soil storage, retention, release and transmission performance of water (Niu et al., 2015). For two similar sandy media, one with a higher ϕ usually has a higher K (Sikandar and Christen, 2012). Kozeny-Carman formula can be expressed as: (Carman, 1937; Kozeny, 1927)
KK−C =
ρw gde2 ϕ3 180μ w (1 − ϕ)2
KAKC =
ρw gde2 (αR e / Rw )1/ m · 180μ w [(Rw / αR e )1/ m − 1]2
(4)
KAKC represents K value calculated by AKC model. Traditionally, AKC model is used to convert rock resistivity collected by geophysical method into K. This method of estimating K by ϕ has been found effective in granular media with low clay content (Flint and Selker, 2003; Urumović and Sr, 2016).
(3)
3.2. LAKC model: a new K-Re relation based on Louis negative exponential model
where KK-C is calculated hydraulic conductivity (m/s); ρw is water density, equal to 1000 kg/m3; μw is dynamic viscosity of water, equal to 0.0014 kg/m·s in this study (Fetter, 1994). Temperature dependence of water density and viscosity are very small (Al-Shemmeri, 2012) and
Although AKC model is useful in estimating K in shallow aquifers (Kazakis et al., 2016; Niwas and Celik, 2012; Sikandar and Christen,
Fig. 2. Hydrogeological profile of the study area. Red sections indicate the tested segments in stratified pumping tests. Dash lines show formation division. 3
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 3. Depth cover and interpreted K values at 44 test segments in stratified pumping test at ZKH5, ZKH2 and ZKH3.
K decays with depth. Researchers tried to find one or several universal attenuation rate(s) for Earth’s crust (Kuang and Jiao, 2014; and references therein), but strictly speaking, A actually depends on compressibility of the medium and should have different values for different lithofacies. Fortunately, the variability of rock compressibility is much weaker than K (Domenico and Mifflin, 1965) so that one can approximately use a single A value within a stratum without introducing remarkable errors. In this study, two individual values of A are separately calibrated for Huanhe formation (K1h) and Luohe formation (K1l).
2012), it does not consider K’s decay with depth and generally fails in thick aquifers where differences in self-weight stress are significant. Here we employ a negative exponential model which was proposed based on statistical data by Louis (1974) and then widely accepted in applications (Anderman and Hill, 2003; Zhang and Franklin, 1993). Louis’ negative exponential model is as follows:
K (h) ≈ K 0·exp(−Ah)
(5)
where K0 is K at surface; K(h) is K at depth h; A is an attenuation coefficient. The larger A is, the faster K decreases with depth. Since AKC model works properly in near-surface sandy aquifers, AKC-based K estimate is taken as K0 in Eq. (5) here, i.e., K0 = KAKC. By combining Eqs. (4) and (5), one immediately obtains a depth-dependent model for K- Re relation i.e. Eq. (6) or (7), which we refer to as LAKC model:
K (h) ≈ KAKC·exp(−Ah)
K (h) =
3.3. Model calibration, cross-validation and error evaluation To investigate the performance of AKC and LAKC models, we calibrate the models against K and Re data from two boreholes and validate the models using data from the third borehole. Cross-validations are carried out for all three distinct calibration-validation combinations of boreholes, i.e. C1: (ZKH2 + ZKH5 | ZKH3), C2: (ZKH2 + ZKH3|ZKH5) and C3: (ZKH3 + ZKH5|ZKH2). A parameter-estimation software, PEST (Doherty, 2008), is used for automatic model calibration. PEST is independent of calculation models and has been widely used in groundwater and soil water transport models (see e.g. Fang et al., 2010;
(6)
ρw gde2 (αR e / Rw )1/ m · ·exp(−Ah) 180μ w [(Rw / αR e )1/ m − 1]2
(7)
While AKC has two parameters α and m, LAKC has one more parameter A to be determined. As mentioned above, A reflects how fast Table 1 The range of particle size (de) used in this study. Lithology
Gritstone
coarse sandstone
Medium sandstone
Fine sandstone
Siltstone
Clay
Particle diameter (de) (mm)
>1
0.5–1
0.25–0.5
0.05–0.25
0.002–0.05
< 0.002
4
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
parameter ranges does not improve the fitting much (not shown). The decay of K with depth has to be considered in these boreholes for better fitting.
Zyvoloski et al., 2003). The ranges and initial guesses of parameters are given according to empirical values in literature (Cai et al., 2017). That is, the ranges of parameters α and m are set to 0.48–1.5 and 1.1–3.0, and initial values are 1.0 and 1.2, respectively, in both models. The initial value and range of A is set to 0.016 m−1 and 0.003–0.03 m−1. Note that A has two independent values in Huanhe (K1h) and Luohe formations (K1l). The root mean square error (RMSE) is employed to evaluate the model performances. RMSE can be calculated by the following Eq. (8):
RMSE =
1 n
4.2. LAKC model Fig. 6 reports the cross-validation results of LAKC model. It can be seen that the estimation results seem reasonable and close to each other for different schemes of data usage. The estimation result suggests that calibrated LAKC model is consistent with all the data. With different value of α and m for three data usage, the RMSE for all data in C1, C2 and C3 data usage are 0.073, 0.124 and 0.079, respectively. Fig. 7 shows profiles of fitted and measured K values in three boreholes for C1 with parameter estimates α = 0.56, m = 1.60, A1 = 0.0155 m−1 and A2 = 0.0070 m−1. The agreement between measured and calculated K profiles confirms the validity of calibrated LAKC model. Similar agreements can be found for the other two groups of estimates thus omitted for brevity. To directly compare the goodness of fits for AKC and LAKC models, scatter plots of measured and calculated K values in two models for data usage C1 are shown in Fig. 8. Coefficient of determination R2 and RMSEs are also included. Fig. 8(a) and 8(b) report the results in calibration and validation stages, respectively. In calibration, RMSEs for AKC and LACK models are 1.410 and 0.068 respectively, and R2 relative to y = x line are −1.038 and 0.944. In validation, RMSEs for AKC and LAKC are 1.030 and 0.087, respectively and the corresponding R2 values are 0.902 and 0.999. Both high RMSE and low R2 indicate that AKC model is not proper to fit the data. A negative R2 in AKC model means that the trend of line y = x is inconsistent with trend in data. In contrast, very low RMSEs and high R2 show that fitted K values by LAKC model are close to measured K values and that LAKC model performs far better than AKC model in both calibration and validation. AIC and BIC values of AKC and LAKC models (number of parameters, p = 2 and 4, respectively) are calculated using Eqs. (9) and (10). AIC of AKC and LAKC are 28.2 and −221.8, respectively; BIC values of them are 31.8 and −214.6, respectively. Clearly, LAKC model has smaller values than AKC model. LAKC model performs much better than AKC model, not because of parameter number. The good performance of LAKC may result from two major reasons, i.e. (1) the ideal geological environment at the site and (2) the significance of gravity impact. The formations mainly consist of sandstones, siltstone and gritstone, and have very low content of clay. Thus, it provides a good situation in which Archie and Kozeny-Carman formulae properly apply. In fact, AKC model works fine in fitting measurements in shallow segments (depth < 100 m). When the entire profile is considered, gravity impact on the profile is primary and overwhelming among those factors which AKC doesn’t take into account. As discussed above, the attenuation coefficient A in LAKC model reflects the variation degree of K and relates to compressibility of media. In this study, the attenuation coefficient A1 in Huanhe formation is larger than A2 in Luohe formation. It implies the decrease of A with increase in depth and compaction degree of strata, which is consistent
n
∑ (Km (i) − K c (i))2 i=1
(8)
where Km is the measured K at the i-th test segment and Kc is the arithmetic mean of calculated K within the i-th test segment according to AKC or LAKC model, and n is the total number of test segments. RMSE reflects the overall deviation of the model from reference (data). The smaller RMSE is, the better the model fits the data. The Akaike information criterion(AIC) (Akaike, 1973) and Bayesian information criterion (BIC) (Schwarz, 1978) were also used to contrast the effects of LAKC model and AKC model in parameter estimation. The AIC and BIC can not only measure the deviation between measured data and model predicted data, but also reflect the instability caused by the different number of model parameters (Enemark et al., 2019). The smaller value of AIC and BIC is, the better the parameter estimation model performs. They are calculated as follows:
AIC = nln(ESS/n )+ 2p BIC = nln(ESS/ n) + p ln n
(9) (10)
ESS is the residual sum of squares between Km and the Kc calculated by AKC and LAKC model, n is the total number of test segments, and p is the number of model parameters. 4. Results and discussions 4.1. AKC model Fig. 4 reports calibration results in cross-validation experiment. In all calibrations, both parameters α and m of AKC model consistently stop at lower bounds of given ranges (i.e. α = 0.48, m = 1.1). Since estimated α and m are the same for C1, C2 and C3, the calibrated AKC models are identical and the resultant RMSEs are all equal to 1.317. It implies that no optimal fitting is found for AKC model within given parameter ranges. Fig. 5 reports fitted K profiles based on AKC model and estimated parameters (α = 0.48, m = 1.1) compared to measured K values. Re profiles are also included for comparison. The fitted K values tend to be much (1–2 orders of magnitude) larger than measured K values, especially at deep segments. The remarkable differences between measured and fitted K values confirm that AKC model fails in fitting the data for given parameter ranges without depth attenuation effect. Note that the given ranges of two parameters include all documented possible values for AKC model so far. Even further extension of
Fig. 4. The cross-validation results of AKC model. Final estimates are shown by solid lines inside boxes while top and bottom of boxes show 95% confidence intervals (CIs). Red line indicates initial values (IV); dash lines are upper and lower bounds (UB and LB) of parameter ranges. 5
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 5. (a–c) K values calculated by AKC model compared to measured K values in three boreholes for data usage C1. Profiles of resistivity Re are also included for comparison. (d) TDS data collected in ZKH3.
Fig. 6. The cross-validation results of LAKC model. Final estimates are shown by solid lines inside boxes while top and bottom of boxes show 95% confidence intervals (CIs). Red line indicates initial values (IV); dash lines are upper and lower bounds (UB and LB) of parameter ranges.
with the finding of Jiang et al. (2011). Parameters α and m are usually determined through core electricity experiments. In this study, the parameters were obtained in parameter estimation and then verified against the empirical values in similar lithotypes in literature. The results show that the estimated values of α and m in the study are close to those in similar lithotypes in literature (Carothers, 1968). In addition, although the Cretaceous Huanhe
Formation and Luohe Formation are two different sedimentary strata, they are lithologically similar. Therefore, the same set of parameters are used in both Huanhe Formation and Luohe Formation in this paper. 4.3. K profiles of 5-cm spacing based on Re data Based on Re data of 5-cm spacing in well log of three boreholes, 6
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 7. K values calculated by LAKC model compared to measured K values in three boreholes for data usage C1.
Fig. 8. Calculated K versus measured K values for data usage C1.
siltstone (e.g. at 103–134 m and 484–524 m in ZKH5, 83–87 m and 468–474 m in ZKH2, 306–313 m in ZKH3) or a thin layer of clay of thickness 0.5 m or less. Besides, the same type of rocks at different depths tend to have different K. For example, the averaged K of medium sandstone in ZKH2 at depth of 90–112 m, 333–376 m and 448–468 m are 0.327, 0.077 and 0.022 m/d, respectively. It suggests once again that apart from lithology, the influence of depth on K must be considered at large depth. A comparison of measured and the log K calculated by previous permeability-depth models, such as the power law model (Manning and Ingebritsen, 1999), exponential model (Saar and Magna, 2004) and Kuang-Jiao model (Kuang and Jiao, 2014), are made here (Fig. 10). The formula K = ρw gk/μw was used to convert permeability k (m2) to hydraulic conductivity K (m/d), and the ρw, g, μw take the same value as
estimated K profiles of 5-cm spacing from AKC model and LAKC models are shown in Fig. 9. For AKC model, α = 0.48, m = 1.1 according to parameter estimation above; for LAKC model, α = 0.56, m = 1.60, A1 = 0.0155 m−1 and A2 = 0.0070 m−1. It is found that despite of scale difference, K profiles from LAKC model (blue line) in three boreholes are all compatible with measured K values at segments (red line) in trend but reveal much more details of heterogeneity. On the other hand, K profiles from AKC model (black line) and measured K profiles are quite divergent, especially the trend in ZKH5 and ZKH2. According to Fig. 9 and lithology records in boreholes, strata with relatively high K generally corresponds to coarse sandstone and medium coarse sandstone, such as at (385–485 m) in ZKH5, (113 m–320 m) in ZKH2, and (120–252 m), (362–513 m) in ZKH3. And depth with low K often corresponds to fine sandstone, argillaceous 7
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 9. Hydraulic conductivity curve of borehole ZKH5/2/3 for data usage C1. Red line: K measured in pumping test; blue line: K calculated by LAKC model; black line: K calculated by AKC model.
contribution of pore water to the electrical resistivity, while KozenyCarman formula invokes an effective porosity ϕe for seepage. ϕe and ϕ are close to each other in shallow sandy aquifers and can be eliminated as intermediate variables in AKC model. But in aquifers of large depth, ϕe and ϕ may be remarkably different and one doesn’t have AKC model (Eq. (4)). If we assume that Kozeny-Carman formula is still valid, and denote
the parameters in AKC and LAKC model. As shown in Fig. 10, these existing permeability-depth models can reflect the attenuation characteristics of K with depth, but the accuracy is far from enough for the study of aquifer heterogeneity. Furthermore, it is difficult to describe the outlier values of K (data points in red circles in Fig. 10), which actually impact groundwater migration significantly. It can be seen that general permeability-depth models are not good enough in aquifer characterization and that local heterogeneity has to be included. The RMSE and R2 for different models are also shown in Fig. 10. The method of calculating average log K in different models is the same as that used in subsections above. The RMSE of exponential model is the lowest in three permeability-depth models, and the RMSE of power law model is the highest because of its spurious large values near the surface, which is an inherent defect of the model (Kuang and Jiao, 2014). The exponential and Kuang-Jiao models seem better at fitting the data except four “outliers” in red circles (ZKH3), yielding moderate R2 values. AKC model does not consider depth and has a result as bad as those of depth-only models or worse. Negative values of R2 which appear in fitting of AKC model and power law model, reflects that the trend of log K calculated by these two models is not consistent with that of measured K. LAKC model has the lowest RMSE and highest R2 among all the models for all data combinations. It shows that LAKC model considering both depth and lithology has higher accuracy than those models considering only one single factor.
the effective porosity as ϕe, then K (h) = LAKC model, one gets:
K (h) =
ρ w gde2
·
ϕe3
180μ w (1 − ϕe)2
ϕe3 ρw gde2 ρ gde2 ϕ3 · · exp(−Ah) = w 180μ w (1 − ϕe)2 180μ w (1 − ϕ)2
. Upon adopting
(11)
Eq. (11) allows one to calculate ϕe and ϕ. Fig. 11 reports the results of ϕe and ϕe/ϕ with depth in ZKH2. The jumps in both subplots at formation interface (depth 475 m) are related to our selection of independent A parameter for Luohe and Huanhe formations. Fig. 11 shows that despite discontinuity at formation interface, both ϕe and ϕe/ ϕ have a decreasing trend. ϕe/ϕ becomes far below 50% for large depth, which is consistent with our analysis of the possible reason of AKC failure above. 5. Conclusions Aquifer heterogeneity has significant influences on groundwater flow field, solute and heat transport. The traditional AKC model which relates K to resistivity Re does not consider the impact of depth on K-Re relation. According to stratified pumping tests at 44 segments in three boreholes, K of shallow strata (< 100 m) is much (1–2 orders of magnitude) larger than that in deep strata. Based on the traditional AKC
4.4. Depth dependency of effective porosity ϕe and ratio ϕe/ϕ A possible explanation of AKC model’s failure is the large difference between total porosity ϕ and effective porosity ϕe in thick aquifers. Archie formula requires a total porosity ϕ to fully account the 8
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
Fig. 10. Comparison of measured log K and log K calculated by models proposed by predecessors (In these models, k is the permeability (m2) and z is the depth (km)). (“All data” represents the all 44 K values measured in ZKH5, ZKH2, ZKH3.)
Fig. 11. The effective porosity ϕe and ϕe/ϕ in ZKH2. 9
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
model and the newly established LAKC model, K within an aquifer of 700 m deep was estimated using Re data from well logs. In the attempt of fitting data from stratified pumping tests using AKC model, it was found that AKC model was unable to fit the data to an acceptable extent, especially for the decay trend in K. Though depth probably also impacts rock resistivity, standard AKC model seems inappropriate in thick aquifers and a depth-based modification to AKC model is necessary. As a depth-modified version of AKC model, LAKC model was able to match all the data very well and presented the decay trend in K with increasing depth. Results from cross-validation suggested that LAKC model calibrated by data from any two boreholes provides reliable K estimates for the third boreholes. According to lithology record and pumping tests in three boreholes, K values of similar lithofacies at different depths are quite different, showing a decay trend with depth. In order to infer K in thick aquifers based on borehole resistivity, one should take into account not only lithology-induced resistivity variation but also depth-led K decay simultaneously. This study demonstrated not only the potential of LAKC model in K estimation based on resistivity, but also the necessity to consider K decay with depth in quantitative description of heterogeneity in thick aquifers. It is also possible to reflect depth impact in AKC parameters (e.g. cementation index m) instead of introducing Louis decay function, which will be explored in our following study. In the future, more geophysical methods (such as the magnetotelluric, transient electromagnetic method, etc.) can be combined to transform resistivity data to K, which will be a promising way to obtain the spatial distribution of K with a larger scale and higher dimension.
Cardenas, M.B., 2007. Potential contribution of topography-driven regional groundwater flow to fractal stream chemistry: residence time distribution analysis of Tóth flow. Geophys. Res. Lett. 34. https://doi.org/10.1029/2006GL029126. Carman, P.C., 1937. Fluid flow through granular beds. Chem. Eng. Res. Des. 75, S32–S48. https://doi.org/10.1016/S0263-8762(97)80003-2. Carothers, J.E., 1968. A statistical study of the formation factor relation. Log Anal. 9, 13–20. Dick, J., Tetzlaff, D., Bradford, J., Soulsby, C., 2018. Using repeat electrical resistivity surveys to assess heterogeneity in soil moisture dynamics under contrasting vegetation types. J. Hydrol. 559, 684–697. Dhakate, R., Singh, V.S., 2005. Estimation of hydraulic parameters from surface geophysical methods, Kaliapani Ultramafic Complex, Orissa, India. J. Environ. Hydrol. 13, 1–11. Doherty, J., 2008. PEST, model independent parameter estimation user manual. Watermark Computing, Corinda, Australia. Domenico, P.A., Mifflin, M.D., 1965. Water from low-permeability sediments and land subsidence. Water Resour. Res. 1, 563–576. Enemark, T., Peeters, L.J.M., Mallants, D., Batelaan, O., 2019. Hydrogeological conceptual model building and testing: a review. J. Hydrol. 569. https://doi.org/10. 1016/j.jhydrol.2018.12.007. Fang, Q.X., Green, T.R., Ma, L., Erskine, R.H., Malone, R.W., Ahuja, L.R., 2010. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions. Soil Sci. Soc. Am. J. 74, 1897–1913. https://doi.org/10.2136/sssaj2009.0380. Fetter, C.W., 1994. Applied Hydrogeology, third ed. Prentice-Hall Inc., New Jersey. Flint, L.E., Selker, J.S., 2003. Use of porosity to estimate hydraulic properties of volcanic tuffs. Adv. Water Resour. 26, 561–571. Frohlich, R.K., 1994. The electric-hydraulic relationship. A geophysical model. Trends Hydrogeol. 1, 347–358. Gleeson, T., Smith, L., Moosdorf, N., Hartmann, J., Dürr, H.H., Manning, A.H., Van Beek, L.P.H., Jellinek, A.M., 2011. Mapping permeability over the surface of the Earth. Geophys. Res. Lett. 38, 93–104. Ingebritsen, S.E., Manning, C.E., 2010. Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism. Geofluids 10, 193–205. Ingebritsen, S.E., Manning, C.E., 2002. Diffuse fluid flux through orogenic belts: implications for the world ocean. PNAS 99, 9113–9116. Jiang, X.W., Wan, L., Wang, X.S., Ge, S.M., Liu, J., 2009. Effect of exponential decay in hydraulic conductivity with depth on regional groundwater flow. Geophys. Res. Lett. 36, 88–113. Jiang, X.W., Wang, X.S., Li, W., Ge, S., 2011. An analytical study on stagnant points in nested flow systems in basins with depth-decaying hydraulic conductivity. Water Resour. Res. 47, 128–139. Jiang, X.W., Wang, X.S., Wan, L., 2010. Semi-empirical equations for the systematic decrease in permeability with depth in porous and fractured media. Hydrogeol. J. 18, 839–850. Kazakis, N., Vargemezis, G., Voudouris, K.S., 2016. Estimation of hydraulic parameters in a complex porous aquifer system using geoelectrical methods. Sci. Total Environ. 550, 742–750. Kozeny, J., 1927. Uber Kapillare Leitung des Wasser im Boden, Wien Akad. Wiss. Sitzungsber. Akad.Wiss. Wien, Math. Naturw Klasse, Abt.II A. Kuang, X., Jiao, J.J., 2014. An integrated permeability-depth model for Earth“s crust. Geophys. Res. Lett. 41, 7539–7545. Louis, C., 1974. Rock Hydraulics BT – Rock Mechanics. Springer Vienna, Vienna, pp. 299–387. https://doi.org/10.1007/978-3-7091-4109-0_16. Louis, I.F., Karantonis, G.A., Voulgaris, N., Louis, F., 2004. The contribution of geophysical methods in the determination of aquifer parameters: the case of Mornos River delta, Greece. Res. J. Chem. Environ 8, 41–49. Manning, C.E., Ingebritsen, S.E., 1999. Permeability of the continental crust: implications of thermal data and metamorphic systems. Rev. Geophys. 37, 127–150. Mizukami, N., Clark, M., Newman, A.J., Wood, A.W., Samaniego, L., 2017. Toward seamless large domain parameter estimation for hydrologic models. Water Resour. Res. 53. Nan, T., 2019, Hydraulic conductivity and resistivity measurements in boreholes in Ordos, China, OSF, https://doi.org/10.17605/OSF.IO/78CQW. Nan, T., Wu, J., 2017. Application of ensemble H-infinity filter in aquifer characterization and comparison to ensemble Kalman filter. Water Sci. Eng. 10 (1), 25–35. https://doi. org/10.1016/j.wse.2017.03.009. Nan, T., Xue, L., Wu, J., 2019. Efficient identification of preferential flow path in heterogeneous media based on stream function. J. Hydrol. 577, 123961. Niederau, J., Ebigbo, A., Marquart, G., Arnold, J., Clauser, C., 2017. On the impact of spatially heterogenous permeability on free convection in the Perth Basin, Australia. Geothermics 66, 119–133. Niu, Q., Fratta, D., Wang, Y.H., 2015. The use of electrical conductivity measurements in the prediction of hydraulic conductivity of unsaturated soils. J. Hydrol. 522, 475–487. Niwas, S., Celik, M., 2012. Equation estimation of porosity and hydraulic conductivity of Ruhrtal aquifer in Germany using near surface geophysics. J. Appl. Geophys. 84, 77–85. Niwas, S., de Lima, O.A., 2003. Aquifer parameter estimation from surface resistivity data. Groundwater 41, 94–99. Ranjram, M., Gleeson, T., Luijendijk, E., 2014. Is the permeability of crystalline rock in the shallow crust related to depth, lithology or tectonic setting? Geofluids 15, 106–119. Rezaei, M., Saey, T., Seuntjens, P., Joris, I., Boënne, W., Van Meirvenne, M., Cornelis, W., 2016. Predicting saturated hydraulic conductivity in a sandy grassland using proximally sensed apparent electrical conductivity. J. Appl. Geophys. 126, 35–41. Robinson, J., Slater, L., Schaefer, L., 2012. Evidence for spatial variability in hydraulic
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 41730856, 41602250, U1503282) and the project “Hydrogeological investigation at 1:50 000 scale in the lakeconcentrated areas of the northern Ordos Basin” of China Geological Survey (Grant No. DD20160293). Data in all figures will be available upon publication of this paper. We declare that we have no conflict of interest. References Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Second International Symposium on Information Theory. Budapest, pp. 261–304. Al-Shemmeri, T., 2012. Engineering Fluid Mechanics. Ventus Publishing ApS, Holland. Anderman, B.E.R., Hill, M.C., 2003. MODFLOW-2000, the U.S. Geological Survey modular ground-water model-three additions to the hydrogeologic-unit flow (HUF) package: alternative storage for the uppermost active cells (SYTP parameter type), flows in hydrogeologic units, and the hydraulic-con. U.S. Geological Survey OpenFile Report 03-347. Archie, G., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 146, 54–62. Aristodemou, E., Thomas-Betts, A., 2000. DC resistivity and induced polarisation investigations at a waste disposal site and its environments. J. Appl. Geophys. 44, 275–302. Arps, J.J., 1953. The effect of temperature on the density and electrical resistivity of sodium chloride solutions. J. Petrol. Technol. 5, 17–20. https://doi.org/10.2118/ 953327-G. Binley, A., Hubbard, S.S., Huisman, J.A., Revil, A., Robinson, D.A., Singha, K., Slater, L.D., 2015. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res. 51, 3837–3866. Bredehoeft, J.D., Belitz, K., Sharp-Hansen, S., 1992. The hydrodynamics of the Big Horn Basin: a study of the role of faults. Am. Assoc. Pet. Geol. Bull. 76, 530–546. Cai, J., Wei, W., Hu, X., Wood, D.A., 2017. Electrical conductivity models in saturated porous media: a review. Earth Sci. Rev. 171, 419–433. https://doi.org/10.1016/j. earscirev.2017.06.013.
10
Journal of Hydrology 578 (2019) 124081
K. Li, et al.
San Andreas fault: a testing of hypotheses. Earth Planet. Sci. Lett. 92, 131–143. Winter, T.C., 1978. Numerical simulation of steady state three-dimensional ground water flow near lakes. Water Resour. Res. 14, 245–254. Whitman, D., Yeboah-Forson, A., 2015. Electrical resistivity and porosity structure of the upper Biscayne Aquifer in Miami-Dade County, Florida. J. Hydrol. 531, 781–791. Yeh, T.-C.J., Liu, S., 2000. Hydraulic tomography : Development of a new aquifer test method 36, 2095–2105. Yeh, T.-C.J., Liu, S., Glass, R.J., Baker, K., Brainard, J.R., Alumbaugh, D., Labrecque, D., 2002. A geostatistically based inverse model for electrical resistivity surveys and its applications to vadose zone hydrology. Water Resour. Res. 38, 1–13. https://doi.org/ 10.1029/2001WR001204. Yin, L., Hou, G., Dou, Y., Tao, Z., Li, Y., 2011. Hydrogeochemical and isotopic study of groundwater in the Habor Lake Basin of the Ordos plateau, NW China. Environ. Earth Sci. 64. https://doi.org/10.1007/s12665-009-0383-z. Yoon, H., Hart, D.B., Mckenna, S.A., 2013. Parameter estimation and predictive uncertainty in stochastic inverse modeling of groundwater flow : comparing null-space Monte Carlo and multiple starting point methods. Water Resour. Res. 49, 536–553. https://doi.org/10.1002/wrcr.20064. Zakari, A., Robert, N., Philippe, N., Jamal, A., 2015. Aquifers productivity in the PanAfrican context. J. Earth Syst. Sci. 124, 527–539. Zha, Y., Yeh, T.J., Illman, W.A., Onoe, H., Mok, C.M.W., Wen, J., Huang, S., Wang, W., 2017. Incorporating geologic information into hydraulic tomography: a general framework based on geostatistical approach. Water Resour. Res. 53, 2850–2876. https://doi.org/10.1002/2016WR019185. Zhang, L., Franklin, J.A., 1993. Prediction of water flow into rock tunnels: an analytical solution assuming an hydraulic conductivity gradient. Int. J. Rock Mech. Mining Sci. Geomech. Abstracts 30, 37–46. Zhao, J., 1998. Rock mass hydraulic conductivity of the Bukit Timah granite, Singapore. Eng. Geol. 50, 211–216. Zheng, C.M., Bianchi, M., Gorelick, S.M., 2011. Lessons learned from 25 years of research at the MADE site. Ground Water 49, 649–662. Zyvoloski, G., Kwicklis, E., Eddebbarh, A.A., Arnold, B., Faunt, C., Robinson, B.A., 2003. The site-scale saturated zone flow model for Yucca Mountain: calibration of different conceptual models and their impact on flow paths. J. Contam. Hydrol. 62, 731–750.
redistribution within an oak-pine forest from resistivity imaging. J. Hydrol. 430, 69–79. Saar, A., Magna, M.O., 2004. Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints. J. Geophys. Res. 109, B04204. https://doi.org/10.1029/2003JB002855. Sanford, W.E., 2017. Estimating regional-scale permeability–depth relations in a fractured-rock terrain using groundwater-flow model calibration. Hydrogeol. J. 25, 405–419. Schwarz, G., 1978. Estimating the Dimension of a Model. Ann. Stat. 6 (2). https://doi. org/10.1214/aos/1176344136. Sikandar, P., Christen, E.W., 2012. Geoelectrical sounding for the estimation of hydraulic conductivity of alluvial aquifers. Water Resour. Manage. 26, 1201–1215. Singh, K.P., 2005. Nonlinear estimation of aquifer parameters from surficial resistivity measurements. Hydrol. Earth Syst. Sci. Discuss. 2, 917–938. Singha, K., Gorelick, S.M., 2006. Effects of spatially variable resolution on field-scale estimates of tracer concentration from electrical inversions using Archie’s Law. Geophysics 71, G83–G91. Snow, D.T., 1968. Hydraulic character of fractured metamorphic rocks of the front range and implications of the rocky mountain arsenal well. Colorado School of Mines Quarterly 63, 167–199. Stober, I., Bucher, K., 2007. Hydraulic properties of the crystalline basement. Hydrogeol. J. 15, 213–224. Urumović, K., Sr, K.U., 2016. The referential grain size and effective porosity in the Kozeny-Carman model. Hydrol. Earth Syst. Sci. 20, 1669–1680. Wang, X.S., Jiang, X.W., Wan, L., Ge, S.M., Li, H.L., 2011. A new analytical solution of topography-driven flow in a drainage basin with depth-dependent anisotropy of permeability. Water Resour. Res. 47, 3101–3106. Wang, X.S., Jiang, X.W., Wan, L., Song, G., Xia, Q., 2009. Evaluation of depth-dependent porosity and bulk modulus of a shear using permeability–depth trends. Int. J. Rock Mech. Min. Sci. 46, 1175–1181. Waxman, M., Thomas, E., 1974. Electrical conductivities in shaly sands: i. Relation between hydrocarbon saturation and resistivity index; ii. The temperature coefficient of electrical conductivity. Spe Journal 12, 213–225. Williams, C.F., Narasimhan, T.N., 1989. Hydrogeologic constraints on heat flow along the
11