HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Vereecken, H. Articles by Harter, T. Search for Related Content PubMed Articles by Vereecken, H. Articles by Harter, T. GeoRef GeoRef Citation Agricola Articles by Vereecken, H. Articles by Harter, T. Related Collections Flow Effective Parameters Geostatistics Upscaling Soil Physics

REVIEWS AND ANALYSES

Upscaling Hydraulic Properties and Soil Water Flow Processes in Heterogeneous Soils

A Review

^a Agrosphere (ICG-IV) Inst. of Chemistry and Dynamics of the Geosphere (ICG), Forschungszentrum Jülich GmbH, D-52425 Jülich, Germany
^b Dep. of Land, Air and Water Resources, 113 Veihmeyer Hall, Univ. of California, Davis, CA 95616-8628
^* Corresponding author (h.vereecken{at}fz-juelich.de)

Received 9 April 2006.

	ABSTRACT

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION CHARACTERIZING SPATIAL... FORWARD UPSCALING APPROACHES THE INVERSE UPSCALING APPROACHES GENERAL CONCLUSIONS AND RESEARCH... REFERENCES

 
This review covers, in a comprehensive manner, the approaches available in the literature to upscale soil water processes and hydraulic parameters in the vadose zone. We distinguish two categories of upscaling methods: forward approaches requiring information about the spatial distribution of hydraulic parameters at a small scale, and inverse modeling approaches requiring information about the spatial and temporal variation of state variables at various scales, including so-called "soft data". Geostatistical and scaling approaches are crucial to upscale soil water processes and to derive large-scale effective fluxes and parameters from small-scale information. Upscaling approaches include stochastic perturbation methods, the scaleway approach, the stream-tube approach, the aggregation concept, inverse modeling approaches, and data fusion. With all upscaling methods, the estimated effective parameters depend not only on the properties of the heterogeneous flow field but also on boundary conditions. The use of the Richards equation at the field and watershed scale is based more on pragmatism than on a sound physical basis. There are practically no data sets presently available that provide sufficient information to extensively validate existing upscaling approaches. Use of numerical case studies has therefore been most common. More recently and still under development, hydrogeophysical methods combined with ground-based remote sensing techniques promise significant contributions toward providing high-quality data sets. Finally, most of the upscaling literature in vadose zone research has dealt with bare soils or deep vadose zones. There is a need to develop upscaling methods for real world soils, considering root water uptake mechanisms and other soil–plant–atmosphere interactions.

Abbreviations: ERT, electrical resistivity tomography • GPR, ground-penetrating radar • PTF, pedotransfer function • TDR, time domain reflectometry

	INTRODUCTION

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION CHARACTERIZING SPATIAL... FORWARD UPSCALING APPROACHES THE INVERSE UPSCALING APPROACHES GENERAL CONCLUSIONS AND RESEARCH... REFERENCES

 
UPSCALING of soil water processes is one of the most important issues in vadose zone research and soil science. Sustainable management of the vadose zone environment as well as the protection of groundwater resources requires an in-depth understanding of the governing soil water flow processes as well as a characterization of its hydraulic properties. There exists, however, a large discrepancy between the relatively small, "local" scale at which soil water fluxes, vadose zone state variables, and vadose zone properties are usually measured, e.g., in soil cores, soil tensiometers, or with soil solution samplers, and the much larger field scale, watershed scale, or regional scale at which hydrologic management decisions are made. For predicting and understanding flow, solute, and energy fluxes at the field, watershed, and global climate modeling scale, vadose zone models must be constructed at scales that are orders of magnitude larger than the scale at which vadose zone properties are determined.

Significant intrinsic (geologic and pedologic) and extrinsic (boundary-driven) spatial variability of vadose zone properties and processes are known to exist at the sub-Darcian pore scale (variations in pore properties across distances of few micrometers), and at the so-called "Darcy" or "local" scale of most laboratory and field instrument methods. Within the vertical soil profile, this scale typically encompasses a few millimeters to centimeters and laterally a few centimeters up to 1 m. Larger scale variations occur due to sedimentary layering vertically (decimeter to meter scale) and between soil mapping units laterally (meter to hectometer scale). Larger scale lateral variability also occurs due to spatial variations in the boundary fluxes across the soil–plant–atmosphere interface at the top of the vadose zone. Such variability is due to variability in precipitation, wind, and moisture flux across a wide range of scales, due to variability in agricultural management practices within and between individual fields, due to land surface geomorphic features, and due to spatial variability in the ecology of the land surface.

In contrast, the grid resolution of watershed, regional, and global climate models is on the order of hectometers to many kilometers laterally. What is the effective mathematical description of vadose zone processes at that scale and how are the parameters of such large-scale models being determined from local-scale information? The pursuit of this question is the focus of so-called upscaling procedures (in the broadest sense). Upscaling techniques provide either effective equations, or effective parameters, or both to describe the vadose zone system's behavior (water, solute, and energy fluxes) at a scale that is relevant for management, assessment, and decision making, while capturing information collected at a much smaller scale and while accounting for the influence of heterogeneities at a scale not resolved by the model (Grayson and Blöschl, 2000; Harter and Hopmans 2004).

In this review, we have defined upscaling as the process that replaces a heterogeneous domain with a homogeneous one in such a manner that both domains produce the same response under some upscaled boundary conditions (Rubin, 2003). The difficulty in upscaling stems from the inherent spatial variability of soil properties and their often nonlinear dependence on state variables. The parameters describing the flow processes in the upscaled homogeneous medium are defined as effective parameters. Often a distinction is made between effective and equivalent parameters. Neuman and Di Frederico (1998) defined effective parameters as those parameters that are used in ensemble-averaged equations (e.g., effective hydraulic conductivity relating the ensemble average flux to the ensemble mean gradient) and equivalent parameters that are present in spatially averaged equations. Practically speaking, effective parameters are an intrinsic property of the homogenized domain (not a function of the particular boundary conditions imposed on the domain), while equivalent parameters are those valid only for a specific set of boundary conditions. In this review we have used the term effective in a more general manner that includes the more specific meaning of equivalent as we refer to the upscaled properties and equations.

In upscaling soil water processes in the vadose zone, four different scales are typically distinguished in the literature (Harter and Hopmans, 2004) (Fig. 1 ). At the pore scale level, water flow is described by either the Stokes or Navier Stokes equations. Description of water flow at this scale requires a complete characterization of the pore geometry of the flow domain, phase boundary conditions, and fluid properties. The next scale is mostly referred to as the local, macroscopic, or Darcy scale. At this scale, Darcy's equation (momentum conservation) and Richards' equation (mass conservation) are typically used to describe soil water processes under variably saturated conditions. The next larger scale is the field scale, which is both a management unit (in agriculture) with its own decision space and also a typical experimental unit to conduct vadose zone research. At that scale, Richards' equation is extensively used to describe and predict variably saturated water flow. Nielsen et al. (1973) were among the first to report on spatial variability of soil hydraulic properties and soil moisture content within an apparently homogeneous field. They found that even seemingly uniform land areas manifest large variations in hydraulic conductivity, while variations in water content, bulk density, and water content are clearly smaller. The study of Nielsen et al. (1973) initiated a large number of field-scale studies exploring the spatial variability of soil water processes and their related properties and state variables (see below). The fourth scale is the regional scale, corresponding to the catchment (watershed) scale in hydrology. At that scale, vadose zone processes are often represented by simple zero-order, black-box, lumped parameter models as part of a complex hydrological system that also accounts for relief, vegetation, weather or climate forces, surface runoff, flood routing, and evapotranspiration, among others.

View larger version (33K): [in this window] [in a new window]   Fig. 1. Forward and inverse upscaling methods to upscale soil water flow processes from the local to the field scale.

We distinguish two major categories of upscaling procedures among those presently available to derive effective hydraulic properties and soil water fluxes. The first group contains methods that derive the parameters of the larger scale process model directly from explicit information about the spatial structure and variability of soil hydraulic properties measured at the small (local) scale. Some of these methods, such as stochastic perturbation and volume averaging, explicitly derive upscaled flow equations with effective hydraulic parameters. Using a forward upscaling model, larger scale properties and fluxes are determined. Subsequently, these approaches will be referred to as forward upscaling (Fig. 1).

The second group, called inverse upscaling, exploits not only information about the controlling system variables (i.e., model parameters, as in the forward approach), but also uses information contained in the measurement of state variables of soil water processes (e.g., soil water tension, water content) and knowledge about boundary fluxes (e.g., totalized recharge, tile drainage). Inverse models estimate model system parameters at the larger scale directly from an assortment of available information or infer information about the small scale (local) spatial structure of the parameter field or both. This involves a model inversion and therefore a priori assumptions about the large-scale process model and about the structure of the small (local) scale parameter distribution.

Some recent approaches expand on the classic inverse approach by further constraining the upscaling problem with secondary, indirect information about primary system and state variables. Such indirect information about, for example, hydraulic properties at the local scale can be obtained from geophysical measurements, remote sensing data, soil maps, and hydrogeological data (so-called "soft data"). This combined inverse modeling approach of using direct and indirect (secondary) measurements of state and system variables will be referred to as data fusion (Yeh and Simunek, 2002).

The review is organized in the following manner: We first give an overview of methods used to characterize spatial variability of local-scale properties. Next, we address the forward upscaling or parameter-driven approach (from the local to the field scale). This includes the three-dimensional stochastic perturbation approach, the concept of the scaleway approach including numerical upscaling, the stream-tube approach, the concept of aggregation, and a very brief overview of other relevant methods such as volume averaging, homogenization, and renormalization. Then we deal with inverse upscaling methods, which make use of large-scale-averaged state values including inverse modeling of large-scale-averaged state values, analysis-of-moment equations that relate spatial moments of parameters and variables, and a geostatistically based estimation of the parameter field from the observed variable fields. This includes a review of data fusion to quantify spatial variability of soil moisture status by using geophysical and remote sensing methods.

	CHARACTERIZING SPATIAL VARIABILITY AT THE LOCAL SCALE

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION CHARACTERIZING SPATIAL... FORWARD UPSCALING APPROACHES THE INVERSE UPSCALING APPROACHES GENERAL CONCLUSIONS AND RESEARCH... REFERENCES

 
The basis for any forward upscaling procedure is a mathematical–conceptual description of the spatial distribution of hydraulic or other vadose zone properties at some small scale at which measurements are taken. Characterization of the spatial variability of properties related to soil water (soil moisture retention characteristic and unsaturated hydraulic conductivity) is typically done at the local scale. For this purpose, small soil cores (centimeter–decimeter scale) are taken in the field and analyzed in the laboratory for their hydraulic properties using standard soil physical techniques. Spatial distributions of hydraulic properties and state variables are determined by taking multiple, closely spaced measurements across a site (Wierenga et al., 1991; Onsoy et al., 2005). A two-step approach is commonly used to describe spatial variability in the vadose zone for flow modeling purposes: Geostatistical methods are used to characterize the structure of spatial variability by use of formal statistical quantities that also define the spatio-temporal relationship between multiple physical variables, e.g., saturated hydraulic conductivity, air-entry pressure, and soil tortuosity factor. Scaling methods are applied before the geostatistical analysis to avoid dealing with multiple, spatially variable but correlated, physical quantities. Scaling methods reduce a multidimensional parameter space into a single-parameter space by taking advantage of the strongly linearizable relationship between soil hydraulic properties and the variable pore-space length scale.

Geostatistical Approaches
The application of geostatistical approaches requires that the hydraulic parameters be described as a function of location. Because it is generally impossible to define or predict parameters at every location in a deterministic sense, their spatial distribution is determined within a geostatistical framework. The geostatistical information is then used as part of an upscaling procedure (e.g., stochastic method) to predict effective hydraulic parameters, upscaled state variables, and fluxes. Geostatistical characterization of a spatially variable quantity is based on the theory of regionalized variables (Matheron, 1971; Journel and Huijbregts, 1978). We will give a brief outline following the description of Russo and Jury (1987a, 1987b). Within a geostatistical framework, a parameter v, e.g., saturated hydraulic conductivity, is considered a random realization of a (three-dimensional) random space function (RSF) V(x), which is defined by a structural or so-called "geostatistical" model (Russo and Jury, 1987a, 1987b) that includes a deterministic and a random, but spatially correlated, component:

[1]

where x is the spatial coordinate vector, m(x) is the deterministic drift function or prior mean, and (x) is the RSF component. A complete characterization of the spatially correlated, random component (x) would involve the definition of an infinite number of probability density functions (pdfs) fⁿ(x₁, x₂, ..., x_n) which represents the n-point pdf of (x) at all locations x₁ to x_n. Since this is not possible, generally two simplifying assumptions are made. The first is the assumption of stationarity, which means that the pdfs are translation invariant: fⁿ(x₁, x₂, ..., x_n) = fⁿ(x₁ + h, x₂ + h, ..., x_n + h). The second is that the first two statistical moments are invariant under spatial translation. The entire RSF is therefore defined by the two-point joint pdf and a complete characterization of the random variable is achieved by describing its first and second spatial moments:

[2]

and

[3]

where C(h) is the spatial covariance of for a lag distance h.

Typically, the spatial structure of the stochastic functions V(x) is defined either by the so-called semivariogram or by the correlogram. The semivariogram is given by

[4]

and the correlation function or correlogram by

[5]

where m is the mean value. The application of stochastic upscaling methods (see below) requires information about the spatial structure of the complete moisture retention characteristic and the unsaturated hydraulic conductivity function. These functions are described by a parameter vector. Its dimension depends on the type of model selected. To apply stochastic perturbation theories, the spatial structure of the individual parameters, but also cross-covariances between the parameters and state variables, needs to be derived. Again, second-order stationarity is typically assumed. Then a cross-covariance function between two stochastic functions V(x) and U(x) is given by

[6]

where the subscripts v(x) and u(x) refer to stochastic functions of two parameters, e.g., air-entry pressure and saturated hydraulic conductivity, and m_i refers to the respective mean values of the stochastic functions.

Geostatistical concepts have been used in many studies to analyze field data of soil hydraulic properties. Most of these studies have focused on deriving the statistical properties of single-parameter distributions (Nielsen et al., 1973; Wierenga et al., 1991; Khaleel and Relyea, 2001), correlations between hydraulic parameters (e.g., Beyers and Stephens, 1983) and semivariograms or covariance functions (e.g., Russo and Bresler, 1981; Jury et al., 1987; Jensen and Refsgaard, 1991; White and Sully, 1992; Russo and Bouton, 1992; Russo et al., 1997; Shouse et al., 1995; Mallants et al., 1996; Vereecken et al., 1997; de Rooij et al., 2004). Some of the studies paid specific attention to the presence of statistical anisotropy in hydraulic parameters (e.g., Russo and Bouton, 1992; Russo et al., 1997), statistical methods to detect and analyze spatial structures (Russo and Jury, 1987a, 1987b; Jury et al., 1987; Ünlü et al., 1990), and the presence of nonstationary behavior of soil hydraulic parameters (Jury et al., 1987; Russo and Bouton, 1992; Russo et al., 1997). It should be noted that information on the spatial correlation of hydraulic properties in the vertical direction is much less common than information on lateral correlation, although the spatial correlation structure of hydraulic parameters in the direction of mean flow has a critical control on flow and transport processes. Moreover, field sites for which spatial correlation lengths of hydraulic properties could be derived are typically located in more arid regions and often in poorly developed soils with a coarser texture (e.g., Entisols) (Ünlü et al., 1990; Russo and Bouton, 1992; White and Sully, 1992; Rockhold et al., 1996). A similar number of studies show variograms of hydraulic parameters with a large nugget, indicating that no spatial structure could be identified (Mohanty et al., 1994; Mallants et al., 1996; Kasteel, 1997; Vereecken et al., 1997; Hammel et al., 1999; Deurer et al., 2000).

Only a few studies analyzed the spatial cross-covariances and semivariances of hydraulic parameters. Ünlü et al. (1990) used the Gardner equation:

[7]

where K(,x) (L T^–1) is the hydraulic conductivity, x the spatial coordinate, t time, (x,t) the pressure head, and (x) a pore size distribution parameter to describe measured moisture retention characteristics and unsaturated hydraulic conductivity obtained from neutron probe measurements from a loam soil. They observed a positive exponential cross-correlation between ln(K_s) and . Vauclin et al. (1994) measured the unsaturated hydraulic conductivities in a loamy soil using the Guelph permeameter and derived the spatial statistics of K_s and in the Gardner equation as well as their cross-covariance. They found a negative spatial correlation between ln(K_s) and ln() at distances up to 24 m.

Mallants et al. (1996) analyzed the spatial variability of the van Genuchten parameters _s, _r, _vg, and n (van Genuchten, 1980) and the saturated hydraulic conductivity using 180 soil cores along a 31-m-long transect in a multilayered soil profile. In addition to the univariate geostatistical analysis, they also derived cross-covariances between the various parameters. They found that the correlation scales of cross-covariances were of a similar magnitude as those pertaining to the direct variograms of cross-correlated parameters, and were typically on the order of 5 to 10 m.

Parameters that describe unsaturated hydraulic properties are not typically measured directly but derived from in situ experimental data using inverse modeling (e.g., data from multistep outflow experiments) or using approximate analytical solutions of the flow equation (e.g., data from infiltration experiments). High sensivity to measurement errors and potential hydraulic model errors lead to uncertainty about the real hydraulic parameter values and to artificial spatial correlation among estimated hydraulic parameters. To what extent the observed variability and spatial correlation structure of estimated hydraulic parameters represents true structural variability versus measurement uncertainty has received little attention. A study by Mertens et al. (2002) demonstrated that parameter uncertainty may be responsible for a significant part of the observed variability of the hydraulic parameters.

In summary, few field data sets presently exist that provide information on the full spatial structure of the moisture retention characteristic and the unsaturated hydraulic conductivity function including the cross-covariances between hydraulic parameters. Data sets that additionally provide the spatial and temporal statistical structure of state variables to validate upscaling theories of soil water processes are even rarer. This is due to the high number of samples needed to characterize the spatial structure, the investment in time and money needed to determine the hydraulic properties despite novel methods like the multistep outflow method, and the costs involved in setting up this type of field experiments.

The Scaling Approach
The scaling approach was introduced in soil science by the pioneering work of Miller and Miller (1956). They derived scaling relationships for matrix potential or pressure head, , and hydraulic conductivity, K, among others, using the microscopic laws for capillary pressure forces and viscous flow in porous media based on a similarity assumption of the pore space structure. Scale factors in general are defined as conversion factors that relate the characteristics of one system to the corresponding characteristics of another system (Tillotson and Nielsen, 1984). These researchers distinguished two approaches in deriving scaling factors of soil water processes: (i) the dimensional analysis, called similitude analysis, which is based on the physical characteristics of a soil, and (ii) the functional normalization approach, which is based on regression analysis. The dimensional analysis approach directly refers to the theory developed by Miller and Miller (1956) and includes fractal-based approaches (Tyler and Wheatcraft, 1990, 1992; Rieu and Sposito, 1991a, 1991b, 1991c; Perrier et al., 1996; Pachepsky et al., 2000). Several examples of applying dimensional techniques to soil properties are given in Tillotson and Nielsen (1984). Sposito (1990, 1998), following up on the work of Miller and Miller, studied the invariance of the Richards equation under scaling transformations. Recently, Williams and Ahuja (2003) and Kozak and Ahuja (2005a) proposed a scaling approach to determine the moisture retention characteristics and hydraulic conductivity curve for a broad range of textural classes based on the work of Gregson et al. (1987) and Ahuja and Williams (1991). Using this approach, Kozak et al. (2005) attempted to scale and estimate evaporation and transpiration of water across soil textures. We briefly review dimensional analysis and the functional renormalization approach. Both approaches have been commonly used to reduce the number of spatially dependent parameters, thereby facilitating geostatistical characterization of heterogeneity in heterogeneous vadose zones.

Scaling of Hydraulic Properties
In the scaling theory developed by Miller and Miller (1956), two soils or porous media are similar when scale factors can be found that transform the behavior of one of the porous media to that of the other (Nielsen et al., 1998). The theory is founded on the concept that the geometry of pore spaces at any location, x, are identical to one another after transforming it to a reference state with a microscopic characteristic length, *. The variability of the pore space can be represented by characterizing the variability in microscopic characteristic length, , and relating it to the hydraulic properties:

[8]

[9]

where _{I =} (x)/* are the scale factors at location x having a mean value of unity and K*() and *() are the reference hydraulic functions describing the reference state of the soil at volumetric water content . Geometric similarity also implies that the scaling factors are identical for both the pressure head and the hydraulic conductivity. Miller similarity is defined at the microscopic scale of the pore space and requires that porous media at different locations differ only in the scale of internal microscopic geometry and therefore have equal porosity. Laboratory studies have shown that the stringent criteria for Miller similarity are met for repacked quartz sand columns (Klute and Wilkinson, 1958; Schroth et al., 1996).

Another approach to scale hydraulic properties was developed by Leverett (1941). He found that experimental soil water retention characteristics from unconsolidated sands could be scaled to one single curve, called the J function, by using the following normalization:

[10]

where S_e is the effective saturation defined as ( – _r)/(_s – _r), with _r and _s the residual and saturated water content, respectively, k the intrinsic permeability, and the surface tension. Through this type of scaling, the moisture retention characteristic of one medium can be related to the moisture retention characteristic of the reference medium containing the same liquid by using

[11]

where the subscript and superscript ref refer to the reference medium. In the case of uniform porosity and using a Leverett scaling relation equivalent to Eq. [9], Eq. [11] yields

[12]

with _L,ref the Leverett scaling factor equal to . Notice that Eq. [12] differs from Eq. [8] only in the definition of the scaling factors. In practice, scaling of the soil water retention curve and application of those scaling factors to the unsaturated hydraulic conductivity is referred to as Miller–Miller type scaling. Scaling of the unsaturated hydraulic conductivity curve (Eq. [12]) and application of those scaling factors to the soil water retention curve is referred to as Leverett type scaling (Oliveira et al., 2006).

Warrick et al. (1977) analyzed the experimental data of Nielsen et al. (1973) using scaling theory and found that the restriction of constant porosity, implicitly contained in the concept of Miller and Miller (1956), could not be maintained due to the large variability in total porosity. This restriction was relaxed by introducing the degree of effective water saturation, S_e. This implies that the water content can be scaled relative to its value at saturation. Figure 2 shows an example of using saturation to scale moisture retention characteristics observed in the field. A total of 154 soil cores were taken to a depth of 1.2 m from four adjacent soil profiles in a Typic Eutrochrept, which has three distinct layers (Kasteel, 1997). Scaling reduced the variation by 81%. The scaling relationships for the hydraulic properties of these so-called Warrick-similar soils (see e.g., Warrick and Nielsen, 1980) are formally similar to the corresponding expressions for Miller similitude after replacing with S_e.

View larger version (21K): [in this window] [in a new window]   Fig. 2. Scaling of measured soil moisture characteristic curves at the field site using effective saturation (data from Kasteel, 1997).

The fractal concept of Mandelbrot (1983) also belongs to the class of scaling methods and has been used to study vadose zone processes. The three characteristics of fractal patterns in the vadose zone are (Perrier et al., 1996): (i) they possess a similar structure across a range of length scales; (ii) their structure is scale independent; and (iii) the structure or pattern cannot be captured entirely by classical geometrical concepts. Tyler and Wheatcraft (1990) showed that the water retention characteristic of a fractal porous medium can be characterized as a so-called Sierpinski carpet. They also demonstrated that the exponent in the Brooks–Corey soil water retention equation (Brooks and Corey, 1964) effectively represents a fractal number that defines the (fractal) distribution of soil pore sizes. Rieu and Sposito (1991a, 1991b, 1991c) further developed a fully self-consistent fractal model of aggregate and pore-space properties for structured soils that includes both a fractal and a nonfractal structure component. A generalized model of both these earlier works was presented by Perrier et al. (1996). Fractal vadose zone structure is consistent with field evidence of power-law variograms of, for example, depth-averaged spatiotemporally varying soil moisture across a 65-ha field site (Green and Erskine, 2004).

The second main group of methods to determine scaling factors is the functional normalization approach, which uses least square regression analysis to derive scale factors from a set of experimental observations. This method is referred to as functional normalization. The objective of functional normalization is to coalesce all measured values of a relationship, e.g., between and or K and , into one reference curve using scaling relations as defined in Eq. [8] and [9]. One single scaling factor is used for each soil sample and for each hydraulic relationship. This yields a frequency distribution of scaling factors. The scaling factors are interpreted in the same way as those derived by dimensional analysis, although they have been obtained without any physical justification. Various methods have been proposed in the literature to determine the scaling factors for matric potential and hydraulic conductivity (Peck et al., 1977; Warrick et al., 1977; Simmons et al., 1979; Russo and Bresler, 1980; Vogel et al., 1991). Hopmans (1987) compared the scaling factors obtained with various methods. These studies showed that the scaling factors for pressure head and hydraulic conductivity have a log-normal rather than a normal distribution. Since the scaling factors obtained using the functional normalization approach are not based on a physical model or certain assumptions about the self-similar structure of the porous medium, there is no constraint on the relation between the scaling factors for the K() and () functions. Warrick et al. (1977) found that the values of needed for scaling differed from those for scaling K. These findings were corroborated in many other field studies (e.g., Ahuja et al., 1984, 1989; Hills et al., 1991).

The previously mentioned methods scale the hydraulic functions independently, with separate scaling factors for the water retention curve and hydraulic conductivity function. Clausnitzer et al. (1992) introduced a method to scale both hydraulic functions simultaneously, yielding a single set of scaling factors that then is consistent with the Miller similitude scaling relationships. They have shown that simultaneous scaling was not always as successful as independent scaling, whereas Hendrayanto et al. (2000) found that simultaneous scaling was the best among five methods tested to scale the hydraulic properties of forest soils.

Vogel et al. (1991) reformulated and extended the concept of similar media to analyze unsaturated or saturated flow in a system of parallel, nonhomogeneous, one-, two-, and three-dimensional soil profiles. Their linear variability concept is based on the functional similarity, which describes the similarity between soil hydraulic properties instead of similarity between internal geometries (Simmons et al., 1979). The total variability in soil hydraulic properties is thought to be composed of a linear and a nonlinear component. The nonlinear component depends on texture and is characterized by the shape of the hydraulic functions. Their analysis also accounts for the time variability of hydraulic properties. Deurer et al. (2000) used soil horizons as "functional similar units". They found that linear and nonlinear scaling factors showed similar spatial structure in the horizontal but not in the vertical direction, with no cross-correlation between the scaling factors.

Oliveira et al. (2006) used data from the Las Cruces Trench site to explore the impact of a single-scaling-factor approach, i.e., Miller–Miller and Leverett scaling, and a multistep stochastic approach on solute transport. In contrast to the multistep stochastic approach, simulation results showed that both single factor scaling techniques, which explicitly assume that the shape factor of the soil moisture characteristic remains constant, seems to underestimate the heterogeneity at the experimental site. Oliveira et al. (2006) further concluded, in accordance with Jury et al. (1987), that both the soil moisture characteristics and the hydraulic conductivity need to be measured to infer scaling factors for both functions separately.

Strengths and Limitations of Scaling
The attractiveness of the scaling approach resides in the fact that the observed variability of scaling factors may be captured in probability density distributions and in spatial-correlation structures. The spatial variability of all of the soil water processes is captured with one, two, or even three scaling variables (see above). Scaling in its various forms, particularly single-variable scaling, has therefore been a convenient approach to analyze the effect of spatially variable hydraulic conductivities on water flow (e.g., Peck et al., 1977; Sharma and Luxmoore, 1979; Hopmans and Stricker,1989; Hopmans et al., 1988; Tseng and Jury, 1993; Roth, 1995; Rockhold et al., 1996; Birkholzer and Tsang, 1997; Kabat et al., 1997; Kim et al., 1997; Shouse and Mohanty, 1998; Wendroth et al., 1999; Deurer et al., 2001). The scaling method has proven its success not only in describing spatial variability of hydraulic properties in many field studies but also as a tool for numerical analyses.

One of the major limitations is the need for detailed information on the variability of soil hydraulic properties requiring extensive measurement campaigns using classical sampling techniques. Despite the substantial progress in soil physical measurement techniques to determine hydraulic functions, determination of all of the hydraulic functions for spatial variability analysis remains time consuming. Also the scaling approach typically requires more than one scaling factor (e.g., two or three—see, e.g., Jury et al., 1987) to describe the observed spatial variation in hydraulic properties. Although convenient, single-parameter scaling leads to significant underestimation of the true variability in the heterogeneous unsaturated flow field and it appears to be unsuitable for upscaling of unsaturated flow or transport processes (Oliveira et al., 2006). Fractal approaches are usually applied to characterize pore- to core-scale hydraulic properties or to quantify or characterize local-scale soil properties (Perrier et al., 2000; Perfect and Kay, 1991; Tyler and Wheatcraft, 1992; Pfeifer and Avnir, 1983; Van Damme and Ben Ohoud, 1990; Friesen and Mikula, 1987; Bartoli et al., 1999; Pachepsky et al., 2000), whereas the application of the fractal approach to analyze soil water processes at the field scale is still limited.

	FORWARD UPSCALING APPROACHES

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION CHARACTERIZING SPATIAL... FORWARD UPSCALING APPROACHES THE INVERSE UPSCALING APPROACHES GENERAL CONCLUSIONS AND RESEARCH... REFERENCES

 
The Scaleway Approach
The scaleway approach is a general conceptual framework rather than a rigorous mathematical framework. In the scaleway approach, upscaling is implemented explicitly and in a stepwise approach across multiple scales (hence, the word scaleway). Spatial variability is considered to exist at multiple scales either continuously (as in fractal systems) or in discrete steps (as in hierarchical systems). The scaleway approach breaks the system into multiple, discrete upscaling steps. At some arbitrary small step related to measurement scale, physicomathematical processes are defined explicitly and the process-relevant structure and topology of the vadose zone is defined exhaustively across the entire domain of interest. The effective properties and processes across that domain become the pixelized "small-scale" processes at the next higher level of upscaling. At the next higher level, a new set of measurements would be needed to define structure and topology at the level across a yet larger domain. The approach was recently presented by Vogel and Roth (2003) as a concept to handle multiscale heterogeneities and to predict flow and transport processes at a specific scale without making assumptions about the nature of the underlying structure because that structure is explicitly considered in this approach. This makes this approach rather flexible when compared with other upscaling approaches (e.g., stochastic perturbation analysis), which are typically limited to analyzing flow processes within specific spatial structures of the heterogeneous parameter field (e.g., joint-Gaussian) and only allow for upscaling across a finite change of scale.

The scaleway approach requires an explicit quantification of the structural properties of the soil at the smaller scale and avoids the need for effective process models at the scale of observation. The scaleway approach aims to bridge the gap between the structure of soils and their function (e.g., transport of solutes) by explicitly characterizing the relevant structural properties. The explicit consideration of the structural properties combined with the forward calculation of the flow process may, however, become numerically very demanding, especially when treating transient flow problems in three-dimensional flow domains. A few studies used high-resolution modeling of flow and transport processes in multi-Gaussian fields (e.g., Ababou and Gelhar, 1988; Desbarats, 1998) and have bridged the gap to stochastic approaches (see below). It is expected that these approaches may become more prevalent with increasing computational capabilities and with the improvement of numerical methods. The use of the scaleway approach is intimately linked to the use of powerful and novel numerical methods and models.

Up to now, this approach has been mainly applied to column-scale processes (Vogel, 1997, 2000; Vogel and Roth, 1998; Kasteel et al., 2000). At this scale, the flow process is typically described using physically based models such as Stokes equation, which is valid at the pore-scale level, or by the three-dimensional Richards equation (Eq. [14]) for variably saturated flow in continuum fields of hydraulic properties that are defined at scales larger than a pore but still smaller than a representative elementary volume. The flow process depends strongly on soil structure, which is determined by the three-dimensional pore geometry. At the pore-scale level, the geometry may be quantified in terms of pore size distributions described by various metrics and the topology of the structure referring to the way the structural units are connected in space. This connectivity can be described by, for example, the Euler characteristic, while the shape of the structural units may be quantified using Minkowski densities and functionals (Vogel, 2002). These Minkowski functionals are measures for quantifying arbitrary binary patterns such as those typically observed in porous media (Vogel et al., 2005; Roth et al., 2005). Similar approaches can be used to characterize the connectivity in continuous fields of hydraulic properties (e.g., permeability) at scales larger than a pore but still smaller than a representative elementary volume.

The need to specifically quantify the structure of the medium and the associated hydraulic parameters requires the use of tomographic methods such as x-ray tomography and image analysis techniques (Tarquis et al., 2003). As these methods do not determine hydraulic parameters directly, proxy relationships between the parameters of interest and the measured variables such as bulk density are needed. These proxy variables can then be used to distinguish structural units. One of the important questions in the scaleway approach is how to derive the () and K() functions. At the core scale, this is solved by an explicit characterization of the pore geometry. Vogel (1997) determined the three-dimensional pore geometry of a soil core using serial sectioning of subsamples obtained from two different structural units. The hydraulic properties can then be estimated from two approaches. The first approach, which is numerically demanding, is based on the numerical solution of the Stokes equation in the full three-dimensional pore geometry. Computationally efficient solutions may be achieved using lattice Boltzmann simulation techniques (e.g., Krafczyk et al., 2000; Tölke et al., 2002). Vogel (1997) followed a second approach that reduces the original pore geometry by defining an equivalent network model consisting of idealized pores described by cylindrical tubes. The network model is generated using the pore size distribution and connectivity derived from the three-dimensional pore geometry of a soil sample. The () and K() are then calculated for this equivalent network using Young–Laplace and Hagen–Poiseuille relationships, respectively. Vogel and Roth (1998) found good agreement between hydraulic properties obtained from this approach and those measured with multistep outflow. An interesting aspect is that in the scaleway approach, connectivity can be integrated in a direct way. Recent studies by Vogel (2000), Zinn and Harvey (2003), and Neuweiler and Cirpka (2005) have demonstrated the importance of accounting for connectivities of structural elements when upscaling water and solute processes, of upscaling from the pore scale to the continuum scale (e.g., Vogel, 2000), and of upscaling heterogeneities at the continuum scale (Zinn and Harvey, 2003; Neuweiler and Cirpka, 2005).

The strengths and limitations of the scaleway concept can only be evaluated in a preliminary fashion, as it has only recently been introduced and only few applications are presently available. It is conceptually very strong but the approach may be limited to small-scale systems. Especially the requirement to explicitly define the structure and the related hydraulic parameters at every point in space may hamper its widespread use; however, the use of proxy variables seems to be an efficient way of circumventing this problem. Application at the field scale would require the use of "tomographic" methods to elucidate the underlying structure (see data fusion approaches, below). Especially the detection of highly connected structures is nearly impossible on the basis of grid observation points but would require a continuous mapping. The potential of hydrogeophysical methods (Vereecken et al., 2003; Rubin and Hubbard, 2005) such as spectral induced polarization (e.g., Titov et al., 2004), seismic methods (e.g., Rubin et al., 1992), ground penetrating radar, GPR (Lambot et al., 2004; Huisman et al., 2003) and electrical resistivity tomography, ERT (Kemna et al., 2003; Vanderborght et al., 2005) may offer the potential to characterize the spatial structure of soil hydraulic properties and flow and transport processes and should be further explored. Insofar as the scaleway concept heavily relies on indirect information about soil structure and topology, it cannot be strictly classified as a "forward upscaling" method, but is an important concept in the application of inverse modeling and especially the data fusion approach.

Stochastic Perturbation Approaches
Theoretical–Conceptual Overview
The scaleway approach requires a complete, deterministic knowledge of the small-scale structure to predict the large-scale behavior. Significant uncertainty typically exists, however, about the exact distribution of hydraulic properties in the unsaturated zone. We introduced geostatistical methods that can be used to statistically describe heterogeneous parameter fields from sparse data measurements. Using such (geo)statistical descriptors of a space variable as input, the so-called stochastic perturbation approaches allow us to statistically quantify upscaled properties, that is, quantify the mean state of the system at the larger scale and the local variance around the mean state at the larger scale, where the variance arises from the incomplete knowledge about the exact structure of the local-scale patterns of the vadose zone.

In the perturbation approach, it is assumed that the local flow q (L T^–1) is described by the Buckingham–Darcy equation:

[13]

where z is the vertical coordinate, which is defined here to be negative in the downward direction, and K(,x) (L T^–1) is the hydraulic conductivity, which we assume here to be an isotropic soil property. In contrast to saturated flow, the hydraulic conductivity in an unsaturated medium is a function of the pressure head. The pressure head field is obtained by solving the Richards equation, which is a mass balance equation:

[14]

where g(x,t) is a sink–source term and C(,x) is the soil water capacity. When the parameters of the local-scale hydraulic properties are treated as random space functions (see above), the flow (Eq. [13]) and mass balance (Eq. [14]) equations are stochastic continuum equations. The solution of the stochastic continuum equations is a stochastic quantity describing a probability distribution for each state variable rather than a deterministic quantity. Using the stochastic equation, moments of the probability distributions of the (dependent) state variables, for example, the water content, pressure potential, and water fluxes, can be computed. Solutions of the stochastic continuum equation often assume a Gardner–Russo parameterization of the hydraulic functions due to its linearity (e.g., Yeh et al., 1985a, 1985b, 1985c; Harter and Yeh, 1998; Russo et al., 1997) described by Eq. [7] and

[15]

where _e is the effective moisture content and m_g is a parameter related to tortuosity. Solutions for other parameterizations have been presented as well. Mantoglou (1992) presented solutions for generalized hydraulic functions. Zhang and Winter (1998) used the Brooks–Corey model (1964) described as

[16]

[17]

[18]

where b is a parameter related to pore size distribution. Tartakovsky et al. (2003) illustrated that the predictability increases when Brooks–Corey or Mualem–van Genuchten parameterizations are used. The parameters [or ln()] and ln(K_s) are assumed to be Gaussian random space functions, whereas _s, _r, and m_g are mostly assumed to be deterministic. A constant m_g implies that the hydraulic functions can be linearly scaled to reference functions (see scaling of hydraulic properties, above). The parameters and variables are expressed in terms of their mean or expected value and a perturbation:

[19]

[20]

[21]

[22]

[23]

where the uppercase characters represent the expected values and the lowercase or primed characters are the perturbations. The perturbations f and a are related to the scaling factors _K and as

[24]

[25]

The hydraulic functions K() and () are written in terms of the expected values and perturbations as:

[26]

[27]

A linearized perturbation equation (i.e., an equation that is linear in the perturbation terms) is obtained by: (i) substituting Eq. [26] and [27] with Eq. [21–23] into Eq. [13] and [14, (ii) neglecting all higher order perturbation terms, and (iii) subtracting the expected value of the equation. For instance, for steady-state flow and constant parameters A and F, the perturbation equation is (Zhang, 2002, Eq. [5.20], p. 226)

[28]

where K_m(H) = exp[F + AH] and J(x) = [H(x) + z(x)] is the mean hydraulic gradient and T refers to transposition. Moment equations are obtained by multiplying the linearized perturbation equation with a perturbation at a location and taking the ensemble mean. Multiplying Eq. [28] by h() and taking the ensemble mean yields

[29]

where C_uv(x,) is the covariance of u at location x and v at location and _x is the gradient at location x. To solve Eq. [29] for the pressure head covariances, the cross-covariances of the pressure heads and the hydraulic parameters, C_fh and C_ah must be derived by solving the following moment equations that are obtained by multiplying Eq. [28] (which is now evaluated at location ) with f(x) or a(x) and taking the ensemble mean:

[30]

[31]

Using Eq. [29] to [31], the pressure head covariance functions can be obtained if the spatial (cross-)covariances of the hydraulic parameters are known and the boundary conditions are defined. Also, the expected pressure head H must be known to solve Eq. [29] to [31]. The value of H can be determined by solving the Richards equation (Eq. [14]) using the expected values of the hydraulic parameters, which is the so-called zeroth order approximation of H: H = ⁽⁰⁾. The advantage of this approach is that H can be calculated directly and that the moment equations become a set of coupled linear equations. The disadvantage is that the effect of spatial variability of the hydraulic parameters on the mean or expected flow variables is not considered.

To describe the effect of the spatial variability of hydraulic properties on the mean flow process, higher than first-order perturbation terms are considered in the expansion of the flow equation (Eq. [13]). Neglecting higher than second-order perturbation terms, the expected value or the ensemble average of the water flux Q_i is

[32]

where E is the expected value. Using Eq. [32], an effective hydraulic conductivity K_eff,ij can be defined as

[33]

It can be shown that when the coordinate axes align with the principal axes of the structural heterogeneity of the hydraulic parameters, the effective conductivity tensor becomes a diagonal matrix (i.e., K_eff,ij = 0 for i j). Using the expected or ensemble averaged fluxes in a mass balance equation, an upscaled flow equation, which describes the mean or expected water contents, fluxes, and pressure heads, is obtained:

[34]

where is the expected water content. In Eq. [34], local sinks or sources are not included but these could be implemented in a straightforward manner. The effective soil water capacity C_eff can be defined and approximated as:

[35]

where E(h) is the covariance between the soil water capacity and matric potential head perturbations. This term has been neglected in most cases since its effect on the large-scale flow is not as important as the effect of the variability of the hydraulic conductivity (Mantoglou and Gelhar, 1987b, 1987c).

Substituting Eq. [33] and [35] into Eq. [34] yields an upscaled flow equation of a similar form as the local flow equation:

[36]

Main Results
Although the form of the upscaled flow equation is equal to the local one, the spatial variability of the local hydraulic properties leads to some important differences between the averaged or upscaled flow process and the local one. The variances of hydraulic parameters and the pressure head as well as their covariances in Eq. [32] and [33] show that the effective hydraulic conductivity is determined by the spatial variability of the local hydraulic parameters. The variance of the pressure heads and the covariances between pressure heads and hydraulic parameters, which are obtained by solving the moment equations (Eqs. [29–31] and the expected water flux Q_i, depend in a nonlinear way on the expected hydraulic gradient J. For transient flow conditions, they also depend on the change of pressure head with time. As a consequence, the effective hydraulic conductivity is a function of the hydraulic gradient and of its history, which implies a hysteretic behavior of the larger scale system.

Yeh et al. (1985a, 1985b, 1985c) studied steady-state unsaturated flow in soils with an isotropic structure and in stratified soils for vertical mean infiltration through an unbounded heterogeneous flow domain. The resulting stochastic flow equation was solved using spectral representation techniques. Two key results were obtained from the theoretical analysis. First, the variance of the mean capillary pressure increases monotonically when and K_s (Eq. [7]) are independent. For a negative correlation between and K_s, which represents a Miller–Miller similar medium, the variance of the mean pressure head is a non-monotonic function of the degree of saturation or the mean capillary pressure. Second, spatial heterogeneity of local hydraulic properties induces an anisotropy of the effective hydraulic conductivity. This is nicely illustrated by the work of Yeh et al. (1985c) using data from a sandy loam and a silty clay loam (Fig. 3 ). The anisotropy ratio of the effective unsaturated hydraulic conductivity calculated from stochastic perturbation theory depends strongly on the mean capillary pressure and increases when the soil gets drier. The direct averaging method leads to a higher anisotropy than in the stochastic theory. This method takes the ratio of arithmetic and harmonic mean conductivity to estimate the effective conductivity anisotropy (for details, see Yeh et al., 1985a). In their derivation for the cases where K_s(x) and (x) are spatially variable, both parameters were assumed to have the same correlation scale and covariance function. They verified their theoretical findings using the field data of Nielsen et al. (1973) obtained on a 150-ha field site and various laboratory experiments. Yeh et al. (1985c) concluded that a major limitation in applying the stochastic method to upscaling lies in the determination of the statistical properties of K_s(x) and (x).

View larger version (26K): [in this window] [in a new window]   Fig. 3. Comparison of the anisotropic ratios obtained from the stochastic theory and the direct average methods for Maddock sandy loam and for Panoche silty clay loam; H is the mean capillary pressure head (from Yeh et al., 1985c; reproduced with permission of the American Geophysical Union).

Mantoglou and Gelhar (1987a, 1987b, 1987c) extended the analysis of Yeh et al. (1985a, 1985b) to transient unsaturated flow assuming that the hydraulic parameter fields are a realization of three-dimensional cross-correlated, second-order, stationary, statistically anisotropic, hydraulic parameter fields. Their analysis showed that the effective hydraulic conductivity (K_eff) and the effective differential moisture capacity (C_eff) depend on the mean values of lnK_s(x), (x), and C(,x), the stochastic properties of the soil properties and the mean flow characteristics such as mean pressure head H, its gradients, and change with time, H/t. The dependence on the mean pressure head gradient and H/t suggests hysteresis of the effective parameters, which was shown to increase with increasing heterogeneity of local-scale parameters. Also, the variances of the local pressure heads and fluxes were found to depend on the larger scale gradients and mean flow conditions (wetting vs. drying; see also Yeh et al., 1985a, 1985b). In addition they confirmed that the anisotropy of the effective hydraulic conductivity depends on the mean soil moisture content.

Harter and Zhang (1999) studied steady-state water flow in heterogeneous unsaturated media and examined the effect of variance lnK_s, ln(), the mean soil water tension, the tortuosity m_g in the Gardner–Russo model (Eqs. [7] and [15]), and soil layering on the coefficient of variation of soil water content. Soil layering was introduced by considering lnK_s to be anisotropic with smaller correlation length in the vertical than the horizontal direction. They showed that soil water content variability increases with soil water tension, as was also found by Yeh (1989) for the case of uncorrelated Gardner parameters. Under dry conditions, the soil water content was found to be very sensitive to the value of ln(), whereas it was found to be sensitive to the value of the tortuosity in the wet range. At low moisture content, small changes in the variability of ln() with respect to _lnKs can lead to substantial changes in soil moisture content variability expressed as a coefficient of variation.

The results obtained in the former studies are based on spectral methods that assume unbounded media in which the mean flow variables are constant or change slowly. This implies that the length scale of the flow domain has to be much larger than the scale of heterogeneity of the relevant hydraulic properties. Russo (1992) relaxed these conditions by proposing an upscaling procedure that combined the stochastic approach of Yeh et al. (1985a) with the use of effective block properties (e.g., unsaturated hydraulic conductivity) having a scale comparable with the scale of the formation heterogeneity defined by, e.g., correlation functions. Mantoglou (1992) extended the three-dimensional stochastic perturbation approach to the study of transient nonstationary flow in bounded domains. Zhang and Winter (1998), Zhang (1999), and Zhang and Lu (2002) used numerical solutions of upscaled moment equations and derived the variance of the pressure heads and flow velocities for nonstationary and nonsteady flow fields.

Experimental Verification of Stochastic Theories
There are only a few studies in the literature dealing with the experimental verification of stochastic theories for soil water processes. Important studies were conducted by McCord et al. (1991), Jensen and Mantoglou (1992), and Wendroth et al. (1999) at the field scale, and Wildenschild and Jensen (1999) and Ursino et al. (2001) at the laboratory scale. McCord et al. (1991) performed hillslope tracer experiments. They derived the anisotropy ratio from simultaneous measurements of the mean hydraulic gradient and the mean water flux, which was estimated from the path of the center of mass of an injected tracer plume. Their analysis illustrated that there was a considerable anisotropy of the hydraulic conductivity, which was caused by layered anisotropy. Simulations of water flow using a state-dependent anisotropy of the hydraulic conductivity reproduced the observed tracer movement the best. Ursino et al. (2001) performed unsaturated tracer experiments in a sand tank that was randomly filled with blocks of 5 by 5 by 0.5 cm of three different sand types that were rotated at an angle of 45° to the horizontal. The tank was uniformly irrigated at the top and tracer experiments were performed for three different flow rates. The deviation of the trajectory of