Modeling 3-D permeability distribution in alluvial fans using facies architecture and geophysical acquisitions

. Alluvial fans are highly heterogeneous in hydraulic properties due to complex depositional processes, which make it difﬁcult to characterize the spatial distribution of the hydraulic conductivity ( K ). An original methodology is developed to identify the spatial statistical parameters (mean, variance, correlation range) of the hydraulic conductivity in a three-dimensional (3-D) setting by using geological and geophysical data. More speciﬁcally, a large number of inexpensive vertical electric soundings are integrated with a facies model developed from borehole lithologic data to simulate the log 10 (K) continuous distributions in multiple-zone heterogeneous alluvial megafans. The Chaobai River alluvial fan in the Beijing Plain, China, is used as an example to test the proposed approach. Due to the non-stationary property of the K distribution in the alluvial fan, a multiple-zone parameterization approach is applied to analyze the conductivity statistical properties of different hydrofacies in the various zones. The composite variance in each zone is computed to describe the evolution of the conductivity along the ﬂow direction. Consistently with the scales of the sedimentary transport energy, the results show that conductivity variances of ﬁne sand, medium-coarse sand, and gravel decrease from the upper (zone 1) to the lower (zone 3) portion along the ﬂow direction. In zone 1, sediments were moved by higher-energy ﬂooding, which induces poor sorting and larger conductivity variances. The composite variance conﬁrms this feature with statistically different facies from zone 1 to zone 3. The results of this study provide insights to improve our understanding on conductivity heterogeneity and a method for characterizing the spatial distribution of K in alluvial fans.


Introduction
Alluvial fans usually house valuable groundwater resources because of significant water storage and favorable recharge conditions. Sedimentary processes forming alluvial fans are responsible for their complex long-term evolution. Usually, the coarsest material (gravel) is deposited in the upper fan, with the gravel passing into sand in the middle of the fan and then into silt and clay in the tail. A high heterogeneity characterizes the deposit distribution because of the shifting over time of the sediment-transporting streams (Zappa et al., 2006;. Hydraulic conductivity distributions in alluvial fans can be assigned according to the various hydrofacies simulated by conditional indicator geostatistical methods (Eggleston and Rojstaczer, 1998;Fogg et al., 1998;Weissmann et al., 2002a, b;Ritzi et al., 2004Ritzi et al., , 2006Proce et al., 2004;Dai et al., 2005;Harp et al., 2008;Hinnell et al., 2010;Maghrebi et al., 2015;Soltanian et al., 2015;Zhu et al., 2016a). However, the geostatistical methods require the stationary assumption; i.e., the distribution of the volumetric proportions and correlation lengths of hydrofacies converge to their mean values in the simulation do-main. The hydrofacies and hydraulic conductivity (K) distributions in alluvial fans are generally non-stationary (Weissmann et al., , 2013Anderson, 2007;Zhu et al., 2016a). Hence, the use of these methods may cause large characterization errors and add significant uncertainty to the predictions achieved by groundwater flow and contaminant transport models (Eggleston and Rojstaczer, 1998;Irving and Singha, 2010;Dai et al., 2014a). Zhu et al. (2016a) adopted a local-stationary assumption by dividing the alluvial fan into three zones along the flow direction of the Chaobai River, China. The zones were properly detected based on the statistical facies distribution. Then, the indicator simulation method was applied to each zone and the simulated hydrofacies distribution in the three zones was used to guide modeling of the K distribution.
Hydraulic conductivity of granular deposits generally varies with grain size, porosity and sorting. Traditional methods for K estimate, e.g., well test, permeability measurements, and grain-size analyses (Niwas et al., 2011), are very expensive, time-consuming and make it difficult to provide representative and sufficient field data for addressing spatial variations of conductivity. Recently, data fusion techniques have been developed for coupled inversion of multiplesource data to estimate K distributions for groundwater numerical modeling. Geophysical data (such as surface electric resistivity and various logging data) are relatively inexpensive and can provide considerable information for characterizing subsurface heterogeneous properties (Hubbard et al., 2001;Yeh et al., 2002;Dai et al., 2004a;Morin, 2006;Sikandar et al., 2010;Bevington et al., 2016). Electric resistivity data have been proven useful to derive sediment porosity distributions (Niwas and Singhal, 1985;Niwas et al., 2011;Niwas and Celik, 2012;Zhu et al., 2016b). Zhu et al. (2016b) simulated the spatial distributions of hydraulic conductivity by combining the interpolated resistivity on the basis of vertical electrical soundings (VESs) and the stochastic simulated facies through the empirical equation, in which the hydraulic conductivity was converted from the porosity data calculated from resistivity measurements and the grain size.
This study proposes a novel approach to reconstruct the three-dimensional (3-D) configuration of conductivity in alluvial fans by combining the hydrofacies spatial heterogeneity provided by a multiple-zone transition probability model with hydrogeological and hydrogeophysical measurements, in particular inexpensive VESs properly calibrated through resistivity logs acquired in a few wellbores. We assume the K distributions are local stationary; i.e., the mean and variance of log conductivity are convergent in each hydrofacies and in each local zone. Therefore, we can compute the log 10 (K) semivariogram in each hydrofacies and in each zone. The spatial structure features of hydraulic conductivity deduced from semivariograms are used during the geostatistical simulation processes of the hydraulic conductivity. The Chaobai alluvial fan (or "megafan" as defined by Leier et al., 2005, andHartley et al., 2010, for very large alluvial fans) in the northern Beijing Plain, China, was selected as study area to test the proposed integrated approach.

Study area
The study area belongs to the Chaobai River alluvial fan (or megafan), in the northern Beijing Plain (northern latitude 40-40 • 30 , eastern longitude 116 • 30 -117 • ), with an area of 1150 km 2 (Fig. 1a). The Chaobai River is the second largest river flowing through the Beijing Plain from north to south. The ground elevation decreases southward with an average 2 ‰ slope. Quaternary sediments were mainly deposited by flooding events with turbulent flow and consist of porous strata containing groundwater. The aquifer system in the alluvial fan can be divided into three zones according to the lithological features ( Fig. 1): an upper fan zone (or zone 1) with coarse sediments (e.g., sandy-gravel aquifers); a middle-upper fan zone (or zone 2), where medium-coarse sediments (e.g., sandy-gravel to sandy-silt aquifers) were deposited; and a fine sediment (e.g., sand and clay multiple aquifers) middle-lower fan zone (or zone 3). Four hydrofacies, including sub-clay and clay (C), fine sand (FS), medium-coarse sand (MS) and gravel (G), were classified based on the interpretations of the cores and textural description of almost 700 boreholes (Zhu et al., 2016a).
The study area is one of the most important regions for the supply of groundwater resource to Beijing. The Huairou emergency groundwater resource region (hereafter EGRR) with an area of 54 km 2 is located in zone 1. The total groundwater withdrawal amounted to 1.2 × 10 8 m 3 in 2003. Several well fields belonging to the "water supply factory" were drilled along the Chaobai River in zone 1 and the upper zone 2. Most of these well fields were built in 1979 with a designed groundwater pumping volume of 1.6 × 10 8 m 3 per year. The average thickness of the exploited aquifer system is approximately 300 m. The long-term over-exploitation of the aquifer system has resulted in a serious drawdown of water levels, which has reduced the exploitable groundwater resources and induced geological disasters, mainly land subsidence, fault reactivation and ground fissures (Cheng et al., 2015;Yang et al., 2015;Zhu et al., 2015). In 2010, the annual groundwater withdrawal at the EGRR and the water factory decreased to 0.86 × 10 8 and 0.65 × 10 8 m 3 , respectively.
The largest cumulative land subsidence from June 2003 to January 2010 was quantified in approximately 340 mm by Zhu et al. (2013Zhu et al. ( , 2015 in Tianzhu County to the south. The characterization of the distribution and spatial variability of the hydraulic conductivity is vital for an optimal use of the limited water resources in this area.

Methodological approach
Presently, a large set of hydraulic conductivity samples can be derived by integrating appropriate relations of various geological data, including hydrogeophysical measurements, borehole lithostratigraphies and hydrogeological information (total dissolved solid (TDS) and groundwater level). These databases can be statistically processed to derive the spatial variation of log 10 (K) for various facies, including clay, fine sand, medium-coarse sand and gravel.
In this paper, the statistical assessment is separately carried out for separated zones, building-up experimental semivariograms that are fitted with exponential models. The optimal parameters of the latter are estimated through a generalized output least-squares criterion. Then, the composite semivariograms are computed using a hierarchical sedimentary architecture Dai et al., 2005) to obtain the K variance in each zone. Finally, the configuration of log 10 (K) is simulated through a multiple-zone sequential Gaussian algorithm with estimated statistic parameters reflecting the K spatial structures in the alluvial fan. Figure 2 shows the steps involved in the developed approach.

Geophysical data
Geophysical data include resistivity logging and vertical electrical soundings. There are six well-electric logs continuously recording the formation resistivity versus depth. Five logs were collected in zone 2 and one in zone 3. Each well log has a lithological description, which helps to relate the resistivity values to the corresponding facies. The average resistivity of G is the largest, with a value of 198 m, and that of C is the smallest with a value of 24 m. Figure 3 compares the outcome of logging data in term of resistivity versus depth and the corresponding stratigraphy, where the groundwater depth is 12 m. The log was acquired in the eastern part of zone 2. The average resistivity from 32.4 to 40.5 m depth, where the sediments are mainly G and MS, is 70.8 m. The resistivity curve shows two evident peaks from 97 to 102 m and between 81 and 84.5 m depth, where the MS is located.
The C resistivity is relatively low due to the good intrinsic electrical conductivity of this facies. For example from 16.5 to 23.5 m depth, where C is the prevalent facies, a low resistivity equal to 27.2 m is recorded. Since a hydrofacies with a smaller grain size has a greater total surface area, the resistivity difference can partially reflect the distributions of particle sizes and the hydrofacies composition. Since the obtained resistivity is the apparent resistivity, we used the resistivity located in the middle of the facies block, where the resistivity is approximate to the real resistivity. Unfortunately, it is unknown to the authors if the logs were calibrated in the field and how the salinity of the formation water, although minimal and almost independent on the site and depth, has been accounted for. On the other hand, the resistivity distributions have good correlations with different hydrofacies along the vertical and horizontal directions. Therefore, in the mathematical framework that follows, we have assumed that the logs have been calibrated and are accurate enough for presenting our work as a proof of concept. VESs using the Schlumberger electrode configuration were carried out by the Beijing Institute of Hydrogeology and Engineering Geology (BIHEG). A number of 113 detecting positions were selected, with a maximum half current electrode space equal to 340 m and the potential electrode space ranging from 1 to 30 m. All the sounding data (1356 VES measurements) recorded the apparent resistivity of the porous medium. These data were inverted to real resistivity using the nonlinear Occam inversion method (Constable et al., 1987), with a low root mean square relative error of 2 %. Figure 4 shows the layered structure fitting model of resistivity and the borehole lithologic observations. The inversed resistivity generally reflects the difference of facies; the thick gravel layer has larger resistivity, whereas the finesand and clay layers have relatively smaller resistivity.

Geological and hydrogeological data
Almost 700 borehole lithologic logs were collected in the study area. The sedimentary deposits show large heterogeneity from the upper to the lower fan zone. In zone 1, the dominant facies is G with a volumetric proportion of 53 %. The volumetric proportion of C is 16 %. In zone 2, the volumetric proportion of C increases to 40%, while that of G decreases sharply to 24 %. In zone 3, the proportion of G decreases further to 6 % and that of C increases to 50 % (Table 1). More detailed information is given in Zhu et al. (2016a). The lithological information in a buffer zone of 200 m around the VES locations has been used to represent the actual facies distribution in the area surrounding the sites of the geophysical acquisitions.
A total of 35 hydrochemistry measurements with a depth from 20 to 270 m were obtained throughout the area. The minimum, maximum and average TDS values are 423 mg L −1 at the depth of 180 m, 943 mg L −1 at the depth of 50 m and 692 mg L −1 , respectively. Generally, the TDS is very low with the higher values measured in the southwestern part of the study area. Because of the relatively small dataset and the observed low variability, in this paper the TDS variation in the vertical direction has been neglected. A TDS map was obtained by interpolating the available records using an ordinary Kriging method with a spherical semivariogram model. A large number of depth of water level measurements were also collected to map the thickness of the unsaturated unit. The TDS and groundwater level at each VES and resistivity log location were derived from the interpolated surfaces.

Hydraulic conductivity estimates from geophysical acquisitions
The hydraulic conductivity K was estimated using the Kozeny-Carman equation: which is widely accepted to derive the hydraulic conductivity from grain size and porosity (Soupious et al., 2007;Utom et al., 2013;Khalil and Santos, 2013;Zhu et al., 2016b). In Eq. (1), d (x,y,z) is the median grain diameter (D 50 , mm) at the location (x, y, z), which was determined according to the lithology information (or lithological descriptions and grain-size distributions), g is gravity, µ the kinematic viscosity (kg/(m s −1 )), δ the fluid density and ϕ (x,y,z) the porosity. ϕ was estimated using Archie's law (Eq. 2), which relates the bulk resistivity of granular medium to porosity: where ρ is the saturated formation resistivity ( m), α the pore-geometry coefficient associated with the medium (0.5 ≤ α ≤ 2.5) and m the cementation factor (1.3 ≤ m ≤ 2.5) (Massoud et al., 2010;Khalil and Santos, 2013). α is set as 1.
In the upper part of alluvial (zone 1 and zone 2), m is set as 1.3 due to the sand being unconsolidated. In zone 3 m is set as 1.7, which reflects slightly cemented sandstones (Niwas et al., 2011); s w is the water saturation, and n the saturation index. The pore-fluid resistivity ( m) ρ w is calculated using the following experimental relation: with TDS in g L −1 , temperature t in • C, and b and β are constant parameters (Wu et al., 2003). For the most common electrolytes, b = −0.95 and β = 0.025. Note that the parameters associated with Eqs. (2) and (3) are site specific and applying these equations to other sites will mean re-adjusting the related parameters. The logarithmically transformed values of the estimated hydraulic conductivity (log 10 (K)) were used for the geostatistical analysis because of its normal distribution (Neuman, 1990). The histograms of log 10 (K) values within each facies are in Fig. 5. There are 102, 2077 and 1716 conductivity samples in zone 1, zone 2 and zone 3, respectively. Considering that Archie's law can only be used for clay-free granular sediments, the K values of C were not estimated in this study. Based on the lithological description information of borehole data, it has been reasonably assumed that clay fraction is negligible in G, MS, and FS facies. The statistics of log 10 (K) for the three facies in three zones are listed in Table 2. The mean log 10 (K) values decrease from zone 1 to zone 3, consistently with the sedimentary transport processes in the alluvial fan. In the upper region (zone 1), high waterflowing energy made the deposits consist mainly of largergrained particles and the coarse-grained sediments are dominant. In the southern part (zone 3), the deposits change to relatively fine-grained particles. The mean log 10 (K) of gravel is greater than 2.4 (lg(m/d)) and that of fine sand is less than 0.2 (lg(m/d)). The lithological information at the depth of the conductivity samples shows that volumetric proportions of FS and MS increase and that of G decreases from zone 1 to zone 3. The results are consistent with the statistic outputs deduced from 694 borehole data by Zhu et al. (2016a).

Semivariogram of hydraulic conductivity
Semivariogram describes the degree of spatial dependence of a spatial random field or stochastic process. It is a concise and unbiased characterization of the spatial structure of regionalized variables, which is important in Kriging interpolations and conditional simulations. The experimental semivariogram can be fitted by an exponential model (e.g., Dai et al., 2014b) wherer k (h ϕ ) and r k (h ϕ ) are the experimental and model semivarograms of log conductivity Y for the kth facies at a lag distance h along the ϕ direction. In this paper we calculate the semivarograms in the vertical and dip directions. N(h) is the number of pair measuring points z o and z p separated by a h lag distance, σ 2 is the variance and λ the correlation range. The variance and range were optimized using the leastsquares criterion, which was solved by the modified Gauss-Newton-Levenberg-Marquardt method (Clifton and Neuman, 1982;Dai et al., 2012). The sensitivity equation method was derived to compute the Jacobian matrix for iteratively solving the gradient-based optimization problem (Samper and Neuman, 1986;Carrera and Neuman, 1986;Dai and Samper, 2004;Samper et al., 2006;Yang et al., 2014;Zhu et al., 2016a). The two sensitivity coefficients ∂r k ∂σ 2 and ∂r k ∂λ are the partial derivatives of the semivariogram with respect to variance and range:

Composite semivariogram of log conductivity
Once the facies semivariograms were obtained in each zone, the composite semivariogram γ (h) could be calculated through the following equation (e.g., Ritzi et al., 2004): where p k and t ki (h ϕ ) are the volumetric proportion of facies k and the transition probability from facies k to facies i in the ϕ direction with a h lag distance, respectively. Equation (8) delineates the composite semivarigoram with respect to the individual facies semivariogram and transition probability. The general shape function and range of the composite semivarigoram can be obtained from individual facies mean length and volumetric proportion with the methods described in Dai et al. (2005). The transition probability t ki (h ϕ ) has an analytical solution as derived by Dai et al. (2007): where δ ki is the Kronecker delta and λ ϕ is the integral scale in the direction of ϕ. A geostatistical modeling tool GEOST (Dai et al., 2014b) modified from the Geostatistical Software Library (Deutsch and Journel, 1992) and TPROGS (Carle and Fogg, 1997) were employed to compute the sample transition probabilities in each zone. The parameters p k and λ ϕ were optimally estimated through a modified Gauss-Newton-Levenberg-Marquardt method. More details are provided by Zhu et al. (2016a). The composite semivariograms for different zones can help us to understand the heterogeneity variations from the upper to lower part of the alluvial fan, as well as the stationary property (local versus regional) of the facies and hydraulic conductivity distributions.

Sequential Gaussian simulation
The sequential Gaussian simulation (SGSIM) is a widely used stochastic simulation method to create numerical model of continuous variables based on the Gaussian probability density function. The process is assumed to be a stationary and ergodic random process (Deutsch and Journel, 1992;Dimitrakopoulos and Luo, 2004). This method can preserve the variance and correlation range observed in spatial samples. SGSIM provides a standardized normal continuous distribution of the simulated variable. With the assumption that the log conductivity distributions are stationary within each zone, we used the SGSIM simulator implemented into GEOST to model the log 10 (K) continuous configuration under a multiple-zone framework. The conductivity of the FS, MS and G facies in each zone was simulated sequentially using the structure characteristics of the semivariograms.
Finally, the 3-D conductivity configuration was derived by combining the stochastic simulated facies (Zhu et al., 2016a) with the SGSIM conductivity distribution and the mean log 10 (K) of the various facies in each zone ( Table 2). The stochastic simulated facies was constructed through the optimized volumetric proportion and mean length of facies in three directions. The mean length in vertical and dip directions were calculated through 694 borehole. The mean length in strike direction was assumed as half as that in the dip direction. During the facies simulation process, borehole data were used as conditional data (Zhu et al., 2016a). In detail, since each cell is characterized by specific facies and zone indices, its conductivity was assigned using the corresponding (in relation to the facies and the zone) 3-D SGSIM outcome in that position. Note that the hydrofacies (e.g., C, FS, MS and G) are defined qualitatively based on the sedimentary structures, borehole lithological descriptions and grain sizes, while the conductivity samples are then deduced from geophysical measurements for each facies at each zone. Since the sub-clay and clay contents from zone 1 to zone 3 are increased due to the changes in the sediment transport conditions, for the same facies we also found this trend and the overall hydraulic conductivities are decreased from zone 1 to zone 3. Since sub-clay and clay are generally characterized by a low hydraulic conductivity value, a uniform K value equal to 0.0001 m day −1 was set to all the C cells.

Variation of log 10 (K) for the various facies
The optimized vertical correlation range and variance of the log conductivity semivariogram (Eq. 5) are listed in Table 3, along with their 95 % confidence intervals. The fitting between the experimental and the model semivariograms is the best in zone 2 because of the abundant samples, whereas the fitting in zone 1 is the worst (Fig. 6). The fitting result of the semivariogram for the G facies is the worst in zone 1. There are two reasons for this: the first is the high variance of the log conductivity of gravel in this zone, and the second is the limited number of samples (102 samples), which makes quite small the pair numbers within each lag spacing. Hence, the computed semivariogram is highly uncertain.
The variance of FS, MS and G in the vertical direction decreases from zone 1 to zone 3. In the upper alluvial fan, sediments were deposited under multiple water-flowing events and with poor sorting. The deposits consist of wide ranges of sediment categories and grain sizes. The variance of G is larger than 1.5, which reflects the high heterogeneity of hydraulic conductivity in coarse deposits. The variances of FS and MS are smaller with values equal to 0.23 and 0.32, respectively. In zone 3, these values decrease to 0.05 and 0.13, respectively, with that of G sharply decreasing to 0.62. In the middle-lower fan zone, the conductivity variation within each facies reduces gradually because the ground surface slope becomes smaller or flat, the sediment transport energy decreases and the deposits within the three facies are well sorted.
Note that the ranges are correlated with the facies structure parameters such as the indicator correlation scale, mean thickness (or length) and volumetric proportion (Dai et al., 2004b(Dai et al., , 2007. The estimated correlation ranges of FS, MS and G along the vertical direction in zone 1 do not show big difference with values equal to 6.0, 8.0 and 6.5 m, respectively. Zone 2 was extended from the fan apex zone (zone 1) with much larger area, which allows for greater preservation potential of finer sediments -such as MS, FS and C -than the more proximal zone 1. Therefore, in zone 2 the volumetric proportions for these three facies increase while that of gravel decreases. The estimated ranges of G and MS are increased. In zone 3, the range difference among the three facies decreases gradually. The range of FS is about 6.0 m, which is twice as much as that of MS. The spatial variation of the structure parameters of three facies causes the large changes of the correlation ranges from zone 1 to zone 3.
Due to the small number of conductivity samples in zone 1, the variance of log 10 (K) along the dip direction is calculated only in zone 2 and zone 3 (Table 4, Fig. 7), as observed along the vertical direction. This phenomenon is possibly a result of a sediment transport energy decrease along the flow direction. Lower energy flow in zone 3 allow for better sediment sorting and weak heterogeneity (or lower variance) in hydraulic conductivity.

Composited semivariogram of log 10 (K)
The composite semivariogram in the vertical direction at each zone is calculated by Eq. (8), using the volume proportions (Table 1) and transition probability (Eq. 9) with the same values of the lag distance used to compute the facies semivariograms (Fig. 8). The values of the optimized variance are 0.68, 0.11 and 0.03 in zone 1, zone 2 and zone 3, respectively. The high-flow energy and the large number of flooding events contributing to sediment deposition are the main causes of the high heterogeneity (largest variance) of the deposits in the upper part of the alluvial fan. The changes of variance between the three zones support the utilization of the local-stationary assumption and simulation of multiple-L. Zhu et al.: Modeling 3-D permeability distribution in alluvial fans using facies architecture Figure 6. Experimental (circle symbol) and model (solid line) semivariogram along the vertical direction for the various hydrofacies in the three zones. Notice that the range in the y axis differs for gravel lithology. zone-based conductivity distributions for the Chaobai alluvial fan.

Configuration of log 10 (K)
The configuration of log 10 (K) in three dimensions is showed in Fig. 9. The distribution of conductivity is generally consistent with that of the facies. Coarse units are more frequently distributed in the upper zone, which makes the average K much larger in this zone than that in the lower part of the alluvial fan. The regions with high conductivity (red color in Fig. 9) in zone 1 are more continuous than that in other parts. The adjacent cells with the smallest conductivity (blue color in Fig. 9) are obviously located mainly in zone 3. The mean conductivity is smaller in the southern part of the study area, where the piezometric drawdowns in the multiple-layer aquifer system were larger and the surface subsidence more serious (Zhu et al., 2013(Zhu et al., , 2015. Note that since we simulated the dip direction along the main water flow direction and, due to the lack of enough data, the strike-directional semivariogram is assumed to be similar to that in the dip direction; the simulated facies in the fan apex did not show a radiating   pattern. More information about simulating the radiating pattern can be found from Carle and Fogg (1997) and Fogg et al. (1998).
Based on the 3-D K configuration, the average value of K in the depth range from 0 to 300 m amounts to 194, 25 and 4 m day −1 in zone 1, zone 2 and zone 3, respectively. These values are comparable with those provided by the Beijing Institute of Hydrogeology and Engineering Geology (2007) based on a number of pumping tests carried out over several years in the study area. In this BIHEG report, the average value of K is > 300 m day −1 in zone 1, between 30 and 100 m day −1 in zone 2 and < 30 m day −1 in zone 3 (Fig. 1b). The fact that the arithmetic average K values are gently smaller than the latter ones is likely due to the fact that the outcome of pumping tests are generally more representative of coarser sediments.
Investigating the stochastic results along the vertical direction, it is interesting to notice that the average K in deep units of zone 1 and zone 2 is smaller than that in the shallow strata. For example, in zone 1 the average K for the cells from 0 to 100 m deep is 295 m day −1 , which is three times as much the value for the depth range between 200 and 300 m. Conversely, no significant variation of K versus depth is observed in zone 3, with only a small decrease of the average K from the deeper to the shallower units. the related variance of the hydraulic conductivity in 3-D domains. In particular, the optimized statistical parameters (e.g., log conductivity variance and correlation range) of semivariograms are estimated using the modified Gauss-Newton-Levenberg-Marquardt method. The Chaobai alluvial fan is used as a case study area. Multiple data, including downhole resistivity logging data, vertical electric soundings, well-bore lithologic logs, TDS measurements and depths to the water table, are integrated to derive a dataset of conductivity values in a 3-D setting. Log conductivity semivariograms fitted with exponential functions were constructed for three facies, including fine sand, medium-coarse sand and gravel, in each of the three zones into which the Chaobai fan is divided to guarantee local stationarity of the statistical process. The composite semivariogram of the three facies has been derived for the two zones where a sufficiently large number of samples are available. The log 10 (K) configuration is simulated using the sequential Gaussian simulation model based on statistic parameters of log 10 (K) and the structure suggested by a 3-D hydrofacies simulation.
For the specific test case, the variance along the vertical direction of fine sand, medium-coarse sand, and gravel decreases from the upper part of the alluvial fan, where the values amount to 0.23, 0.32 and 1.60, to the lower portion of the Chaobai plan with values of 0.05, 0.126 and 0.62, respectively. This behavior reflects the higher transport energy in the upper alluvial fan that causes a poor sediment sorting.
In the middle alluvial fan, the transport energy decreases and the sediments tend to be relatively well sorted. The variance of the gravel is larger than that of other lithologies. The different flow energy significantly affected the coarse sediments in the vertical direction. Along the dip direction, the variance of three facies (gravel, medium-coarse sand and fine sand) in the middle fan is larger than that in the lower fan. The composite variance of log 10 (K) in the vertical direction shows that the large heterogeneity in the upper fan (with a value of 0.68) decreases in the lower zone.
The distribution of hydraulic conductivity is consistent with that of the facies. Hydraulic conductivity is much larger in the upper zone than that in the lower part of the alluvial fan. This result provides valuable insights for understanding the spatial variations of hydraulic conductivity and settingup groundwater flow, transport and land subsidence models in alluvial fans.
Concluding, it is worth highlighting that we depicted an original method to detect the variance and configuration of conductivity by fusing multiple-source data in 3-D domains. The proposed approach can be easily used to statistically characterize the hydraulic conductivity of the various alluvial fans, which worldwide are strongly developed to provide high-quality water resources. We are aware of some restrictions in the dataset available at the date for the Chaobai alluvial fan, for example the assumed uniform distribution of TDS versus depth and the relatively small number of the conductivity samples in the upper fan zone. A more accurate description of the semivarigrams in the dip and lateral directions will be included in our future study to improve the developed 3-D permeability field. Moreover, our assumption that the logs are well calibrated might be another source of uncertainty that can be reduced in our forthcoming work. Nonetheless, the proposed methodology will be re-applied in the near feature as soon as new information will become available, thus allowing one to improve the estimation accuracy of spatial statistics parameters and the configuration of hydraulic conductivity in this Quaternary system so important for the Beijing water supply.

Data availability
The geophysical measurements, borehole lithostratigraphies and hydrogeological information in the northern part of the Beijing Plain can be partly accessible by contacting Beijing Institute of Hydrogeology and Engineering Geology.
Author contributions. Lin Zhu, Huili Gong and Zhenxue Dai derived the method of spatial variance and 3-D configuration of conductivity, performed data analysis and wrote the draft manuscript. Gaoxuan Guo collected the geological and geophysical data, and discussed the results. Pietro Teatini discussed the results, and reviewed and revised the manuscript.
Competing interests. The authors declare that they have no conflict of interest.