Vertical trend analysis of equivalent hydraulic conductivity in alluvial fan gravel deposits considering open void connectivity

: This study demonstrated groundwater flow simulations to investigate a vertical trend of equivalent hydraulic conductivity of alluvial fan gravel deposits in Sapporo, Japan, considering open void connectivity. Equivalent hydraulic conductivity was defined according to Darcy’s Law for a cube of 10 m in size, consisting of one million cells assigned among fully packed (without open voids), loosely packed (with less-connected voids) or very loosely packed (with well-connected voids) deposits. The stochastic generation was performed under each configuration in terms of target depth sections for vertical trend analysis, and horizontal variogram ranges (random, high, and low connectivity, and no open voids) for open void connectivity. The logarithmic average of 100 equivalent hydraulic conductivities was calculated in each configuration, and the vertical trends were determined. The simulation results showed that the equivalent hydraulic conductivity increased when the open void frequency was large in the shallow zone and the connectivity of the open voids was assumed. In particular, the high connectivity assumption was needed to match the in situ trend with a decay exponent of 0.05 m –1 . Modeling the vertical trend with such a large decay exponent was essen-tial to obtain realistic solutions of the groundwater flow and transport system in the alluvial fan.


INTRODUCTION
Information about vertical trends in hydraulic conductivity (K) can contribute to an improved understanding of groundwater flow and transport systems over a range of scales (Saar and Manga, 2004;Jiang et al., 2009;Cardenas and Jiang, 2010;Zlotnik et al., 2010). K of geological materials in a vertical direction generally varies as a result of reductions in porosity through physical and chemical processes (Ingebritsen et al., 2006). Of these, the compaction under buried pressures has an important influence on K especially in unconsolidated deposits. Jiang et al. (2010) found that K decreased exponentially with depth (Z) as the pressure increased in both unconsolidated sediments and Correspondence to: Yoshitaka Sakata, Division of Human Environmental Systems, Faculty of Engineering, Hokkaido University, Kita-13 Nishi-8 Kita-ku, Sapporo, Hokkaido 060-8628, Japan. E-mail: y-sakata@eng. hokudai.ac.jp consolidated rocks, as follows: where K 0 is K at the ground surface and A is the decay exponent, measured in m -1 . Similarly, from permeability tests of consolidated rocks in previous studies, Jiang et al. (2009) reported the range of A was between 10 -3 and 10 -4 m -1 .
Despite the various studies about consolidated rocks, few have examined the vertical trends of K in unconsolidated sediments. One practical reason is that sedimentary formations have traditionally been conceptualized as geologically homogeneous, horizontally stratified layers, and it is generally assumed that the layers are not sufficiently thick to contain vertical trends. However, this assumption may not be valid for some topographic types, such as alluvial fans where unconsolidated gravel deposit layers may range from tens of meters to > 100 m thick. In addition, the monotonic gravel sequences cannot be clearly defined especially in the deep zone from limited and disturbed core samples.
Hydraulic conductivity K of unconsolidated sediments is conceptualized to be a function of effective grain size related to total porosity (Bear, 1988). In a sedimentary basin, effective grain size is variable not only in the horizontal direction, for example coarsening toward the upstream direction, but also in the vertical direction due to changes in the depositional environment. Also, in an alluvial fan the grain size often increases as the sampling depth decreases because of progradation or the fining trends toward the surface caused by recession (Einsele, 2000). Thus, vertical changes in grain size will influence vertical trends of K in alluvial fans (Neton et al., 1994). Furthermore, K of unconsolidated gravel deposits is closely related to the open voids between gravel grains, where preferential flow paths form as fractures in consolidated rocks. Sedimentary studies found the gravel deposits with open voids in trenches and outcrops, and K of these deposits was several orders of magnitude higher than in other deposits (Jussel et al., 1994;Heinz et al., 2003;Zappa et al., 2006;Ferreira et al., 2010). The field evidence of open voids in unconsolidated gravel deposits was also inspected as loose packing of fine materials in undisturbed core samples. For example, in the Toyohira River alluvial fan, Sapporo, Japan, Sakata and Ikeda (2013) classified the undisturbed cores into three open void indicators and indicated that the frequency of loose packing parts was strongly related to high values of slug test results at the sampling depths. The previous study also determined the vertical trend of core-scale K with the decay exponent of A = 0.11 m -1 . Such a large decay exponent plays a role of controlling the groundwater flow and transport near the ground in the alluvial fan (Sakata, 2014).
Vertical trends of such core-scale K may not, however, be directly applicable to regional studies of which unit scale (mesh size) was larger, for example over ten meters, than the borehole diameter for sampling and testing. The K assigned to each mesh in regional models should correspond to the equivalent value as a result of the spatial distribution of small-scale K values that are not random, but are correlated (Carrera and Mathias, 2010). Especially in alluvial fan gravel deposits, the connectivity of open voids strongly influences the equivalent K, because core-scale K with open voids is much higher than without. The connectivity can be generated in an assumption of horizontal and vertical correlations as variograms (or covariances) in geostatistical modeling (dell'Arciprete et al., 2012). The vertical variograms can be determined even from core samples in a single borehole. On the other hand, the determination of the horizontal variograms is almost impossible in practice with limited data.
In this study, the vertical trends of equivalent K by generating open void connectivity are examined in the alluvial fan gravel deposits in Sapporo, Japan. The equivalent K is calculated by using numerical flow simulation in a cube that consists of one million cells, assigned to one of three open void indicators. A stochastic approach is used to realize the distributions of open voids in the cube using assumed configuration of seven depth sections and four scenarios of open void connectivity: no open voids, and random, low, and high connectivity. The averages of the equivalent K of the cube are calculated for seven depth sections and its vertical trend is determined from the averages for each scenario. The simulated trends for different scenarios are compared with those derived from in situ measurements, and the significance of considering open void connectivity for groundwater modeling in the alluvial fan gravel deposits is discussed.

SITE DESCRIPTION
The study site is located on the alluvial fan of the Toyohira River, at approximately 43°N and 141°E on Hokkaido, the most northerly island of Japan. Figure 1 shows the regional and the local topographic maps of this study, and locations of wells. Sapporo, the capital city of Hokkaido with a population of about 1.9 million, is located on the alluvial fan. The subsurface geology of the fan was described in previous studies (e.g. Oka, 2005). This study briefly described the summary as below. The fan formed beside the Toyohira River as a result of deposition during the Quaternary. The fan has a radius of ~ 7 km and extends over an area of ~ 31 km 2 . The river flows northwards from mountains comprising Tertiary volcanic rocks, through the fan, to the lowlands near Ishikari Bay. The fan is divided into a western Holocene section and an eastern Pleistocene section (Daimaru, 1989). Originally, Sapporo City was mainly confined to the Holocene section because the water table was closer to the surface and the groundwater quality was better than on the Pleistocene section. Before the 1970s, most wells in the Holocene fan deposits were drilled to approximately the same depth as the water table to facilitate water extraction. Since then, the water table has been depleted by ~ 4 × 10 4 m 3 per year because of on-going groundwater extraction and underground constructions such as subways. Hu et al. (2010) also indicated that the groundwater resources in the alluvial fan are presently under threat, including problems relating to climate change, as in other alluvial fans of Japan.

IN SITU MEASUREMENTS FOR COMPARISON
Pumping tests are generally performed in water wells to obtain the prespecified discharge rates and the resultant changes in water level. The hydraulic properties (transmissivities and storage coefficients) can be estimated from the water level variations during the tests. In the alluvial fan, however, available data of previous pumping tests are limited in public, and most of them provided only the measurements, without data analysis of hydraulic properties. Among them, this study accumulated pumping test data from 13 wells with analyzed transmissivity values of the alluvial fan gravel deposits. The data were obtained from the Hokkaido Regional Development Bureau, which is responsible for managing the Toyohira River. Although the number was limited, the test data cover almost the entire Holocene fan (Figure 1). The well screens of the test boreholes were at different depths, and the transmissivities were not constant because K was spatially variable. The pumping test results and the average hydraulic conductivities (K p ) over each screen, along with the center depths of the well screens (Z p ), are summarized in Table SI, and the plots are presented in Figure 2. The center depths of the well screens ranged from 5.05 to 52.50 m. The common logarithm values, log K p , which ranged from -4.18 to -2.51 log m s -1 and varied by more than one order of magnitude, indicated that the alluvial fan gravel deposits were generally highly permeable. As shown in Figure 2, K p decreased logarithmically with depth, and so the parameters of Equation (1) were determined by least squares regression; K 0 and A are 10 -2.51 m s -1 and -0.052 m -1 , respectively, with a coefficient of determination (R 2 ) of 0.73. The value of A was one or two orders of magnitude greater in the gravel deposits than in consolidated rocks, as described above (Jiang et al., 2009). The vertical trend was therefore important for the groundwater flow and transport systems in the fan relative to consolidated formations. However, the in situ value of A was almost half of the value (0.11 m -1 ) in the previous core analysis (Sakata and Ikeda, 2013) at the mid-fan (near the pumping test locations 1 and 4 in Figure 1). The earlier study assumed the difference of A was related to the open void connectivity in the gravel deposits, indicating typically high values (over 10 -1 m s -1 ) of core-scale K near the ground, probably due to the local open voids at the sampling depths. Thus, the high A (= 0.11 m -1 ) value of corescale K was influenced strongly by the several high values of K. On the other hand, the typically high values of K p were not measured, because the open void connectivity was limited in the larger volumes for testing.
In Figure 2, the colors of the plots indicate the horizontal distances (km) from the fan apex (the southern edge of the fan); 43.015°N, 141.350°E. In most sedimentary basins, there are fining of grain sizes in the downstream direction of the river, resulting in the horizontal trend of K (Freeze and Cherry, 1979). However, Figure 2 indicates that the horizontal distances were not obviously related to K p as

MODELING AND ANALYSIS METHOD
This study demonstrated to calculate equivalent K (K e ) of the gravel deposits including open voids through the numerical simulation. For this purpose, this study defined the cube of 10 × 10 × 10 m in size, consisting of one million cells, each of which measures 0.1 × 0.1 × 0.1 m. K e of each cube was calculated according to Darcy's Law, as in Desbarats (1987), by dividing simulated discharges by the cross-section of the model and the hydraulic gradient, which was assigned as a boundary condition. Each cell was assigned stochastically to one of three open void indicators, as defined in our earlier study (Sakata and Ikeda, 2013), namely high packing (HP), loose packing (LP), or very loose packing (VLP). The previous study (Sakata and Ikeda, 2013) also derived core-scale K of each open void indicator as a function of the effective grain size for HP in porous media, or of length fractions (frequencies) for LP and VLP in fractured media. Photographs, explanations, and the permeability of the cores from the earlier study are summarized in Table I. The coefficients of core-scale K models in Table I were determined to match the slug test results at the same sampling depths. In the previous study, the effective grain size was the size of 20% cumulative weights by the sieve tests for undisturbed cores. The effective grain size showed the upward fining trend of depth, indicating the vertical trend of K according to the porosity change. The length fractions of LP and VLP in each depth was measured by the inspection of the undisturbed cores. The previous study assumed the LP and VLP as gravel deposits with open voids because the core-scale K of LP and VLP was higher than K of HP, especially when the length fractions increased. The resulted proportionality of VLP, at 1.87, was much higher than that of LP, at 0.0167, indicating that VLP contributed to preferential groundwater flow paths in the gravel deposits. This was probably because the open voids were well-connected in VLP, relative to LP. In other words, K e were dependent primarily on the frequency and the connectivity of VLP.
The HP, LP, and VLP cells were stochastically generated for different configurations in terms of depth sections for vertical trend analysis, and of horizontal variogram ranges of LP and VLP for open void connectivity. This study defined seven depth sections for analysis, as follows: 1) Z 1 = 0-10 m, 2) Z 2 = 5-15 m, 3) Z 3 = 10-20 m, 4) Z 4 = 20-30 m, 5) Z 5 = 30-40 m, 6) Z 6 = 40-50 m, and 7) Z 7 = 50-60 m. The generation was undertaken using the algorithm for sequential Gaussian simulation (Deutsch and Journel, 1998). All configurations, i.e. effective grain sizes (d e ), frequencies (p) of HP, LP and VLP, and vertical variogram ranges (λ v ) of HP, LP, and VLP for a spherical model were determined from the core analysis data of the previous study (Sakata and Ikeda, 2013). The configurations were summarized in Table SII. The core data were obtained at two boreholes near the pumping test locations 1 and 4. The parameters, d e and p values, were averages of the sieve test VERTICAL TREND OF HYDRAULIC CONDUCTIVITY data and the length fractions of LP and VLP, respectively, in each depth sections. Another parameter, λ v values in the spherical variograms of each open void indicator, was determined through indicator transformation (Deutsch and Journel, 1998) to match experimental variograms of core data with the theoretical spherical model in the least square. In order to address the uncertainty of horizontal variogram ranges in limited borehole data, this study used a practical approach to assume horizontal to vertical anisotropy of variogram range λ h /λ v , as proposed in Deutsch (2002). This study also assumed that the braided river condition in Deutsch (2002) was the same as those in the alluvial fan, i.e. the value of λ h /λ v ranged from 20 to 100. Therefore, four scenarios of open void connectivity were used in this study, namely no open voids (p of HP = 1, p of HP and LP = 0), random connectivity (λ h = 0), low connectivity (λ h /λ v = 20), and high connectivity (λ h /λ v = 100). Three examples of the stochastic models are presented in Figure 3. In Figure 3a from the shallowest layer, Z 1 and λ h /λ v = 100, the LP and VLP cells were connected horizontally throughout the cube, and the connected zones acted as preferential flow paths in the groundwater flow simulation. In Figure 3b from a deeper layer, Z 4 , the LP and VLP cells were connect-ed locally but the connected zones did not extend throughout the whole cube. In Figure 3c from a deeper layer with a random connectivity assumption, the densities of LP and VLP cells were high throughout the cube, but not as well connected.
The K of each cell as HP, LP or VLP was assigned individually by a function, listed in Table I. K of HP generally decreased with depth from Z 1 to Z 7 according to the fining d e , but was assumed to be constant within each depth section. The L values of LP and VLP in the functions corresponds to the vertical length of open voids as the width of preferential paths. This study counted the number of individual 0.1 m cells that were vertically connected in the stochastically generated fields. In this simulation, constant head boundaries of 20 and 19 m were assigned to the upstream and downstream sides of the cube, respectively, resulting in the mean hydraulic gradient of 0.1 (m/m), as the highest value of the water table around the losing section of the river discharge (Sakata, 2014). However, the calculation results of K e were rarely dependent on the head boundary condition (the comparison was omitted due to paper length restriction). In each cube, the steady discharge of water was calculated by simulating the flow among the Most pore spaces between gravel grains are fully packed by fine sediments. A lack of fine sediments is partially evident, but the apertures are limited within the sample diameter.  (Harbaugh, 2005). K e was calculated by dividing the simulated discharges by the hydraulic gradient (0.1 m/m) and the cross-sectional area (100 m 2 ). This study applied the stochastic approach to realize 100 spatial patterns of K for each configuration of seven depth sections and four open void connectivity assumptions. The author considered that the approach was more effective to address the uncertainty of open void connectivity than another conventional approach to specify a spatial pattern of K by matching calculated hydraulic heads with measured heads. The logarithmic average of K e from 100 iterations at each depth (K − e ) was calculated and confirmed to evaluate its convergence with the small standard deviation (σ) of log K e . After the confirmation, the vertical trends of K − e were determined by the regression method in Equation (1)

RESULTS AND DISCUSSION
The simulation results (i.e. log K − e and standard deviation σ of log K e ) for the seven depth configurations and the four open-void connectivity scenarios were determined as summarized in Table SIII. were the largest, from -2.49 to -3.70 m s -1 in each depth section, when the high connectivity was assumed. These comparisons indicate that K − e was dependent on the effective grain size, the open-void frequencies and the open void connectivity. Resulted σ of log K e was less than 0.1, and its relative ratio to log K − e was less than 1% in most configurations except the high connectivity assumption at Z 1 (4%) and Z 3 (2%). The small values of the relative ratio indicated that the K − e values almost converged in each configuration by calculating 100 stochastic iterations.
The vertical trend of K − e were determined as the regression lines in Figure 4. The R 2 values ranged from 0.66 to 0.72 and increased when the open-void connectivity was assumed, particularly when the connectivity was high. Similar to K − e , the K 0 values also increased when the open voids were generated and their connectivities were assumed, indicating that K − e was more variable at shallower depths according to the occurrence of open voids and the connectivity, than at deeper depths where there were limited open voids. As a result, K 0 and A reached maximum values of 10 -2.47 m s -1 and -0.051 m -1 , respectively, when the high open-void connectivity was assumed. The simulated trends were similar to the in situ measurements (K 0 = 10 -2.51 m s -1 and A = -0.052 m -1 ), which showed that the parameters of the in situ and simulated trends were almost equal only when the high open-void connectivity was assumed. Even though in situ data for comparison were limited and the actual connectivity of the open voids also remains uncertain, the agreement suggested that the simulation results were reasonable and that the vertical trends were controlled by the open-void connectivity, as well as the grain size and the open void occurrence.
For the simulated and in situ trends, A was approximately -0.05 m -1 and was one order of magnitude greater than in the consolidated rocks. In most sedimentary basins, vertical trends with such high decay exponent are key for interpreting groundwater flow and transport on a regional scale. In the Toyohira River alluvial fan, for example, the groundwater is mostly recharged by the precipitation from the surface and by the seepage from the riverbed (Mukai et al., 2008). When the decay exponent is very large, groundwater infiltration to the deep zones is inhibited and the groundwater flows mainly in horizontal directions within shallow high-permeability aquifers, as shown in the numerical flow simulation (Sakata et al., 2016). The vertical trend is probably more significant for the transport system into the deep zone, for example of infiltration of dense contaminants such as heavy oil in the developed area of the midfan. A future study will consider how the vertical trend in the permeability of the alluvial fan gravel deposits influence groundwater transport.

CONCLUSIONS
In this study, vertical trends in equivalent hydraulic conductivity of the alluvial gravel deposits, Sapporo, Japan, for a cube of 10 × 10 × 10 m in size at different depths under four assumptions of open-void connectivity were evaluated from numerical flow simulations. The open void connectivity was generated through the stochastic simulation of three open void indicators. The simulations showed that equiva-VERTICAL TREND OF HYDRAULIC CONDUCTIVITY lent hydraulic conductivity increased when the open voids frequency and the open void connectivity were assumed, especially under the high connectivity assumption. The equivalent hydraulic conductivity also decreased with depth as the open voids were rarely seen in deeper zones. This simulation study revealed that the simulated vertical trend was also dependent on the connectivity between open voids generated in the cube. Of the four scenarios, the trends for the high connectivity between open voids were similar to the in situ trends. The agreement indicated that the high connectivity between open voids were needed to be considered when analyzing vertical trends in the equivalent hydraulic conductivity of the alluvial fan gravel deposits. The decay exponent was 0.05 m -1 , at least 10 times higher than for consolidated rocks, indicating that the groundwater flow and transport systems are controlled by the vertical trends of equivalent hydraulic conductivity in the alluvial fan.

ACKNOWLEDGMENTS
The author sincerely thanks the Hokkaido Development Bureau, Japan, for granting permission to use the core data. The author is also grateful to our colleagues, Prof. Ryuji Ikeda and Prof. Kazuhisa A. Chikita, for their invaluable advice and guidance. The analysis work was supported by our graduate student Mr. Maeda Shingo with dedicated efforts. The manuscript was improved by constructive and suggestive comments from two anonymous reviewers.