Spatial Variability of Hydraulic Properties and Sediment Characteristics in a Deep Alluvial Unsaturated Zone

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Identifyingandanalyzingvariability in soil properties is the first step in assessing vadose zone flow dynamics and in predicting the fate of solute transport in soils. Variability in soil properties is a critical element across wide areas of research including the improvement of agricultural practices, environmental protection in agricultural areas (Robert et al., 1996), environmental protection at potential waste discharge sites, land–atmosphere interactions, and global climate change (Green et al., 2007). It has been shown, theoretically and in field experiments, that the spatial variability of soil properties can significantly impact the amount of solute leaching in soils and that solute concentrations may vary significantly across short distances as a result of soil heterogeneity (e.g., Lund et al., 1974; El-Kadi, 1987; Harter and Yeh, 1996; Russo et al., 1997; Desbarats, 1998; Minasny et al., 1999; Bagarello et al., 2000; Coutadeur et al., 2002). This may lead to large amounts of solute being leached quickly in some portions of the soil profile, while others retain the solute for very long periods of time.

Studies have recently focused on quantitatively assessing the variability of soil physical properties between, within, and across morphologically defined soil series taxonomic units (Makkawi, 2004; Iqbal et al., 2005; Herbst et al., 2006). Duffera et al. (2007) conducted two mixed-model analyses and principal component analysis to describe the field-scale horizontal and vertical spatial variability of soil physical properties and their relations to soil map units in typical southeastern U.S. Coastal Plain soils. Their results indicated that some of the soil physical properties such as soil texture, soil water content (θ), and plant-available water showed significant horizontal spatial structure and were captured by soil map units. Other variables such as bulk density (ρ_b), total porosity (ϕ), and saturated hydraulic conductivity (K_s) did not show much spatial correlation in the field and were unrelated to soil map units. Iqbal et al. (2005) used geostatistical analysis and constructed semivariogram functions of soil physical properties (e.g., ρ_b, K_s, and θ). They used the structured semivariogram functions in generating fine-scale kriged contour maps and indicated that a sample spacing of 400 m provided an adequate measure to define the spatial structure of soil texture and a 100-m sampling range was adequate for soil hydraulic properties and bulk density.

Most of these studies have focused on the variability in the root zone horizons, which constitutes only the upper 1 to 2 m of the soil. In many agricultural areas, particularly in arid and semiarid regions, groundwater levels may reach up to 30 m deep or more. Few studies have surveyed soil properties to such depths or monitored solute leaching to a deep water table (e.g., Onsoy et al., 2005; Baran et al., 2007). Also, most of the intensively used agricultural areas in semiarid regions around the world are located in large to very large basins with surficial geology dominated by continental, geologically young, unconsolidated deposits of typically very heterogeneous structure. Our current understanding of the variability of hydraulic properties and their impact on solute fate and transport below the root zone is therefore limited and based on greatly simplified models.

In this study, we characterized the variability of the geologic (sediment) and hydraulic properties throughout a 16-m-deep, alluvial vadose zone consisting of unconsolidated, alluvial deposits typical of the alluvial fans of the eastern San Joaquin Valley, California. In a novel approach, we used geologic characterization to replace the soil series description common in other spatial variability studies of soil hydraulic properties. The study was implemented at a research orchard at the University of California Kearney Agricultural Center. The Kearney site provides a unique, extensively sampled and characterized field site with a well-controlled, long-term fertilization research experiment that was completed just before our intensive deep vadose zone sampling campaign (Onsoy et al., 2005). Specifically, this study (i) determined the hydraulic properties of the deep unsaturated zone and their relationship to sedimentary facies and texture, and (ii) statistically and geostatistically analyzed these hydraulic properties. The results provide the basis for an analysis of flow and solute transport in deep alluvial sediments, which will be the main focus of a subsequent study.

Details of the field site characterization efforts have been described in Harter et al. (2005) and Onsoy et al. (2005). Briefly, the Kearney site, a former ‘Fantasia’ nectarine [Prunus persica (L.) Batsch var. nucipersica (Suckow) C.K. Schneid.] orchard, is located on the east side of the San Joaquin Valley (Fig. 1 ), approximately 30 km southeast of Fresno, CA, at the University of California Kearney Research Center. The site is about 0.8 ha (2 acres) and is located on the Kings River alluvial fan, a highly heterogeneous sedimentary system consisting of coarse channel deposits, coarse to fine overbank deposits, fine floodplain deposits, paleosols, and fine eolian deposits. Sedimentary layers exposed to the surface for a sufficient amount of geologic time have developed soil profiles with distinguishable horizons. The type of sedimentary layering, the paleosols encountered, and the range of soil textural classes present at this site are rather typical for many areas in the San Joaquin Valley that have deep vadose zones (Weissmann et al., 2002). Similar alluvial conditions are also found in the Salinas Valley and in the desert basins of southern and southeastern California. As in many surrounding areas, groundwater levels at the Kearney site are significantly deeper than the root zone. Since 1970, water levels have fluctuated between approximately 11 and 20 m below the surface. In 1997 (the time of sampling), the unsaturated zone was approximately 16 m thick.

The orchard was planted in 1975 and had four cultivars of nectarines. A fertilization experiment was conducted at the orchard with five levels of fertilization applied in a randomized complete block design. Details of the orchard geometry and the fertilization experiment can be found in Johnson et al. (1995) and Onsoy et al. (2005).

During 1997, on completion of the fertilizer experiment, three subplots were selected for detailed sampling and intensive data analysis (Fig. 2 ). Approximately 900 m of geologic material were obtained from 62 continuous soil cores drilled to the water table (∼16 m), with 18 to 19 cores collected at each of the three subplots (Fig. 2). An additional north–south transect throughout the entire orchard, consisting of six cores spaced 12 m apart, was sampled to obtain estimates of heterogeneity at the scale of the entire orchard.

While water content and NO₃ distributions were analyzed from samples of all 62 boreholes (Onsoy et al., 2005), only 19 of the 62 boreholes were used for the analysis of soil hydraulic properties and laboratory texture analysis. A complete sedimentologic description by color, texture, grain size and roundness of sands and gravels, and sediment structure was performed on all cores. Texture was identified using field estimation methods of the Soil Conservation Service (1994); a Munsell color chart was used to identify color. Cross-bedding, mottling, clay coatings, aggregate presence and size, and cementation or concretions were identified in the sedimentologic description. Individual sediment beds were identified and logged based on the aggregate of these descriptors.

Hydraulic characterization was performed on 120 undisturbed core samples taken at various depths from the 19 core locations. Hydraulic characterization included determination of saturated hydraulic conductivity, grain size distribution, and measurement of the dependence between unsaturated hydraulic conductivity, moisture content, and soil water pressure. Additional measurements such as bulk density and sand, silt, and clay fractions were also included. Soil moisture was measured in the field on disturbed core samples taken adjacent to the undisturbed core samples.

Saturated hydraulic conductivity was measured using the constant-head method (Klute and Dirksen, 1986). The Division of Agriculture and Natural Resources analytical laboratory determined soil texture based on the percentages by weight of sand, silt, and clay (hydrometer method, ASTM, 1985). Bulk density was obtained gravimetrically from the undisturbed cores. The soil water retention and unsaturated hydraulic conductivity relations are basic elements necessary for the simulation and prediction of flow and transport in the vadose zone. A multistep outflow technique (Eching and Hopmans, 1993) was used to determine these relationships.

The principle of the multistep outflow technique is to observe the water outflow from an initially saturated soil core sample along with soil water suction changes in that sample at increasing steps of dryness. The method has two components: (i) implementation of a laboratory experiment, and (ii) computer analysis of the laboratory experiment to determine the hydraulic parameters of the unsaturated hydraulic conductivity function and of the soil water retention curve via inverse modeling.

For the laboratory experiment, a 10-cm-long, saturated, undisturbed sample was placed into a pressure–suction chamber under atmospheric pressure conditions. During the experiment, the air pressure was increased in several discrete steps during the course of several days (typical for sands) to several weeks (typical for clays). Each stepwise increase in air pressure forces water to flow out of the soil core sample until the soil water suction in the pores matches the applied air pressure. Using high-precision instrumentation, we monitored how quickly the soil pressure inside the core changed in response to each pressure step and we monitored the outflow rate from the core with time. The core was instrumented with a tensiometer at its center measuring the soil water suction. A burette connected to the core captured the outflow. Soil pressure and outflow were recorded automatically using pressure transducers connected to the tensiometers and burettes, and the data were sent to a computer. After completion of each experiment, the measured data were converted into meaningful units using laboratory-derived calibration curves (Tuli et al., 2001).

To streamline the implementation of the multistep laboratory experiments, the 120 samples were arranged into 12 sets (or runs) of 10 samples (or cells) running in parallel. The implementation of a single set (10 parallel laboratory multistep experiments) typically took 3 to 6 wk including setup and take-down, depending on the texture of the samples. Coarse-textured samples are typically faster to run than fine-textured samples due to their faster response to pressure changes.

The multistep outflow experiments were successfully completed for 118 undisturbed cores. Due to a variety of experimental complications and errors, however, the multistep outflow data for 21 soil cores were unusable, resulting in 97 viable samples for the inverse modeling process.

Among the 97 samples, transient data were unavailable for all the samples in Runs 7 and 8 (20 samples). Due to transducer failure, seven samples in Run 4 also had unusable transient data and thus the total number of data sets was reduced to 70 samples. For the 27 samples with missing transient data, there existed handwritten data for the equilibrium conditions between pressure steps during the outflow experiment. Implementation of the inverse modeling for these 27 samples and the remaining 70 samples with transient data was described in detail in Tuli et al. (2001).

After obtaining van Genuchten parameters for all soil samples, statistical and geostatistical analyses were performed on the hydraulic parameters. Distribution goodness-of-fit was tested using a standard Kolmogorov–Smirnov test (Stephens, 1974). Transformations were applied to those parameters found not to follow a normal distribution (see below). Differences in population means of the vector of soil hydraulic properties (hydraulic conductivity, shape parameters, residual and saturated water content) among the various classes of sediments were tested using multivariate analysis of variance (MANOVA; Bray and Maxwell, 1982). Multivariate ANOVA allows testing of group vector means rather than testing group means of a single variable as in traditional ANOVA. The Wilks' lambda multivariate test was used to determine statistical significance. The test generates an F value for which the null hypothesis (that group means are not different) was here rejected at P values <0.05 unless otherwise indicated. Multivariate ANOVA requires that the hydraulic properties are normally distributed and assumes that variances across all sediment classes are similar (homoscedasticity). Homoscedasticity of each hydraulic parameter across sediment classes was tested using Levene's test (Milliken and Johnson, 1984). Where the Wilks' lambda test indicated that significant overall differences existed between sediment classes, two post hoc tests were performed to determine which sediment classes were statistically similar and which were statistically different. The first post hoc test is the Tukey honestly significant difference (HSD) test, which is used for unequal sample sizes (between groups or sediment classes). Alternatively, the Newman–Keuls (another post hoc test) is used to make a pairwise determination of whether parameter means between groups (sediment classes) are significantly different. In determining whether differences are statistically significant, both post hoc test methods implicitly account for the fact that, in the post hoc test, multiple pairwise comparisons are made. The Newman–Keuls is a modification of the HSD test that does not require a homoscedastic data set. All statistical analyses were performed using STATISTICA Version 6.1 (Statsoft, 2004).

For the geostatistical analysis of the spatial parameter distribution, we computed standard semivariograms for the different van Genuchten parameters and the results were fitted to a spherical semivariogram model. The software package GSLIB (Deutsch and Journel, 1992) was used and the results were used to evaluate the sampling spacing in the horizontal and vertical directions.

The material obtained in the borehole cores is exclusively composed of unconsolidated sediments. The sediments range in grain size from clay to pebble and cover a wide spectrum of silty, sandy, and loamy sediments in between. The colors of the sediments range from grayish brown to yellowish brown, more randomly to strong brown (no significant reduction zones). The thickness of individual sediment beds varies from <5 cm for some clayey or silty floodplain and alluvial overbank materials to >2.5 m for large sandy deposits associated with buried former stream channels. Both sharp and gradual vertical transitions are present between texturally different units. Based on field descriptions of texture, color, and degree of cementation, five major sediment units were identified: (i) sand, (ii) sandy loam, (iii) silt loam to loam, (iv) silt to silty clay, and (v) paleosol. The relative occurrence of each category as a percentage of the vertical profile length are 17.2% sand, 47.8% sandy loam, 13.8% silt loam to loam, 8.3% silt to silty clay, and 12.9% paleosol. The following is a brief description of the main features of each sediment unit.

The sand (S) is quartz rich and contains feldspar, muscovite, biotite, hornblende, and lithic fragments. Cross-bedding at the scale of a few centimeters could be observed occasionally within fine-grained sand, showing reddish-brown layers intercalated with gray-brown ones. The dominant color of the sand is a light gray to light brown as the brown hue increases with increasing loam content. The mean thickness of sand layers is 1.7 m; however, it can reach as much as 2.5 m. Very coarse sand and particles up to pebble grain size (up to 1 cm) could be observed occasionally at the bottom of sand units but were not present in all cores.

Sandy loam (sL) sediments have usually light olive to yellowish brown color. Some of these sediments are considered to be weakly developed paleosols because of their stronger red-brownish color, root traces, and the presence of aggregates. Mean bed thickness is 50 cm; however, it can be as much as 2 m thick. The sorting is moderate to good. Thin clay or silt layers (0.5–1 cm) occasionally occur in sandy loam units.

Silty loam and loam (tL-L) is usually slight olive brown to brownish gray in color. The bed thickness is at a scale of a few centimeters to tens of centimeters. Fine-grained sediments often show sharp contacts between the units. Changes from one unit to the next exist at small distances. Cross-bedding can more frequently be observed within silty sediments than in fine sands. Root traces and root-shaped mottles are quite common.

Silt and silty clay (T-C) is the finest textured sediments observed in distinct layers. These are distinguished from the tL-L category by the apparent absence of sand. The main color is brownish gray to olive brown. Fine, <1-mm-thick root traces and mottles are also quite frequent in these fine-textured sediments. Much of these occur between 8- and 13-m depth where we observed vertically and laterally quite heterogeneous and relatively thin bedding of varying, mostly fine textures. In this region, the fine-textured units are often thinly laminated to massive and have a mean thickness of 12.8 cm, while the mode is only about 3 cm due to the thickness distribution being highly positively skewed. Many of the sediments in this category have the appearance of glacial rock flour. A thick, clayey silt bed even extends to 50-cm thickness and was observed across the site in most of the cores. The average spacing between these fine-textured strata is 56 cm, with 20% spaced <10 cm apart, 16% spaced 10 to 20 cm apart, and 10% spaced 20 to 30 cm apart.

Paleosols (P) could be recognized in different stages of maturity and play an important role in the geologic interpretation of the profiles. The paleosols separate distinctly different sedimentary deposition regimes or stratigraphic sequences (Weissmann et al., 2002) and were generally observed across the site in all boreholes. They show a brown to strong brown, slightly reddish color, exhibit aggregates, ferric nodules and concretions, and few calcareous nodules, and are identified in drilling logs as hard, cemented layers. They display a sharp upper and a gradual lower boundary, as is typical for paleosols (Retallack, 1990). Clay content decreases downward in the paleosols. The thickness of the paleosol horizons ranges from 50 cm to about 2 m.

Major stratigraphic units (lithofacies) observed in the core profiles were used to construct a large-scale geologic framework for the research site based primarily on the similarities in the sequence of sediment units between individual profiles (Fig. 2). The framework consists of a vertical sequence of lithofacies, where each lithofacies is characterized by its textural, color, and sediment structure composition. The deepest parts of the cores from 15 to 15.8 m display a strong brownish colored, partly clayey paleosol hardpan (P2). From a depth of 12 to 15 m below the surface, the main textural units are sandy loam to fine sandy loam (SL3). Coarse sand and gravel or fine-grained sediments are occasionally right on top of the paleosol. These sediments show a remarkable wetness due to proximity to the aquifer water table. Sediments between about 8- and 12-m depth are vertically and laterally quite heterogeneous with relatively thin bedding (thickness of centimeters to decimeters), consisting mainly of clayey, silty, and loamy material (C-T-L). Between 6 and 8 m below the surface, a sand layer (S2) is found with laterally varying thickness averaging 1.7 m. A weak, mostly eroded paleosol is developed on top of the sand unit. Sandy loam with intercalated sand and clayey and silty material (SL2) is found at a depth about 4 to 6 m below the surface. A 0.2-m- to >1-m-thick paleosol hardpan (P1) occurs at a depth of about 3 to 4 m. Variable sedimentary structures, predominantly consisting of sand (S1), sit on top of the hardpan and have a mean thickness of 0.75 m. Occasionally, the P1 layer is covered by a very thin clay layer (C) with a thickness of <0.25 m. This clay layer was not found in many of boreholes, however, and was not included in any further analysis. Sandy loam and subordinated loamy sand and loam (SL1) are present from the top of that sandy layer to the land surface with an average thickness of 2.5 m (Fig. 3 ).

Stratigraphically, the sediments of this unsaturated zone profile represent (from top to bottom), the Modesto, Riverbank, and Upper Turlock formations, where the upper hardpan represents the Riverbank Paleosol, separating the younger Modesto formation from the older Riverbank formation. The lower hardpan, just above the water table, represents the Upper Turlock Lake Paleosol (Weissmann et al., 2005). While spatially extensive, the absolute elevations of the upper and lower boundaries of the individual facies appear to vary significantly in space, indicating a small but significant slope in the facies boundaries.

Equivalent to using soil horizons in root-zone studies, the extensive and detailed sedimentologic data set was used to guide the selection of undisturbed sediment core locations for later hydraulic analysis with the multistep outflow experiments. The selected core locations were distributed across the various facies and were classified by their respective facies membership. Clear polyethylene terephthalate glycol acetate liners for the soil core collection allowed an in-field determination of the undisturbed core sample location such that a core sample would be located within a facies; however, a core sample did not necessarily consist of a single massive sedimentary layer. The 10-cm cores, in some cases, included multiple identifiable sediment layers, especially in the finer textured facies. Where possible, undisturbed cores were taken from relatively homogeneous, massive sediment structures. The measured texture classification represents the bulk 10-cm core and was obtained only after the multistep experiment was performed (laboratory classification). In some cases, primarily for samples representing the fine-textured facies, field identification of texture and laboratory determination of texture differed somewhat. Laboratory analyses were generally higher in sand content than field descriptions (Table 1 ). It is thought that this reflects the presence of a significant amount of very fine sand in the fine-textured samples, but also the nonuniformity in the texture of individual multistep outflow sediment cores.

To determine, whether sediment classification accounts for some of the observed spatial variability in hydraulic properties, we statistically tested all three sediment classification schemes developed as part of this study (Table 1):

1. Classification Type I is a field description of the sediments as described above. This classification led to five classes [sediment field (five groups)]: sand (S), loamy sand and sandy loam (sL), silt loam and loam (tL-L), silt and clay (T-C), and paleosols (P).

2. Classification Type II is a laboratory analysis of sediment texture, which led to four classes [sediment lab (four groups)]: sand (S), loamy sand (lS), sandy loam (sL), and silt loam and loam (tL-L). The lab analysis did not identify any samples as only clay or silt.

3. Classification Type III is the eight-lithofacies designation [lithofacies (eight groups)] described above: SL1, S1, P1, SL2, S2, C-T-L, SL3, and P2.

The lithofacies classification of the sediment samples could also be simplified by combining multiple, vertically separated facies with similar properties into a single class (e.g., SL1, SL2, and SL3 all become SL), which yielded a Classification Type IV [sediment (four groups)] with four groups or classes: S, SL, P, and C-T-L.

Inverse modeling to obtain the van Genuchten parameters from the multistep outflow experiments was repeated for three sets of initial parameter values. Repeated optimization provides a measure to evaluate the uniqueness of the optimized parameter set. The three optimized parameter sets, corresponding to the three sets of starting values, provided very similar albeit not identical results (Fig. 4 ). We typically observed some small differences in the simulated output between the three different optimal data sets. The final van Genuchten parameters obtained through optimization of either the transient or steady-state data were selected specifically by comparing the mass balance error of the computed flow simulation and the sum of the squared residual of the measured vs. simulated data. The mass balance error for most of samples was found to be <2%.

For 27 of the 97 samples, no transient data were available from the multistep outflow experiments; instead, hydraulic parameters other than K_s were fitted by considering only the steady-state water content and tension at the end of each step. For these 27 samples only, the K_s values used for the statistical and geostatistical analyses were obtained from direct measurement of saturated flow such that a complete data set was available for the analysis. The correlation between measured log K_s and fitted log K_s for the remaining 70 samples is significant but not strong (Pearson's r = 0.31). More importantly, their means and variances are not statistically different (P < 0.05), thus avoiding statistical bias.

For the 97 soil samples, the average K_s was 5.07 cm h⁻¹ and the variance of log K_s was 0.93, which is remarkably similar to the value reported by Russo and Bouton (1992) at the Bet Dagan trench site, which is 2.5 m deep. Optimized saturated hydraulic conductivity showed characteristics of a lognormal distribution, with values varying across more than four orders of magnitude. Saturated hydraulic conductivity values were found to range between 60.5 cm h⁻¹ (associated with a sand sample) and 0.0077 cm h⁻¹ (associated with a silt loam sample). The van Genuchten shape factors α and n also showed characteristics of a lognormal distribution, with maximum and minimum values of α of 0.078 and 0.0013 cm⁻¹, and were also associated with sand and silt loam samples, respectively. Maximum and minimum values of n were 7.47 and 1.12 and were associated with sand and hard pan samples, respectively.

Despite the lack of living roots or agricultural modification of the measured sediments in the deep vadose zone, similarities are apparent between the range and variability of the optimized sediment hydraulic parameters for our site and those reported previously in several studies investigating root-zone hydraulic properties only. The van Genuchten parameters K_s, α, and n were previously reported to be lognormally distributed by Hopmans et al. (1988) at the Hupsel watershed. The Bet Dagan trench data (Russo and Bouton, 1992) also suggested that α was lognormally distributed. The distribution type for these parameters therefore appears to not be affected by the absence of plant roots or agricultural practices in these unsaturated sediments.

The hydraulic parameters for each sample are not entirely independent of each other: hydraulic conductivity (log K_s) is significantly correlated to log α and log n. The two shape parameters of the van Genuchten model, log α and log n, are also significantly correlated, as well as the two transformed bounds of the water content, θ_r′ and θ_s′. There are no significant correlations between the shape parameters and the water content parameters. While significant, the linear correlation coefficients (Pearson's r) are relatively low (0.2–0.4, Table 2 ), except between log K_s and log α, where r = 0.7. The latter confirms findings at a field site involving root-zone samples only, where significant but generally mild correlations were found between saturated hydraulic conductivity and shape parameters (de Rooij et al., 2004).

We next consider the relationship between hydraulic parameter spatial variability and sediment characterization, which can be obtained from field characterization, possibly enhanced by laboratory analysis of grain size distributions. The lithofacies model constructed from the detailed field geologic descriptions (Classification Type III) provides a framework for the macroscopic distinction of various sedimentary units that are—in the framework of identifying macroscopic variability—equivalent to the common distinction made between different soil units identified in soil surveys. The lithofacies concept, however, operates within a three-dimensional sedimentary context rather than a two-dimensional soil map. Similar to the approach applied by Allen-King et al. (1998), we used the facies-based approach as a tool to separate macroscopic, large-scale variability from microscopic, Darcy-scale variability of sediment hydraulic parameters given by the van Genuchten model. In particular, we were interested in determining the degree to which field-identifiable sedimentary structure constrained laboratory-derived hydraulic properties commonly used for unsaturated flow modeling. While visual inspection of soil properties has been found to be a weak indicator of spatial variability in hydraulic properties within a given soil (e.g., Nielsen et al., 1973), macroscopic lithofacies is commonly used to distinguish hydraulically distinct groundwater flow units (Davis et al., 1993; Weissmann and Fogg, 1999). Does the spatial variability of hydraulic parameters between field-identifiable sediment or lithofacies types explain a significant portion of the overall variability of the hydraulic parameters? To consider the power of various sediment classification schemes, we also considered whether the various sediment classification schemes differed markedly in capturing the between-sediment variability of the hydraulic parameters.

The mean, its standard error, and especially the standard deviation of the hydraulic properties show significant differences between lithofacies, particularly for log K_s, log n, and log α (Table 3 , Fig. 6 ). The largest values of these parameters were found in the S1 and S2 lithofacies. This outcome is consistent with the grain size distributions, which indicates that the S1 and S2 facies contain coarse-textured sand materials and were found to be significantly different from the other facies due to the near-complete absence of any fine-textured material. In contrast, the lowest K_s values were not surprisingly found in the C-T-L facies. Consistent with its fine-textured materials, the C-T-L facies retains more moisture than other facies, thus it has the highest mean and variance of residual water content. For saturated water content, both the S2 and C-T-L facies have a higher mean than the other facies, but the latter also exhibits the largest variation among all facies. The SL3 has the smallest mean for both residual and saturated water content (Table 3).

Considering the sediment classification obtained directly from the field description (five groups, Classification Type I), the results are mostly consistent with those for the lithofacies (Classification Type III) characterization (Fig. 6): the S facies has the highest log K_s, log n, and log α. The two water content parameters are highest in the finest textured T-C facies, but also high in the tL-L facies. Saturated water content is lowest in the paleosol (P) facies. For the sediment classification based solely on a laboratory analysis of the grain size distribution (Classification Type II), the saturated water content was slightly higher in the coarsest (S) and finest textured groups (tL-L). For the other parameters in this classification, relative differences among means are consistent with other classifications (Table 3).

For the lithofacies classification (Classification Type III), within-facies variations of the hydraulic properties in the C-T-L facies were nearly equal to those found within the entire data set. This further confirms that the C-T-L unit, which forms the thickest stratigraphic layer, is also the most variable unit at this site. The C-T-L unit also represents the largest fraction of the data set (31 samples, almost one-third of the entire data set). When combining sedimentologically similar lithofacies (Classification Type IV), however, the SL group is as large as the C-T-L group, while the remaining third of the data is comprised of 16 S samples and 19 P samples.

One-way MANOVA shows that differences in mean hydraulic properties between the various sediment classes (groups) are statistically significant for all four investigated sediment classification schemes (Table 4 ). The significance of the sediment grouping does not rest in a single sediment class alone. Post hoc tests using pairwise comparison of grouped means for each hydraulic parameter generally showed that there are two to three, and—in the case of the lithofacies classification especially for log α—as many as four homogeneous groups across sediment classes, with significant differences in the parameter means between homogeneous groups (Table 5 ). Many of the parameters are associated with multiple homogeneous groups. For log α, the two S facies and three SL facies associate varyingly with each other and with the two P and the C-T-L facies. Importantly, the structure of the homogeneous groups shown in Table 5 is not consistent between parameters. For the lithofacies characterization, each parameter has a distinctly different association of homogeneous groups (Table 5). Considering the vector of all five hydraulic parameters and their varying homogeneous group associations, the post hoc tests support the hypothesis that each of the identified lithofacies classes is characterized by its own unique vector of hydraulic parameter means.

Sediment classifications with fewer classes provide slightly simpler associations of homogeneous groups. The most persistent contrast in hydraulic properties is that between the log K_s, log α, and log n means of the S (or S1 and S2) class(es) and of the finer textured classes. For the field sediment classification, for example, the S facies for all three of these parameters is significantly different from the four remaining otherwise homogeneous, finer textured facies. The collapsed lithofacies classification (four groups of Classification Type IV) also shows this association pattern.

The hydraulic properties were found to be only partially homoscedastic across sediment groups (using Levene's test), with varying combinations of parameters that show inhomogeneous variance between groups, depending on the selected sediment classification schemes. To control the lack of complete homoscedasticity, we also tested with the more robust Newman–Keuls test. The latter is more conservative in highlighting differences between groups. It results in slightly more homogeneous associations than the Tukey HSD (less contrast, Table 5). We also applied a completely nonparameteric test of significance, the Kruskal–Wallis test, which assumes neither a particular underlying distribution nor homoscedasticity. That test provided similar results to the Tukey HSD test, confirming that the HSD test, in this case, was not strongly sensitive to the lack of homoscedasticity. For the lithofacies Classification Type III scheme, all three tests confirmed that each parameter associated differently across sediment groups. This underscores the earlier observation that the individual parameter group means within lithofacies are all significant in explaining some of the larger scale variability. Similar results were obtained when applying the Newman–Keul or Kruskal–Wallis tests to Classification Types I and II. It was therefore preferable to generate statistics of hydraulic properties separately for each of these major lithofacies, even if they fall within similar textural groups.

This site exhibits exemplary properties suitable for conducting hierarchical geostatistical analysis of hydraulic parameters (Barrash and Clemo, 2002; Ye et al., 2005). A complete hierarchical analysis, however, requires an even more extensive data set than the one available here. For example, Barrash and Clemo (2002) conducted hierarchical geostatistical analysis for 4699 porosity data at five fluvial deposit units in the Boise Hydrogeophysical Research Site. Ye et al. (2005) applied the analysis to 1344 moisture content data from five layers at the Hanford Site in the state of Washington. Onsoy et al. (2005) used nearly 1000 core data to determine the hierarchical (geo)statistical distribution of water content and NO₃ at the Kearney site.

Cost and time constraints in performing the multistep outflow experiment did not allow us to have such an extensive set of data as required for hierarchical geostatistical analysis. Our data set consisted of only 97 data points for the eight different facies. Variogram analysis was performed, instead, on the hydraulic parameter data set as a whole since there were not enough data points in each lithofacies to perform hierarchical geostatistical analysis. Lithofacies boundaries were considered to be texture discontinuities across which hydraulic properties were uncorrelated. For the geostatistical analysis, only data pairs with both points belonging to the same lithofacies were considered.

The sill of the horizontal semivariogram was similar or close to that of the vertical direction for all of the van Genuchten parameters, suggesting that complete variability is generally observed across the vertical section of the facies (no zonal anisotropy, Fig. 7 ). Spatial continuity in the horizontal direction, however, is much larger than that in the vertical direction.

Horizontal spatial continuity for all parameters lies in the range between 5 and 8 m. Vertical spatial continuity is in the range of 0.5 to 1.5 m but is also limited by the average facies thickness.

The nugget/sill ratio was calculated and used to determine whether the sampling scheme was fine enough to capture the spatial correlation of the parameters. Iqbal et al. (2005) suggested that a range of c₀/c between 0 and 25% represents a strong spatial dependence, while a range between 25 and 75% represents moderate spatial dependence, and weak spatial dependence is suggested for ranges of c₀/c >75% and indicates that the distance between sampling points is not large enough to capture the spatial continuity of the parameters.

Our analysis shows that the nugget/sill ratio of the horizontal direction for all parameters was <25%, indicating adequate sampling frequency in the horizontal direction for detection of spatial structure. Moderate spatial dependence in the vertical direction is observed, however, suggesting inadequate sampling frequency in that direction. This was attributed to the occasional existence of thin layers of clays within lithofacies units (see above and Fig. 3). Similar conclusions were made by Ward and Gee (2002) at the Hanford Site, where the vertical correlation structure for hydraulic parameters was weak and attributed to the vertical distance between measurement points (typically about 3 m) being larger than the possible correlation length in the vertical direction, estimated at 1 m.

Understanding the heterogeneity of flow processes in the deep vadose zone is critical to assessing the fate and transport of point and nonpoint source pollutants to aquifers in many semiarid and arid regions with a deep water table in unconsolidated sediments. The deep vadose zone, however, lacks the descriptive information of larger scale variability in soil properties that is commonly available for (shallow) soil root-zone studies from soil surveys.

In this study, extensive continuous sediment core information was successfully used to characterize the textural and structural variability throughout the thick unsaturated zone. Eight major lithofacies units were identified at the site. Subfacies structures have also been identified at the millimeter, centimeter, and decimeter scales, particularly in the finer grained sedimentary facies units, which is believed to cause a high variability of hydraulic parameters at these facies.

Our work shows that field characterization of major lithofacies provides a hierarchical framework of soil hydraulic property variability equivalent to, for example, soil mapping units. Two major sources of variability were identified: large-scale variability between the major, field-identifiable sedimentologic facies, and the smaller scale variability within individual sedimentologic facies. The MANOVA results show that the means of the hydraulic parameter vector differ significantly between individual facies, ascertaining the significance of individual field-identifiable facies as a source of hydraulic variability. The standard deviation for log K_s, for example, is 0.93 overall and ranges from 0.6 to 1.0 within individual facies. The spread between the means of log K_s is >1.5 for the lithofacies classification.

Different sediment classification schemes for delineating larger sedimentologic units provide somewhat different results, but we observed qualitative similarities in the structure of hydraulic property means between these schemes. All schemes, including the eight-group lithofacies scheme, provided classes that were hydraulically significantly different from one another. Importantly, sediment texture and structure (e.g., cementation) information was significant in selecting appropriate facies units. Furthermore, significant differences existed in the vector of hydraulic property means even for texturally similar but geologically separated units, therefore justifying a layered lithofacies characterization of the vadose zone as the larger scale hierarchical level.

Within major facies, smaller scale variability was shown to contribute significantly to the overall hydraulic property variability within the unsaturated zone. Saturated hydraulic conductivity and the van Genuchten shape parameters α and n were lognormally distributed, similar to other well-characterized soil and aquifer sites. The data for θ_r and θ_s followed neither normal nor lognormal distributions and were transformed using hyperbolic arcsine and log-ratio transformations, respectively, for the statistical and geostatistical analyses, which are applied only to normally distributed random variables. Vertical small-scale variability was attributed to cross-bedding, thin layering, and intercalated clay and silt beds observed especially in the finer textured facies. The measurement method (multistep outflow experiments) has a support scale of 10 cm (core length), however, which may locally exceed the observed thickness of both the cross-bedding observed in the sand and the intercalated silt and clay layers. Hence, additional variability probably exists at a smaller observable scale than that of the hydraulic characterization implemented here. Horizontal continuity of hydraulic properties (the range of the variogram) was approximately one order of magnitude larger than the vertical range.

Our study provides a novel and extensive quantitative assessment of hydraulic parameter variability in alluvial sediment material and a framework for characterizing hierarchical spatial variability in a deep vadose zone. We expect that similar properties and similar variability exists in facies at other alluvial fan sites, at least within the San Joaquin Valley, even if the facies assemblage or facies sequences are different (Weissmann et al., 1999).

Importantly, the data set provides an opportunity to apply and test stochastic modeling and other upscaling methods (e.g., Harter and Yeh, 1996, 1998; Harter and Zhang, 1999). Highly heterogeneous flow conditions are possibly prevalent at the site, with strong fingering or preferential flow paths channeling much of the water flow and solute transport through a relatively small portion of the unsaturated domain (Onsoy et al., 2005, Harter et al., 2005). Future work using stochastic simulation methods is needed to further elucidate and field validate such findings.

This work was performed with funding from the California Fertilizer Research Program and the University of California Kearney Foundation. We thank in particular Katrin Heeren, Adelphi-Consult GmbH, Berlin, who performed the geological analysis of the cores during drilling, and Gary Weissmann, University of New Mexico, for his invaluable contributions in interpreting the site stratigraphy. We are thankful to Jim MacIntyre and Michelle Denton for implementing and inverse modeling the multistep outflow data. Geoprobe Systems in Salina, KS, generously provided direct push coring equipment for this project.