Measurement and Modeling of Unsaturated Hydraulic Conductivity

The unsaturated zone plays an extremely important hydrologic role that influences water quality and quantity, ecosystem function and health, the connection between atmospheric and terrestrial processes, nutrient cycling, soil development, and natural hazards such as flooding and landslides. Unsaturated hydraulic conductivity is one of the main properties considered to govern flow; however it is very difficult to measure accurately. Knowledge of the highly nonlinear relationship between unsaturated hydraulic conductivity (K) and volumetric water content () is required for widely-used models of water flow and solute transport processes in the unsaturated zone. Measurement of unsaturated hydraulic conductivity of sediments is costly and time consuming, therefore use of models that estimate this property from more easily measured bulk-physical properties is common. In hydrologic studies, calculations based on property-transfer models informed by hydraulic property databases are often used in lieu of measured data from the site of interest. Reliance on database-informed predicted values with the use of neural networks has become increasingly common. Hydraulic properties predicted using databases may be adequate in some applications, but not others. This chapter will discuss, by way of examples, various techniques used to measure and model hydraulic conductivity as a function of water content, K(). The parameters that describe the K() curve obtained by different methods are used directly in Richards’ equation-based numerical models, which have some degree of sensitivity to those parameters. This chapter will explore the complications of using laboratory measured or estimated properties for field scale investigations to shed light on how adequately the processes are represented. Additionally, some more recent concepts for representing unsaturated-zone flow processes will be discussed.


Introduction
The unsaturated zone plays an extremely important hydrologic role that influences water quality and quantity, ecosystem function and health, the connection between atmospheric and terrestrial processes, nutrient cycling, soil development, and natural hazards such as flooding and landslides.Unsaturated hydraulic conductivity is one of the main properties considered to govern flow; however it is very difficult to measure accurately.Knowledge of the highly nonlinear relationship between unsaturated hydraulic conductivity (K) and volumetric water content () is required for widely-used models of water flow and solute transport processes in the unsaturated zone.Measurement of unsaturated hydraulic conductivity of sediments is costly and time consuming, therefore use of models that estimate this property from more easily measured bulk-physical properties is common.In hydrologic studies, calculations based on property-transfer models informed by hydraulic property databases are often used in lieu of measured data from the site of interest.Reliance on database-informed predicted values with the use of neural networks has become increasingly common.Hydraulic properties predicted using databases may be adequate in some applications, but not others.This chapter will discuss, by way of examples, various techniques used to measure and model hydraulic conductivity as a function of water content, K().The parameters that describe the K () curve obtained by different methods are used directly in Richards' equation-based numerical models, which have some degree of sensitivity to those parameters.This chapter will explore the complications of using laboratory measured or estimated properties for field scale investigations to shed light on how adequately the processes are represented.Additionally, some more recent concepts for representing unsaturated-zone flow processes will be discussed.

Hydraulic conductivity measurement
The most direct and most generally reliable measurements of unsaturated hydraulic conductivity are from steady-state flow methods.These methods are seldom applied, however.In the simple gravity-driven implementations they have serious drawbacks: limitation to the wettest soil conditions, and slowness, sometimes requiring months for a single K measurement.The Steady-State Centrifuge (SSC) method, used to measure K() for the samples presented here, extends the range to lower steady-state K measurements by at least three orders of magnitude, and allows at least six points of the unsaturated K relation with  to be characterized for a pair of samples in about five days (Nimmo and others, 2002).
The steady state centrifuge (SSC) method used to measure K() on 40 samples from the Idaho National Laboratory is the Unsaturated Flow Apparatus1 (UFA) version (Conca and Wright, 1998;Nimmo et al., 2002) of the method originally developed by Nimmo et al. (1987).The core samples were sub-cored in the laboratory using a mechanical coring device into a 4.9-cm-long, 3.3-cm-diameter retainer designed specifically to fit into the buckets of the UFA centrifugal rotor.The SSC method requires that steady-state conditions be established within a sample under centrifugal force.Steady-state conditions require application of a constant flow rate and a constant centrifugal force for sufficient time that both the water distribution and the water flux within the sample become constant.When these conditions are satisfied, Darcy's law relates K to  and matric pressure () for the established conditions as follows: where q is the flux density (cm/s), C is a unit conversion factor of 1 cm-water/980.7 dyne/cm 2 (1 cm of water is equal to a pressure of 980.7 dyne/cm 2 and 1 dyne is equal to 1 g cm/s 2 ),  is the density of the applied fluid (g/cm 3 ), is the angular velocity (rad/s), and r is the radius of centrifugal rotation (cm).If the driving force is applied with the centrifuge rotation speed large enough to ensure that d/dr <<  r, i.e., matric-pressure gradients that develop in the sample during centrifugation are insignificant, the flow is essentially driven by centrifugal force alone.The flow equation then simplifies to 2 () The  threshold for which the d/dr gradient can become negligible depends on the hydraulic properties of the medium of interest.Nimmo et al. (1987) presented model calculations showing that the d/dr gradient becomes negligible at relatively low speeds for a sandy medium and at higher speeds for a fine-textured medium.Speeds high enough for this purpose also normally result in fairly uniform water content throughout the sample, permitting the association of the sample-averageand  values with the measured K.
After achieving steady flow at a given q,  is measured by weight and  is measured by nonintrusive tensiometer (Young and Sisson, 2002) or the filter paper method (Fawcett and Collis-George, 1967) in cases where suctions exceed 800 cm.These measurements along with the computation of K, yields a triplet of data (K,) for the average water content within the sample.Repeat measurements with different q and in some cases different rotational speed give additional points needed to define the K(), and  characteristics.There are five samples for which ) was not measured.Three of the samples had gravel at the top surface which prohibited sufficient contact with the tensiometer and the other two simply lack that measurement.K sat was measured using either the standard benchtop falling head method (Reynolds and others, 2002) or the falling head centrifuge method (Nimmo et al., 2002) where a is the cross sectional area of the inflow reservoir (cm2 ), L is the sample length (cm), A is the cross sectional area of the sample (cm 2 ), t is time (initial and final), z is the height of water above the plane in which the sample rotates (initial and final height, cm), and b r is the position of the bottom of the sample (cm).

Hydraulic conductivity estimation
When direct measurements of K() are not obtainable it is possible to estimate using more easily measured properties such as particle size distributions.To estimate K), and K sat parameters (measured and modeled) were combined with Mualem's (1976) capillary-bundle model, one of the most widely used K) models available.Mualem's model infers a poresize distribution for a soil from its  curve based on capillary theory, which assumes that a pore radius is proportional to the  value at which that pore drains.Mualem's model conceptualizes pores as pairs of capillary tubes whose lengths are proportional to their radii; the conductance of each capillary-tube pair is determined according to Poiseuille's law 2 .In this formulation, K() is defined as where K r () is relative hydraulic conductivity.To compute K() for the whole medium, the conductance of all capillary-tube pairs is integrated as Where L is a dimensionless parameter interpreted as representing the tortuosity and connectivity of pores with different sizes, usually given the value 0.5, of saturation),  r is residual water content and sat is saturated water content and is the retention curve with matric pressure expressed as a function of water content.The parameters used in Mualem's model as described above were obtained in several ways for the data presented here (fig 1) measured water retention data were fit with the Rossi-Nimmo (1994) model, 2) measured water retention data were fit with the van Genuchten (1980) model, and 3) water retention data were estimated based on particle size distributions and bulk density using the Rosetta model (Schaap et al., 1998).

Water retention models and hydraulic conductivity estimation
The choice of the water-retention model used to produce the parameters required for hydraulic conductivity estimation also has an effect on the resulting parameters.The Rossi-Nimmo junction (RNJ) model is one that was chosen here to fit the  measurements because this model is more physically realistic over the entire range of  from saturation to oven dryness ( d ) than other parametric models (for example Brooks andCorey, 1964 andvan Genuchten, 1980) that include the empirical, optimized residual water content ( r ) parameter, which is not well defined.According to capillary theory the largest pores are associated with  values near zero and drain first, followed by drainage of successively smaller pores as  approaches  r .With the Brooks and Corey and van Genuchten models there is an asymptotic approach to  r meaning that if it is taken to be >0, the number of small pores approaches infinity at a water content above zero, which is physically unrealistic.The curve represented by the RNJ model does not have a parameter analogous to  r ; the curve goes to zeroat a fixed value of calculated for the conditions of  d .The RNJ model, like many other parametric water retention models, can be analytically combined with the capillary-bundle model of Mualem (1976) to estimate K() (Fayer and others, 1992;Rossi and Nimmo, 1994;Andraski, 1996;Andraski and Jacobson, 2000).
The RNJ model consists of three functions joined at two points (defined as i and  j , fig. 2): a parabolic function for the wet range of  i ≤  ≤ 0), a power law function (Brooks and Corey, 1964) for the middle range of  j ≤  ≤  i ), and a logarithmic function for the dry range of  d ≤  ≤  j ).This model has two independent parameters: (1) the scaling factor for ( o ), and (2) the curve-shape parameter (Sometimes, o is associated with  at which air first enters a porous material during desaturation (referred to as "air-entry pressure"), but, actually, air begins displacing water in the largest pores at a higher (less negative) than  o as evidenced by the departure of  from saturation to the right of  o on the  curve (fig.2).Unlike the model of Brooks and Corey, which holds  fixed between = 0 and the air-entry pressure, the RNJ model produces a smooth curve near saturation, represented by a parabolic function, that allows the pore-size distribution (the first derivative of the  curve) to be represented more realistically.The curve-shape parameter indicates the relative steepness of the middle portion of the  curve, described by the power-law function.Larger  values cause the drainage portion of the curve to appear steeper.
Fig. 2. Example of water-retention  curve showing the components of the curve-fit model developed by Rossi and Nimmo (1994).
The parabolic function applies for  i ≤  ≤ 0, and is represented by: where  sat is saturated water content expressed volumetrically and c is a dimensionless constant calculated from an analytical function involving the parameter (described below)which also is dimensionless.
The power law function applies for  j ≤  ≤  i , and is represented by: The logarithmic function applies for  d ≤  ≤  j , and is represented by: The dependent parameters are calculated as follows: 1 and For convenience, a  d value of -1  10 7 cm-water (the pressure at which the curve goes to zero ) was used in the model fits for all core samples.This is a reasonable value for a soil dried in an oven at 105 o -110 o C under typical laboratory conditions (Ross and others, 1991;Rossi and Nimmo, 1994).
The RNJ model is integrable in closed form for use in the Mualem (1976) hydraulicconductivity model as described below (Rossi and Nimmo, 1994).Relative hydraulic conductivity, the ratio between the unsaturated and saturated conductivity can be expressed as: where and The measured water-retention data also were fit with the empirical formula of van Genuchten (1980) which has the form: where , n, and m are empirical, dimensionless fitting parameters.Using measured and values,  and n parameters are optimized to achieve the best fit to the data.The parameter m is set equal to 1-1/n in order to reduce the number of independent parameters allowing for better model convergence and to permit convenient mathematical combination with Mualem's model (van Genuchten, 1980) as follows Where K sat is saturated hydraulic conductivity and L is a dimensionless curve-fitting parameter.
Most widely-used unsaturated flow and transport models use the van Genuchten model rather than the Rossi-Nimmo junction model to represent .The van Genuchten equation is parameterized by  sat ,  r , , and n, where the scaling parameter for  is  (analogous to  o ) and the curve-shape parameter is n (analogous to ).

Hydraulic property databases
Property transfer models (PTMs) are another way to estimate unsaturated zone hydraulic property data such as  and K().PTMs, which can be based on simple or complex relationships among variables of interest, serve the purpose of estimating hydraulic properties from more easily measured bulk properties such as particle size distribution and bulk density.Published databases of hydraulic properties, such as those of Holtan et al. (1968), Mualem (1976), Nemes et al. (1999), Wösten (1999), are often used in studies when direct measurements are not possible or when data for a large number of samples are required, such as in development and testing of new models and theories or in comparative or regionally extensive analyses.Some PTM development and testing studies use these published databases (Vereecken et al., 1989 and1990;Schaap et al., 1998;Hwang and Powers, 2003), while others use unpublished data or data presented only in parameterized or graphical form (Gupta and Larson, 1979;Arya and Paris, 1981;Wagner et al., 2001).Schaap and Leij (1998) evaluated the effect of data accuracy on the uncertainty of PTMs and concluded that the performance of a PTM depends strongly on the data being used for calibration and testing, however, estimated properties may be sufficient depending on the application for which they are used.Desirable features of a database include (1) high reliability and precision of measurements, (2) high quality, minimally disturbed samples, (3) a large and diverse sample population, (4) consistency in measurement techniques across the data set, (5) a full suite of hydraulic and bulk property data for each sample, and ( 6) ease of use.The database of Perkins and Nimmo (2009) used for the examples in this chapter presents a data set that, though smaller in sample number than many published databases, is ideal in several ways.Sample collection and preparation techniques were selected to ensure minimal sample disturbance, measurements were performed by the same laboratory techniques under highly controlled conditions, and measurements were done by a limited number of researchers.Samples are from diverse geographic, climatic, and geomorphic environments, and the data were originally generated for various research purposes including recharge estimation, simulation of variably saturated flow and solute transport, theoretical studies of porous media, and property transfer model development.Samples were from various soil depths, including many from below the root zone.Other published data sets commonly include samples from shallow depths; about 80% of the samples in the data set of Nemes et al. (1999) come from depths less than 1 m.Additionally, published databases often contain measurements done on disturbed agricultural soils.The data used here for illustrative purposes include bulk density ( b ), particle density ( p ), particle-size distribution, saturated hydraulic conductivity (K sat ), hydraulic conductivity as a function of water content (K()), and water content as a function of matric potential ()).

Data analysis
The data used to illustrate the effect of parameterization on K() estimation and numerical flow simulations is from the database of Perkins and Nimmo (2009) described above.Specifically they are from a core sample from the unsaturated zone at the Idaho National Laboratory (INL) in eastern portion of the Snake River Plain.The medium is sandy in texture with a K sat of 3.90 x 10 -3 cm/s and a porosity of 0.42.Additionally, measurements from 40 INL samples were used to evaluate the error in K() produced by each estimation technique.

Error calculations
The root-mean-square error (RMSE), also referred to as the standard error of the estimate, is used here as a goodness-of-fit indicator between measured and predicted values of K().The parameters for predicting K() were obtained using water retention data fit with the RNJ model, the van Genuchten model, and retention parameters predicted by the Rosetta model.The RMSE is calculated as: where y j is the measured value, ŷ j is the predicted value of the dependent variable, and n is the number of observations.Smaller values of RMSE indicate that the predicted value is closer to the measured value of the variable.K() values span several orders of magnitude which, in effect, unequally weights points in the RMSE calculation, therefore the values were logarithmically transformed prior to calculation.The number of K() points measured for each sample was between three and ten.

Parameter testing with numerical simulation
Parameterized and K() curves, representative of the modeled media, are required input for numerical flow and solute transport simulations, therefore numerical simulations were run in order to assess the influence of the input parameters on modeled results.
Utilizing parameterized unsaturated-hydraulic properties (K() and ) flow was simulated using the U.S. Geological Survey variably-saturated two-dimensional transport model (VS2DT) (Lappala and others, 1983;Healy, 1990;Hsieh and others, 1999) in order to assess the effect of the chosen input parameters.The model was modified to allow for the use of the Rossi-Nimmo water-retention parameters (Healy, personal communication 2006) in additional to the van Genuchten model parameters.VS2DT solves the finite difference approximation to Richards' equation (Richards, 1931) for flow and the advection-dispersion equation for transport.The flow equation is written with total hydraulic potential as the dependent variable to allow straightforward treatment of both saturated and unsaturated conditions.Several boundary conditions specific to unsaturated flow may be utilized including ponded infiltration, specified fluxes in or out, seepage faces, evaporation, and plant transpiration.As input, the model requires saturated hydraulic conductivity, porosity, parameterized unsaturated hydraulic conductivity and water-retention functions, grid delineation, and initial hydraulic conditions.Three simulations based on soil core properties are presented here.Relations between pressure head and water content were represented by functions developed by Rossi and Nimmo and van Genuchten, in all cases using Mualem's model to calculate relative hydraulic conductivity.Parameters were also obtained with the Rosetta model, as described earlier.The 2-by 2-m domain was discretized into 1-by 1-centimeter (cm) grid blocks with a boundary condition chosen to simulate 60 minutes of infiltration at a constant rate of 0.01 cm/s over a 25 cm section at the top left of the domain.Initial hydraulic conditions were specified as uniform water content (10% volumetrically).

Discussion
Errors were calculated for K() for the 40 samples and ranged from 0.21 to 9.45 with the best overall performance achieved using the RNJ model for fitting the measured water retention data where 47.5% of the samples has the lowest RMSE values.The maximum RMSE value for the RNJ model is at least a factor of 4 lower than the other parameterizations.On average, the Rosetta estimations were slightly better than those achieved by fitting the measured water retention data using the van Genuchten model (27.5% vs. 25.0% of the samples having the lowest RMSE values respectively).Table 1 shows the values for all 3 parameterizations for each sample.
For several of the van Genuchten K() estimates, the RMSE values were unusually large (fig. 3).This occurred in cases where few data points were available and the data clustered within a small range in .The van Genuchten model is not physically realistic for the entire range of  from saturation to  d because the model uses  r as an optimized parameter.For the particular cases where nodata are available in the dry range and the measured points slope steeply within a small range in , the curve asymptotically approaches  r starting from near  sat .On the resulting K() curve, K decreases sharply with little change in .The) curve represented by the RNJ model goes to zero at a fixed value of  calculated for the conditions of  d ; therefore, even when few data points are available, the   The model domain, which was 2 m wide and 2 m deep, had uniform hydraulic properties (sandy textured material) and constant infiltration over a 50 cm section at the top left corner for 60 minutes at a rate of 0.01 cm/s.Changes in water content were output at the end of the infiltration period at 11 observation points (fig.5) to assess lateral and vertical water movement.
Fig. 5. Observation points for water content output at the end of the infiltration period.
The Rosetta parameterization, one of the most commonly used methods, had a K() RMSE of 8.06, the highest of the 3 cases.The wetting front was very diffuse compared to the other models and reached the 100 m vertically and laterally to 30 cm beyond the infiltration boundary within the infiltration period (fig.6).With the water retention curve fitting as the basis for obtaining parameters, the maximum water content is known from the curve.The Rosetta model only has textural information therefore the porosity is estimated and the water contents never reach the true known saturation value, which is slightly higher than the model predicts.The van Genuchten parameterization had a K() RMSE of 1.76.The wetting front was very sharp and progressed vertically to 60 cm with some wetting laterally to 10 cm beyond the edge of the infiltration boundary (fig.6).Saturation was reached more quickly than the other 2 cases.The Rossi-Nimmo parameterization had the lowest K() RMSE of the 3 cases at 0.98.The wetting front progressed vertically about 60 cm and laterally to 10 cm beyond the edge of the infiltration boundary (fig.4).The results using known retention curves with the Rossi-Nimmo and van Genuchten models to predict K() are very similar and presumably closest to reality.This analysis shows variable simulation results due to parameterization even though the conceptual model is highly simplified and homogeneous.The effects of soil structure, layering, preferential flow, variable water inputs, nonuniform initial moisture conditions, and other mechanisms that are often found to dominate field situations are not considered here.These effects would further complicate the model uncertainty as influenced by the chosen parameter estimation method in ways that may be difficult to anticipate.The extremely simplified type of analysis presented here assumes that flow through the unsaturated zone proceeds as described by the Richards' equation.Many studies show that there are processes, such as preferential flow, that are not adequately described by Richards' equation (Perkins et al., 2011;Tyner et al., 2006;Köhne and Gerke, 2005;Jaynes et al., 2001;Yoder, 2001;Iragavarapu et al., 1998, Flury et al., 1994).There are recent studies suggesting that hydraulic conductivity and water retention may be important for modeling diffuse matrix flow; however, they may not capture some of the dominant field-scale processes such as preferential flow.These processes may be better represented by the addition of

Fig. 1 .
Fig. 1.Steps used in the estimation of unsaturated hydraulic conductivity.

Fig. 3 .
Fig. 3. Comparison of the RMSE values predicted hydraulic conductivity (log K) based on parameters from curves fits to water retention data with the van Genuchten and Rossi-Nimmo models, and the Rosetta model.The box indicates the interquartile range (25 th to 75 th percentile), the red line is the median value, the whiskers indicate the values that lie within 1.5 of the interquartile range, and the points indicate outliers.Simulation results from the VS2DT model illustrate the effect of parameterization on model results.Figure 4 shows the hydraulic conductivity curves used in the model plotted with measured values.

Fig. 4 .
Fig. 4. Measured and estimated hydraulic conductivity curves used in the VS2DT numerical simulations.

Table 1 .
relation remains somewhat realistic and, in turn, allows for a better estimate of K(.The RNJ model has a much narrower range in values than the other models and also produces no extreme outliers.RMSE vales for 40 samples from the INL with water retention parameters obtained from curve fitting with the van Genuchten and Rossi-Nimmo models, and estimated with the Rosetta model.