Published in Vadose Zone Journal 3:786-795 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: RESEARCH ADVANCES IN VADOSE ZONE HYDROLOGY THROUGH SIMULATIONS WITH THE TOUGH CODES
Role of Rock Heterogeneity on Lateral Diversion of Water Flow at the Soil–Rock Interface
Niclas Bockgård* and
Auli Niemi
Dep. of Earth Sciences, Uppsala Univ., Villav. 16, SE-752 36 Uppsala, Sweden
* Corresponding author (niclas.bockgard{at}hyd.uu.se)
Received 22 August 2003.
 |
ABSTRACT
|
|---|
Groundwater recharge is an important, yet often a highly uncertain upper boundary condition in groundwater flow models for fractured rock. The net infiltration at ground surface can often be estimated from water budgets, but the redistribution of water between a soil layer and the underlying bedrock is much more difficult to quantify. To address this question, we studied the role of lateral flow diversion at the soil–rock interface and how it influences the percentage of net infiltration that becomes inflow into the bedrock. Numerical experiments were performed with different net infiltration conditions, soil–rock hydraulic conductivity contrasts, hydraulic gradients, and most notably, rock heterogeneity. The rock was first represented as a deterministic homogeneous continuum, then as a stochastic heterogeneous continuum and analyzed with Monte Carlo simulations. The heterogeneous simulations predicted a large variation in the inflow into the rock between different stochastic realizations (from 50 to 100% of the net infiltration). Furthermore, at high hydraulic gradients in the bedrock, the mean flux into the rock from the heterogeneous simulations was smaller than the corresponding results from the homogeneous medium. The heterogeneous simulations showed formation of local unsaturated zones below the groundwater table. These occurred in cases with large hydraulic gradients at locations with large hydraulic conductivity contrasts, with low-conductive regions overlying high-conductive ones. This phenomenon was not captured with the homogenous models, which thereby may overestimate flux into fractured rock in cases where hydraulic gradients are large.
|
INTRODUCTION
|
|---|
GROUNDWATER RECHARGE is an important factor in the evaluation of groundwater resources, such as assessing the impact of underground construction work or water withdrawal activities. Recharge is defined as the flux of water through the water table into the saturated groundwater system. It forms the upper boundary condition for groundwater flow and transport models, making the results from these models very sensitive to applied recharge. A key issue is that groundwater recharge is difficult to quantify. Recent reviews include de Vries and Simmers (2002) and Alley et al. (2002). Net infiltration, defined as the difference between infiltration and evapotranspiration, can often be estimated from an overall water budget, at least as a large-scale temporal mean value. However, the redistribution of water at depth, such as between a soil layer and the underlying bedrock, is much more difficult to quantify. Three different situations are possible:- The soil layer is unsaturated and the water table is located in the bedrock.
- The water table is located in the soil and the bedrock is saturated.
- Perched groundwater in the soil is overlying an unsaturated zone in the upper bedrock.
In the two latter situations, often only a small part of the net infiltration will eventually enter the underlying bedrock. In terms of the water budget of an overlying high-conductive soil layer, a small flux from the soil into the underlying low-conductive bedrock may be insignificant in relation to other fluxes and is therefore often neglected. In terms of the water balance and flow processes of the underlying low-conductive, low-storativity bedrock, however, the value of this upper boundary condition and related uncertainties are significant.
The present study investigates lateral diversion of water flow at the interface between a soil layer and underlying fractured bedrock and how it influences the percentage of net infiltration that becomes inflow into the bedrock. The objective is to gain understanding and insight into (i) to which extent the infiltration at the land surface reaches the bedrock in various hydrogeological situations and (ii) how this is influenced by fracture-induced heterogeneity of the bedrock. Numerical simulations are performed in a typical soil–rock profile consisting of a soil layer overlying fractured bedrock. First, a deterministic homogeneous representation is used for the rock to gain insight into the mean behavior. Second, a stochastic heterogeneous description with Monte Carlo analysis is used to investigate how the processes and mechanisms of interest are influenced by the fracture-induced heterogeneity.
We are aware of only a limited number of earlier studies addressing the distribution of water between a soil layer and underlying bedrock. Harte and Winter (1995) simulated the effects of hydraulic conductivity contrast and topography on recharge to a bedrock aquifer for a hypothetical profile of glacial deposits and fractured crystalline bedrock in a humid climate. They concluded that the hydraulic conductivity of overlying soil deposits is important in controlling recharge rate to the bedrock, while topography and large-scale trends in hydraulic conductivity of the rock affect the size of the bedrock recharge area. For the same geological setting, Tiedeman et al. (1998) simulated the distribution of groundwater flow between soil and bedrock on a basin scale. They found that the distribution of groundwater flow is strongly controlled by the transmissivity of the soil and bedrock and, as can be expected, by the position of the water table. Large-scale flow is simulated in these studies, and the effect of fracture-induced, small-scale heterogeneity is not taken into account.
An intensively instrumented field site designed for studying groundwater recharge processes in both soil and the underlying fractured bedrock was described by Gburek and Folmar (1999). The depth to bedrock is less than 0.5 m and the water table is typically 6 m below land surface. Monitoring includes, for example, percolation through lysimeters in the uppermost 1 to 2 m of the bedrock and measurement of groundwater levels. It was observed that groundwater levels always responded 1 to 2 h after the start of percolation. Furthermore, the responses of different lysimeters to precipitation were very different, which raises the question of the role of heterogeneity on the recharge process. Lee and Lee (2000) used time series analysis of groundwater level measurements in soil and rock in combination with rainfall data to characterize the recharge patterns of a bedrock aquifer. They related the different responses of a shallow well in soil and a deep well in rock to the recharge to the soil zone and the flow through the fracture network, respectively. However, neither the quoted modeling studies, nor the analysis by Lee and Lee (2000) address the role of the small-scale heterogeneity at the soil–rock interface, the topic of our research.
It should be pointed out that in our study the water input at land surface is given as a constant net infiltration. Thus, the possible effects that subsurface processes can have on evapotranspiration are not taken into account. This is a justified assumption given the humid climate conditions considered here, at least in a preliminary analysis like this. In humid areas, the evapotranspiration, and hence the net infiltration, is often controlled by meteorological factors, or when it is limited by soil moisture, it is the topmost soil layer, the root zone, that is involved. The diversion at the soil–rock interface is therefore likely to have only a small influence on evapotranspiration and thereby the net infiltration. In arid areas, however, the soil water storage would have a large effect on evapotranspiration, and several mechanisms exist that can draw water from large depths. Capillary rise can be significant in relation to percolation, and trees and shrubs are reported to extract water from depths of more than 50 m in loose deposits. Vapor transport also forms a considerable flux at sites with a large vertical temperature gradient, for example, due to seasonal temperature variations in deserts. In general, groundwater recharge in arid areas is much more sensitive to the near-surface conditions than in more humid areas (de Vries and Simmers, 2002) and the recharge depends greatly on the ability of water to escape evapotranspiration by rapid percolation. De Vries and Simmers (2002) provided a detailed description of the climatological, geological, and biological processes affecting groundwater recharge.
|
MATERIALS AND METHODS
|
|---|
Conceptual Model
The model domain is a two-dimensional profile of a soil layer overlying fractured bedrock (Fig. 1)
. It represents a small part of a large-scale flow system, more specifically the very upper part of a ridge. The ridge is a recharge area, with essentially vertical inflow into the subsurface, recharging the large-scale flow system. Groundwater in the overlying soil layer is partly discharging also locally, to a nearby stream on the hill slope. We model the flow through the unsaturated and saturated zones of this soil–bedrock system and vary the net infiltration rate, the soil–rock hydraulic conductivity contrast, the hydraulic gradient, and the heterogeneity of the rock material. We allow the model to find the proper position of the water table and then observe how much of the net infiltration that enters the bedrock in various situations. Depending on the particular scenario, the water table can thus be situated either in the soil or in the underlying bedrock.
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 1. The model domain and types of boundary conditions used in the simulations. The slope of the land surface and soil–rock interface is 5%. The model domain is a small part of a large-scale flow system (right).
|
|
The boundary conditions are assigned such that the flow in the bedrock is predominantly vertical, representing a groundwater recharge zone in the bedrock, while in soil water is also allowed to flow laterally, discharging locally to a nearby surface water stream. To represent this situation, no-flow boundary conditions are applied to the left boundary of the domain (symmetry boundary) as well as to the right boundary of the rock zone (Fig. 1). A specified head boundary condition is used at the lower boundary, the value of which is varied to study the effect of different hydraulic gradients in the rock that could arise from the large-scale topography or from different underground constructions, most typically from a tunnel. A small gradient (i.e., a large positive head on the lower boundary) represents a relatively flat topography, while a large gradient represents a steep topography or possibly a tunnel in the rock. A tunnel located directly below the model would result in free drainage, which is equivalent to an atmospheric pressure condition at the lower boundary of the model. The reference datum, z = 0, is the bottom boundary of the model. The head on the boundary is varied between h = 0 m (e.g., corresponding to a tunnel opening at the bottom of the model) to h = 10 m. Lateral drainage of the soil zone is simulated using a head-dependent flux boundary condition according to
 | [1] |
where qx (m s–1) is the lateral flux over the boundary, h(z) is the hydraulic head at the boundary, h0 is a reference head, and C (s–1) is saturation-dependent hydraulic conductance. The reference head is fixed to a level 2.5 m below the level of the soil–rock interface at the boundary, and the hydraulic conductance is chosen to correspond to the conductance of a 50-m-long column of soil material. This boundary condition could be seen as drainage to a nearby stream with constant head, located downstream from the boundary. The upper boundary condition is a spatially uniform and constant net infiltration. Rates of 200 and 150 mm year–1, typical for conditions in Sweden, were used. Since net infiltration was used, evaporation and transpiration were not simulated.
We do realize that the results are highly sensitive to the boundary conditions imposed. They were chosen to represent typical conditions in Sweden in a bedrock recharge zone, but more scenarios could obviously be included. The slope of the soil–rock interface and the head value controlling the head-dependent flux for horizontal drainage in the soil could be systematically varied as well, like the lower boundary condition. To limit the number of possible cases, however, this was not done in this study. Since the primary objective of our study was to look at the significance of rock heterogeneity on the amount water inflow that enters the bedrock, we considered this a justified choice because the effect of boundary conditions will be similar in the homogenous and heterogeneous cases.
Medium Properties
The hydraulic properties were selected to represent those of a medium-textured soil and typical fractured crystalline rock as encountered in Sweden. The properties of the soil correspond to those of a sandy-loamy till at the NOPEX Site north of Uppsala, Sweden (Lundin et al., 1999). The soil material was modeled as a homogenous porous medium. The properties are summarized in Table 1.
For the rock, properties typical for fractured crystalline rock (e.g., Niemi et al., 2000, Öhman and Niemi, 2003) were used (Table 1). The fractured rock was first described as a homogeneous medium with deterministic properties. Next, the effect of fracture-induced heterogeneity was considered by modeling the rock as a heterogeneous stochastic continuum. In the latter case the saturated hydraulic conductivity, Ks, was assumed to be lognormally distributed and spatially correlated. A 0.5-m support scale was selected for the stochastic analyses. We wish to clarify that the support scale (Neuman, 1987) is the scale in which the heterogeneity characteristics are introduced into a stochastic model; in other words, it is (i) equal to the size of the hydraulic conductivity zones of the model, often also the same as the size of the numerical elements, and (ii) it should also be equal to the scale in which the field measurements underlying the input statistics have been measured (such as packer interval of hydraulic tests). When using a certain support scale in a stochastic continuum model, the underlying assumption is also that the data at this scale fulfills the criterion of continuum behavior (for more details see e.g., Niemi et al., 2000). The 0.5-m support scale was used here because it is (i) a reasonable scale where one can assume the continuum approximation to be valid for the conductivity characteristics and (ii) we have data available, both concerning the saturated and unsaturated characteristics.
On the basis of data from other fractured rock sites (Fig. 2)
, a standard deviation of log10(Ks) of 1.5 and a correlation length of 2 m can be considered reasonable for the 0.5-m support scale used. Inspection of the standard deviation of log(K) as a function of the scale of measurement in Fig. 2a appears to show no clear dependency. For small-scale data (0.1–3 m) standard deviation values between 0.5 and 2.5 are observed and the value 1.5 was chosen from the middle of this range. The correlation length vs. scale data (Fig. 2b) shows scatter as well, but most small-scale observations fall in a range that justifies choosing a 2-m correlation length for the 0.5-m support scale. An exponential isotropic semivariogram model with zero nugget and a range of 2 m were used in the modeling.
Models for Unsaturated Fractured Rock
While models for unsaturated flow in porous soils are relatively well established, those for unsaturated rocks are far less so. Good recent reviews for the latter include Wanfang et al. (1997) and Pruess et al. (1999a). The available physically based models can be divided into discrete fracture models and continuum models. Discrete fracture models are usually best suited for small-scale analyses and for upscaling the local, small-scale information, while continuum models are used in larger scales. A special case of continuum models are dual-continuum models that treat matrix and fractures as two separate but interacting continua. Unsaturated flow in fractured rock can be dominated by the characteristics of the rock matrix, by the characteristics of the fractures, or a combination of the two, in which the matrix characteristics dominate the flow in some saturation ranges and fracture characteristics dominate in others.
For our analysis we use the so-called fracture continuum approach (Birkholzer et al., 1999, Finsterle, 2000). In doing so we assume (i) flow is dominated by the fractures and (ii) the rock behaves as a continuum at the 0.5-m support scale of our stochastic analysis. For the continuum assumption to be valid, the rock has to be relatively densely fractured, an assumption that is reasonable since we are considering the very uppermost part of the bedrock. The fact that the flow is dominated by the fractures, allowing one to neglect the matrix contribution, is quite well established for this type of rock in saturated conditions (e.g., Niemi et al., 2000, Öhman and Niemi, 2003). In terms of the unsaturated conditions there are no data available from Scandinavian crystalline rocks, but we assume the models developed for the unsaturated zone of the relatively densely fractured Yucca Mountain to be applicable. We therefore adopt the fracture continuum approach as well as the appropriate parameter values used by Birkholzer et al. (1999) in their fracture continuum model for the fractured tuff at Yucca Mountain (Table 2). Table 2 gives examples of some of the few data that presently exist concerning parameters for unsaturated fracture flow. The data in Table 2 come from the Yucca Mountain and Box Canyon sites in the USA and from the Grimsel site in Switzerland, some of the few sites where researchers have performed detailed investigations concerning the unsaturated flow characteristics for fractured media.
To describe the capillary pressure–liquid saturation and hydraulic conductivity–liquid saturation relations the van Genuchten (1980) characteristic curves are used for both the soil and the rock:
 | [2] |
 | [3] |
where Se is the effective saturation, defined as
 | [4] |
While the van Genuchten model is commonly accepted for soils, its applicability for fractured media has been questioned because of the fundamentally different pore structure of fractures compared with the porous material for which the model was originally developed (Liu and Bodvarsson, 2001). However, because of a lack of good alternatives, the model is used for fractured media as well. For example, in all the studies listed in Table 2, the characteristic curves are given in terms of the van Genuchten parameters. Also in our research, we used them for both the soil and the rock.
In the heterogeneous simulations, capillary strength parameter 
is correlated to the heterogeneous saturated hydraulic conductivity according to Leverett's (1941) scaling rule:
 | [5] |
Here ref is the value of at the reference saturated hydraulic conductivity, Kref = 10–8 m s–1. With Eq. [5] we take into account the fact that small fractures generally exhibit stronger capillarity than large fractures with wider apertures. The parameter m in Eq. [2] and [3] can be related to the spatial distribution of fracture apertures and is assumed to be a constant. This is equivalent to assuming that the variability of fracture apertures is statistically homogeneous (Birkholzer et al., 1999).
Simulations
The numerical model TOUGH2 (Transport Of Unsaturated Groundwater and Heat, Pruess et al., 1999b) was used to solve the governing equations for the unsaturated–saturated flow system. The stochastic conductivity fields were generated as unconditioned sequential Gaussian simulations with the GSLIB software package by Deutsch and Journel (1998). Different realizations were generated by changing the seed number, while keeping all geostatistical parameters constant.
In the heterogeneous case, the size of the numerical elements were taken as equal to the support scale of the underlying heterogeneity (i.e., 0.5 m). This is a rather common practice in stochastic hydrology, both saturated and partially saturated, and while a finer (i.e., "subsupport scale") numerical discretization could also be used to study the numerical accuracy in the heterogeneous simulations, it was not considered necessary in this preliminary study. Further studies should consider the effect of finer meshes. In the homogeneous case the model domain is discretized into square elements with a side length of 0.25 m. In our opinion it is not critical to have the same grid spacing in the homogeneous and heterogeneous cases. If grid effects were to be studied in the homogeneous cases, then an even finer mesh might be needed than the 0.25 m used here. However, the grid effects are likely to be of second-order importance and would be better addressed in a separate future study. Harmonic mean weighting was used at element interfaces for both the absolute and relative permeabilities, as recommended by Pruess et al. (1999b) for steady-state simulations.
In the homogeneous simulations, both the hydraulic conductivity of the rock and the specified head at the lower boundary were systematically varied. Three different values were used for hydraulic conductivity of the rock (K = 10–6, 10–7, and 10–8 m s–1), and for each of these the specified head at the lower boundary was varied from 0 to 10 m. Two different net infiltration rates (150 and 200 mm yr–1) were used. The output variables of interest are the head field and the flux over the soil–rock interface.
In the heterogeneous simulations one mean conductivity value was used for the rock (K = 10–8 m s–1), and for this the lower boundary head was varied, as in the homogeneous case. The net infiltration is 200 mm yr–1. For each set of boundary conditions 100 realizations were simulated. The statistical convergence of the output statistics from the Monte Carlo simulations was usually reached after about 50 realizations. The mean and standard deviation for head values usually converged after 30 to 50 and 40 to 50 realizations, respectively. Convergence of the mean and standard deviation for flux over the soil–rock interface took place after 25 to 55 and 20 to 50 realizations, respectively.
|
RESULTS AND DISCUSSION
|
|---|
Homogeneous Simulations
Simulated fluxes from soil into the bedrock for the different homogeneous cases are summarized in Fig. 3
. The figure shows the flux into bedrock as a function of the specified head at the lower boundary (essentially a measure of the hydraulic gradient) as obtained for the different bedrock conductivities and for the two net infiltration rates. In most cases the lowest vertical hydraulic gradient that can be simulated is that corresponding to a specified head of 9 m at the lower boundary. With higher specified head values at the lower boundary, starting from a head equal to 10 m, there is ponding at the upper boundary of the model (i.e., pressures greater than atmospheric pressure are induced in elements at the ground surface). Since the model cannot handle surface runoff, these cases are not included in the results and further analyses. In one of the cases (the higher net infiltration of 200 mm yr–1 in combination with the lowest bedrock conductivity K = 10–8 m s–1), ponding took place already with a specified head of 6 m at the lower boundary.
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 3. Flux into the rock domain as a function of the specified hydraulic head at the lower boundary. Results are for homogeneous simulations with different net infiltration rates and various hydraulic conductivities of the rock.
|
|
Inspection of the results in Fig. 3 indicated the following:- Depending on the specific hydraulic conditions imposed, the flux into the bedrock is 44 to 100% of the net infiltration.
- The lowest percentage (44%) corresponds to the case with lowest vertical hydraulic gradient (highest specified head), lowest conductivity of the rock, and the lower net infiltration. As can be expected, the flux into the bedrock increases with increasing vertical hydraulic gradient and increasing hydraulic conductivity of the rock.
- With the higher net infiltration (200 mm yr–1), larger hydraulic gradients are needed to take all the net infiltration (100%) into the bedrock.
Obviously, in all the cases where 100% of the net infiltration is going into the bedrock, the water table is also situated in the bedrock. Figure 4a
shows the position of the water table for one example conductivity contrast and for different values of lower boundary head. It can be seen that for the case with 100% of the net infiltration going into the rock, the water table is below the soil–rock interface (lower dashed line in Fig. 4a), while for the other cases it is above. Also for these latter cases, the higher the imposed hydraulic gradient in the rock (the lower the hydraulic head at the lower boundary), the lower the elevation of the water table.
View larger version (20K):
[in this window]
[in a new window]
|
Fig. 4. Position of the water table for a net infiltration of 200 mm yr–1 and rock hydraulic conductivity of 10–8 m s–1 based on (a) the homogeneous simulations and (b) the mean of the heterogeneous simulations. Lines 2 through 6 correspond to various specified values for hydraulic head at the lower boundary (units in meters). The upper boundary of the model and the soil–rock interface are indicated with broken lines.
|
|
Heterogeneous Simulations
The heterogeneous simulations were performed for one of the base cases of the homogeneous ones, corresponding to a mean hydraulic conductivity of the rock equal to 10–8 m s–1 and a net infiltration equal to 200 mm yr–1. Hydraulic head at the lower boundary of the model was varied as in the homogeneous cases. Frequency distributions of the flux into the bedrock for the different boundary head conditions are shown in Fig. 5
. The flux is given as normalized with respect to the net infiltration, so a value equal to 1 corresponds to 100% of the net infiltration going into the bedrock. In Fig. 6
the results are shown along with the results from the corresponding homogeneous case (dashed line). The results from the individual realizations are shown as dots, and the geometric mean of all the realizations is shown as a solid line. Inspection of the results indicates the following:- There is a large variation in the results among different realizations. For the same base case the flux into bedrock can be from 50 to 100% of the net infiltration, which is a larger variability than that observed due to different boundary conditions in the corresponding homogeneous case (70–100%).
- The mean flux into rock decreases with increasing boundary head (decreasing vertical hydraulic gradient), but the effect is not as pronounced as in the corresponding homogeneous case. In other words, the dependency on vertical hydraulic gradient is not as strong when heterogeneity is taken into account.
- At high vertical hydraulic gradients, the mean flux from the heterogeneous realizations is smaller than the corresponding homogeneous value.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 5. Frequency distribution of the total flux for the soil–rock interface for (a) 2 m, (b) 3 m, (c) 4 m, (d) 5 m, (e) 6 m, and (f) 7 m specified hydraulic head at the lower boundary of the model.
|
|
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 6. The total flux through the rock domain as a function of specified hydraulic head at the lower boundary. The dots correspond to individual realizations, the solid line to the arithmetic mean of the realizations. The error bars indicate the standard deviation. The dashed line shows results from the homogeneous simulations.
|
|
To analyze the reasons for the difference we can look at the solutions from heterogeneous simulations in more detail. Figures 7 and 8
show characteristic results from two heterogeneous simulations. Figure 7 represents a realization with a large flux into the rock, and Fig. 8 represents one with a small flux. In both cases we see that unsaturated zones occur also below the water table. This happens in cases with high vertical hydraulic gradients and in locations with large hydraulic conductivity contrasts where low-conductive regions overlie high-conductive ones (i.e., tight intact rock overlying large fractures). This phenomenon is not captured in the homogenous simulations where the entire rock domain below the water table remains saturated. The fact that the heterogeneous model can generate these unsaturated zones and the homogeneous model cannot is also a likely explanation of the difference between the mean value from the stochastic realizations and the corresponding homogeneous value, as observed in Fig. 6.
To further evaluate the cause of the differences among the individual realizations, the geometric mean of local unsaturated hydraulic conductivities (KrelKabs) is plotted as a function of total flux into the rock (Fig. 9)
. It is well accepted that for two-dimensional saturated flow the geometric mean of the saturated hydraulic conductivity is a good measure of effective (total) hydraulic conductivity through the domain, and therefore of total flux. For partially saturated flow a similar measure would be the geometric mean of the local unsaturated hydraulic conductivities (note that we use "unsaturated hydraulic conductivity" for KrelKabs instead of the commonly used "effective hydraulic conductivity" to avoid confusion with "effective hydraulic conductivity" in the sense of total hydraulic conductivity as used in stochastic theory). Figure 9 shows that while there is a correlation between the geometric mean of unsaturated hydraulic conductivities and the flux, it is very weak, indicating that some other factors besides simply the mean-value effective hydraulic conductivity also influence the flux.
Comparing Fig. 7 and 8, representing cases with high and low total flux, respectively, we find that in the case with low flux there is a more even distribution of high-conductive flow paths, but in the case with high flux there is one very pronounced high-conductive pathway (corresponding to 18% of the total flux). Thus, in addition to the geometric-mean effective hydraulic conductivity of the rock, occurrence of high-conductive connected pathways may also influence the flux into the bedrock at the scale in question. Using the terminology of stochastic hydrology we explain this by the fact that the stochastic field does not fulfill the criterion of an infinite two-dimensional stochastic field, an underlying requirement for the geometric-mean effective hydraulic conductivity to be valid (Gelhar, 1993). It is likely that the variance of the results in Fig. 6 would decrease with increasing domain size (Gelhar, 1993), and the results would be less sensitive to individual realizations. The mean value (solid line in Fig. 6), however, is likely not to be so sensitive to the domain size. For example, Ewing and Gupta (1993a)(1993b) analyzed the effect of domain size on both the variance and the effective mean value hydraulic conductivity by Monte Carlo analysis and percolation theory. According to their results, the mean value did not change with increasing domain size as long as the correlation length was less than about two-thirds of the network length, a criterion well fulfilled in this study for the variograms of unsaturated hydraulic conductivities.
Finally, it should be noted that effective hydraulic conductivity (in the sense of total hydraulic conductivity) of partially saturated media is also a function of capillary pressure (Gelhar, 1993). Therefore, given the strong heterogeneity in water saturation (see Fig. 7 and 8), and subsequently the capillary pressures in our simulated fields, it is not surprising that the correlation in Fig. 9 is not more pronounced.
|
CONCLUSIONS
|
|---|
We modeled the steady-state distribution of water flux in a system where soil overburden is overlying fractured bedrock. We looked at different net infiltration rates, hydraulic gradients, soil–rock hydraulic conductivity contrasts, and most notably, the effect of rock heterogeneity. The rock was first represented as a homogenous deterministic continuum, then as a stochastic continuum with a spatially correlated hydraulic conductivity field.
The results from the homogeneous simulations showed that depending on the specific hydraulic conditions imposed, the flux going into the bedrock varied from 44 to 100% of the net infiltration. The lowest percentage (44%) corresponded to the case with the lowest imposed vertical hydraulic gradient, lowest rock conductivity, and the lower net infiltration. The percentage increased as the magnitude of these variables increased. The range obtained was in good agreement with results from earlier studies. Tiedeman et al. (1998) found from their modeling of the Mirror Lake groundwater basin that about 60% of the flux through the groundwater system involved flow into the bedrock. In their case the hydraulic conductivity of the soil was one order of magnitude greater than that of the rock, and the study involved the entire basin, whereas we looked at a recharge area only. Their applied recharge rate was somewhat higher (300 mm yr–1) but of the same order of magnitude as the net infiltration we used (150–200 mm yr–1). Harte and Winter (1995) in turn obtained fluxes through the bedrock of between 50 and 90% of a somewhat smaller total recharge rate (30–80 mm yr–1) and similar hydraulic conductivity structure compared with that of Tiedeman et al. (1998).
Next we looked at the effect of rock heterogeneity by modeling the rock as a heterogeneous stochastic continuum. The results showed that there was a larger variation in the flux into the bedrock between different stochastic realizations (from 50 to 100%) than due to different boundary conditions (from 70 to 100%). The flux into the bedrock was correlated, but only weakly, to the geometric mean of unsaturated hydraulic conductivities. Given the strong heterogeneity in simulated capillary pressures, the fact that effective hydraulic conductivity of partially saturated media is also a function of capillary pressure (Gelhar, 1993) partly explains the weak correlation. Another explanation may be the occurrence of well-connected high-conductive pathways. It is likely that the variability between the different realizations would decrease with increasing domain size, but because the scale of the domain and the scale at which heterogeneity was varied are both based on typical settings in Scandinavia, we can conclude that some uncertainty in predicting fluxes into bedrock will remain, even when the statistical characteristics of the rock are known.
The heterogeneous simulations also showed formation of local unsaturated zones below the groundwater table. These occurred in the case of large vertical hydraulic gradients and in locations with large hydraulic conductivity contrasts, with low-conductive regions overlying high-conductive ones. This phenomenon was not captured with the homogenous models, where the entire rock domain below the water table always remained saturated. Subsequently, in cases where hydraulic gradients are large the homogeneous models may overestimate flux into bedrock, even in attempts to obtain estimates of the mean, expected-value behavior. Thus, when modeling flow through a fractured rock surface, it is always especially important to account for the heterogeneity effects when gradients are large.
The studies presented here should be considered preliminary estimates of how rock heterogeneity may influence the lateral diversion of water flow at the interface between a soil layer and underlying fractured bedrock. Future studies should consider the effects of different types and degrees of heterogeneity, as well as the role of varying the boundary conditions, especially the head-dependent boundary that influences the base case of the lateral diversion in our model. Furthermore, the behavior is likely to be affected if the three-dimensionality is taken into account. All transient effects are also likely to be important, including rainfall and snowmelt events, which are typically highly transient and season dependent. Finally, well-instrumented field sites are necessary to further investigate relevant processes, such as the occurrence of local unsaturated zones in systems with strong hydraulic gradients.
|
ACKNOWLEDGMENTS
|
|---|
The authors would like to thank Allan Rodhe for valuable comments. We appreciate the comments and suggestions from Stefan Finsterle and three reviewers that greatly improved the manuscript. The work described in this paper was carried out in a project supported by the Geological Survey of Sweden (SGU) (contract 03-1139/2000).
|
REFERENCES
|
|---|
- Alley, W.M., R.W. Healy, J.W. LaBaugh, and T.E. Reilly. 2002. Flow and storage in groundwater systems. Science 296:1985–1990.[Abstract/Free Full Text]
- Andersson, J.-E., R. Nordqvist, G. Nyberg, J. Smellie, and S. Tirén. 1989. Hydrogeological conditions in the Finnsjön area, SKB AR 89-24. Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm.
- Armitage, P., D. Holton, N.L. Jefferies, B.J. Myatt, and P.M. Wilcock. 1996. Groundwater flow through fractured rock at Sellafield. Nuclear Sci. and Technol. Ser. EUR 16870 EN. E.C. Publ., Eur. Commission, Brussels.
- Birkholzer, J., G. Li, C.F. Tsang, and Y. Tsang. 1999. Modeling studies and analysis of seepage into drifts at Yucca mountain. J. Contam. Hydrol. 38:349–384.
- Deutsch, C.V., and A.G. Journel. 1998. GSLIB: Geostatistical software library and user's guide. 2nd ed. Oxford Univ. Press, New York.
- De Vries, J.J., and I. Simmers. 2002. Groundwater recharge: An overview of processes and challenges. Hydrogeol. J. 10(1):5–17. doi:10.1007/s10040-001-0171-7
- Doughty, C. 2000. Numerical model of water flow in a fractured basalt vadose zone: Box Canyon site, Idaho. Water Resour. Res. 36:3521–3534.
- Ewing, R.P., and S.C. Gupta. 1993a. Modeling percolation properties of random media using a domain network. Water Resour. Res. 29:3169–3178.
- Ewing, R.P., and S.C. Gupta. 1993b. Percolation and permeability in partially structured networks. Water Resour. Res. 29:3179–3188.
- Finsterle, S. 2000. Using the continuum approach to model unsaturated flow in fractured rock. Water Resour. Res. 36:2055–2066.
- Fischer, U., B. Kulli, and H. Flühler. 1998. Constitutive relationships and pore structure of undisturbed fracture zone samples with cohesionless fault gouge layers. Water Resour. Res. 34:1695–1701.
- Gburek, W.J., and G.J. Folmar. 1999. A ground water recharge field study: Site characterization and initial results. Hydrol. Processes 13:2813–2831.
- Gelhar, L.W. 1993. Stochastic subsurface hydrology. Prentice Hall, Englewood Cliffs, NJ.
- Gimmi, T., M. Schneebeli, H. Flühler, H. Wydler, and T. Baer. 1997. Field scale water transport in unsaturated crystalline rock. Water Resour. Res. 33:589–598.
- Harte, P.T., and T.C. Winter. 1995. Simulations of flow in crystalline rock and recharge from overlying glacial deposits in a hypothetical New England setting. Ground Water 33:953–964.[ISI]
- Lee, J.-Y., and K.-K. Lee. 2000. Use of time series data for identification of recharge mechanism in a fractured bedrock aquifer system. J. Hydrol. (Amsterdam) 229:190–201.
- Leverett, M.C. 1941. Capillary behavior in porous solids. Trans. SPE 142:152–169.
- Liedholm, M. 1991. Technical note no. 19. In M. Liedholm (ed.) Conceptual modeling of Äspö Technical Notes 19-32. SKB Progress Rep. 25-90-16b. Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm.
- Liu, H.-H., and G.S. Bodvarsson. 2001. Constitutive relations for unsaturated flow in a fracture network. J. Hydrol. (Amsterdam) 252:116–125.
- Lundin, L.-C., S. Halldin, A. Lindroth, E. Cienciala, A. Grelle, P. Hjelm, E. Kellner, A. Lundberg, M. Mölder, A.-S. Morén, T. Nord, J. Seibert, and M. Stähli. 1999. Continuous long-term measurements of soil–plant–atmosphere variables at a forest site. Agric. For. Meteorol. 98–99:53–73.
- Neuman, S.P. 1987. Stochastic continuum presentation of fractured rock permeability as an alternative to REV and fracture network concepts. p. 533–561 In I.W. Farmer et al. (ed.) Proc. 28th U.S. Symposium on Rock Mechanics, Brookfield, VT. 1987. A.A. Balkema, Lisse, The Netherlands.
- Niemi, A., K. Kontio, and A. Kuusela-Lahtinen. 2000. Hydraulic characterization and upscaling of fracture networks based on multiple-scale well test data. Water Resour. Res. 36:3481–3497.
- Norris, S., L.E.F. Bailey, M.M. Askarieh, and G.E. Hickford. 1997. An assessment of the post-closure performance of a deep waste repository at Sellafield. Vol. 2. Nirex Sci. Rep. S/97/012. Nirex, Harwell, UK.
- Öhman, J., and A. Niemi. 2003. Upscaling of fracture hydraulics by means of an oriented correlated stochastic continuum model. Water Resour. Res. 39(10):1277. doi:0.1029/2002WR001776.
- Pruess, K., B. Faybishenko, and G.S. Bodvarsson. 1999a. Alternative concepts and approaches fore modeling flow and transport in thick unsaturated zones of fractured rock. J. Contam. Hydrol. 38:281–322.
- Pruess, K., C. Oldenburg, and G. Moridis. 1999b. TOUGH2 user's guide. Version 2.0. Rep. LBNL-43134. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Tiedeman, C.R., D.J. Goode, and P.A. Hsieh. 1998. Characterizing a ground water basin in a New England mountain and valley terrain. Ground Water 36:611–620.[ISI]
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
- Wanfang, Z., H.S. Wheater, and P.M. Johnston. 1997. State of the art of modeling two-phase flow in fractured rock. Environ. Geol. 31(3–4):157–166.
TOP
ABSTRACT
INTRODUCTION
