# Spatial variability of soil surface properties and consequences for the annual and monthly water balance of a semiarid environment (EFEDA experiment)

Isabelle Braud, R. Haverkamp, J.L. Arrue, M.V. López

► To cite this version:

Isabelle Braud, R. Haverkamp, J.L. Arrue, M.V. López. Spatial variability of soil surface properties and consequences for the annual and monthly water balance of a semiarid environment (EFEDA experiment). Journal of Hydrometeorology, 2003, 4 (1), pp.121-137. hal-02581329

HAL Id: hal-02581329

https://hal.inrae.fr/hal-02581329v1

Submitted on 9 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Spatial Variability of Soil Surface Properties and Consequences for the Annual and Monthly Water Balance of a Semiarid Environment (EFEDA Experiment)

ISABELLE BRAUD

LTHE, Grenoble, and CEMAGREF, UR Hydrologie-Hydraulique, Lyon, France

RANDEL HAVERKAMP

LTHE (CNRS UMR 5564, INPG, UJF, ORSTOM), Grenoble, France

J. L. ARRÚE AND M. V. LÓPEZ

Consejo Superior de Investigaciones Científicas, Estación Experimental de Aula Dei, Zaragoza, Spain

(Manuscript received 25 October 2001, in final form 28 June 2002)

## ABSTRACT

During the second phase of the European International Project on Climatic and Hydrological Interactions between Vegetation, Atmosphere, and Land Surface (ECHIDA) Field Experiment in a Desertification Threatened Area (EFEDA) the spatial variability of the soil water retention and hydraulic conductivity characteristics of layers at 2–12- and 17–27-cm depth was characterized. A simplified method, based on particle size distribution and simple infiltration tests, was used. It provided these characteristics at the nodes of a 1-km grid over 10 × 10 km² around the town of Tomelloso (Castilla–La Mancha, Spain).

A total number of 78 sample points were used to address the problem of soil surface properties variability and its consequences on the monthly and annual water balance. The Simple Soil Plant Atmosphere Transfer model (SiSPAT) 1D Soil–Vegetation–Atmosphere Transfer (SVAT) model was run with a 1-yr climatic forcing for the 78 soil profiles until equilibrium was reached. As no runoff was generated, the spatial variability of the water budget components only concerned soil evaporation, transpiration, and deep drainage. It was found that (i) the choice of the type of boundary condition at the bottom of the soil profile was greatly influencing the final variability, (ii) the variability of transpiration was the largest in situations of water stress for the vegetation, and (iii) soil evaporation was the most sensitive component when plants were well supplied with water.

Various aggregation methods of soil surface parameters (use of the arithmetic mean, median of the parameters, or parameters associated to the average soil texture of the Clapp and Hornberger classification) were assessed. The use of median parameters in a single 1D simulation was found to provide the best agreement with the average of the 78 simulations performed for each grid cell using locally measured soil properties. The use of average soil texture parameters led to a significant bias, especially in the case of water stress.

## 1. Introduction

In atmospheric and general circulation models, surface and soil hydraulic properties are often assumed to be homogeneous over meshes of 100–10 000 km². However, these characteristics and, in particular, soil hydraulic characteristics are known to be highly variable in space (e.g., Nielsen et al. 1973; Russo and Bressler 1981; Vauclin et al. 1994; Mallants et al. 1996; Lewan 1996; Bell et al. 1980; Haverkamp et al. 1996). The practice of assuming homogeneous properties can lead to errors in the calculation of surface fluxes. As these

fluxes constitute the lower boundary condition of most large-scale atmospheric and/or climatic models, the prediction reliability of these models is strongly reduced (Shao and Henderson-Sellers 1996; Henderson-Sellers 1996, see special issue of *Global Planetary Change*, vol. 13, no. 1–4). The introduction of some of the variability aspects into general circulation models (GCMs) using statistical dynamic approaches has been shown to improve significantly their performance (e.g., Entekhabi and Eagleson 1989; Famiglietti and Wood 1991, 1992; Avissar 1992, 1995). A large effort was also dedicated to the definition of “effective” or “aggregated” parameters (e.g., Raupach and Finnigan 1995), supposedly to provide the same mean flux as that obtained by resolving explicitly the spatial variability. When soil moisture was fairly homogeneous, Noilhan and Lacarrère (1995) ob-

<footnote>*Corresponding author address:* Isabelle Braud, CEMAGREF, UR Hydrologie-Hydraulique, CP 220, 3bis Quai Chauveau, 69336 Lyon Cédex 9, France.
E-mail: braud@lyon.cemagref.fr</footnote>
tained a reasonable agreement between average evaporation fluxes derived from a 3D atmospheric model and the equivalent 1D simulation. When variability of soil hydraulic properties was large, some discrepancies could be observed between average fluxes and those derived using average soil parameters (e.g., Braud et al. 1995b; Kim and Stricker 1996; Kabat et al. 1997; Braud 1998; Boulet et al. 1999). Milly and Eagleson (1987) and Kim et al. (1997) showed that the effect of soil parameter variability was largest in the case of runoff generation and that, in general, effective parameters could be more easily defined for total evaporation.

In this paper the regional scale (10 × 10 km²) relevant for atmospheric models is considered. First, a dataset collected in central Spain in Castilla-La Mancha is presented. Water retention and hydraulic conductivity parameters were measured at the nodes of a 1-km grid over 100 km². For the estimation of soil hydraulic properties, existing methods can be categorized as being either predictive or based on direct experimental measurement techniques. Observations can be made on samples in the laboratory or in situ at the location of interest. Methods based on direct observations are often difficult to implement and time consuming. Predictive methods employ information on textural and structural properties such as particle size distributions, organic matter, and/ or dry bulk density in order to estimate the hydraulic properties (e.g., Clapp and Hornberger 1978). The dataset presented here uses a simplified in situ method, aiming at minimizing time and human resources needed to estimate soil hydraulic properties when the number of samples is large. In the following, the soil water retention and the hydraulic conductivity curves are represented by the Brooks and Corey (1964) model:

$$ 
\left\{ 
\begin{array}{ll}
\frac{\theta}{\theta_s} = \left( \frac{h_{bc}}{h} \right)^{\lambda} & \text{for } h \le h_{bc} \\
\theta = \theta_s & \text{for } h_{bc} \le h \le 0
\end{array}
\right. 
\qquad (1) 
$$

$$ \frac{K}{K_s} = \left( \frac{\theta}{\theta_s} \right)^{\eta}, \qquad (2) $$

where h (m) is the soil water pressure, θ (m³ m⁻³) is the volumetric moisture content, and K (m s⁻¹) is the soil hydraulic conductivity. At each location, there are five unknown parameters: two shape parameters λ (−) and η (−) defined by the soil textural properties; and three scale parameters, that is, the saturated water content θs (m³ m⁻³), the saturated hydraulic conductivity Ks (m s⁻¹), and the scale parameter for the pressure hbc (m), all three strongly related to the soil structural properties (Haverkamp et al. 1998a). The residual water content is assumed to be 0.

The dataset was used to study at regional scale the variability of the components of the annual and monthly water budget of a vineyard representing the typical vegetation of the area. Several choices for the lower bound-

ary condition of the soil profiles were considered in the soil-vegetation-atmosphere transfer model (SVAT) used for this study. Different rules for aggregating surface soil parameters were investigated in order to analyze the feasibility of an effective lumped parameterization.

## 2. Materials and methods

### a. The study area

The field measurement program was conducted in the framework of the European International Project on Climatic and Hydrological Interactions between Vegetation, Atmosphere, and Land Surface (ECHIVAL) program. This program was dedicated to the study of the interaction between the soil, the vegetation, and the atmosphere in various regions of the globe. One of these experiments, the ECHIVAL Field Experiment in a Desertification Threatened Area (EFEDA), focused on the semiarid environment (Bolle et al. 1993). The study area was situated near Tomelloso, in the Castilla-La Mancha region located in central Spain (Fig. 1). An area of 10 × 10 km² (39°09'36"N, 2°58'03"W, for the top left corner, and 39°04'42"N, 2°51'46"W, for the bottom right corner) was selected. The area represented a typical grid element within a GCM. Determination of aggregation/disaggregation rules for surface parameters and fluxes was the main focus of the experimental and modeling work. The soil of the study area was fairly uniform. However, a field campaign was conducted in June 1994 in order to document the spatial variability of soil hydraulic properties over the whole area. The latter was covered with a 1-km sampling grid mesh (Fig. 1) resulting in 100 measurement locations. The altitude varies from 669 m above sea level at grid point site A0 (Fig. 1) to 835 m at site G9, increasing in the northwest-southeast (NW-SE) direction.

The average slope of the area was less than 2%. The area covered mostly cultivated land (92 sites) with the following land use distribution: vineyards (49 sites); fallow land, generally part of a cereal-fallow rotation scheme (23 sites); barley (11 sites); wheat (5 sites); olive trees (3 sites); and chickpeas (1 site). Only eight sites were lying at abandoned vineyards (five sites) and shrubs (three sites). Soil profiles were composed of an arable top layer with a depth of approximately 30 cm overlying a deep calcareous crust layer down to 4–5 m. Hydraulic properties of the crust layer were characterized using an internal drainage method at site TOM 6 (E4, Fig. 1). Details are given by Haverkamp et al. (1996). A significant proportion of stones and limestone porous crust fragments were encountered in the arable top layer. On the average 35% of the soil surface was covered by stones and rock fragments of 5–10 cm in diameter.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Location of the study area and map of the silt + clay content (%) of the study area located in the Tomelloso region (Castilla-La Mancha, Spain). The sampling grid is also pointed. Full line gives the altitude of the area.</figcaption></figure>

b. Derivation of water retention and hydraulic conductivity parameters

1) SHAPE PARAMETER λ OF THE RETENTION CURVE

The particle size distribution can be modeled using the following function (Haverkamp et al. 1998a):

$$F(d) = \left[ 1 + \left( \frac{d_g}{d} \right)^N \right]^M \quad \text{with} \quad N \frac{2}{1-M} \quad (3)$$

where d (m) is the particle size diameter, d_g (m) is the particle size scale parameter, and M (-) is the particle size shape parameter.

When such data are available, M can be directly identified from particle size distribution data. If only the soil texture is available (i.e., clay, silt, and sand content), Zammit (1999) proposed a “cartography” of this parameter within the soil texture triangle based on the analysis of a 750 soil database called Grizzly (Haverkamp et al. 1998c).

The shape parameters λ of the water retention curve can be derived from the knowledge of soil texture and, more specifically, from the knowledge of the values for M and N. The following relationship was fitted by Braud et al. (2002, manuscript submitted to Eur. J. Soil Sci., hereafter BRAU) on data from the Grizzly soil database (Haverkamp et al. 1998c):

$$\lambda = MN(0.382 + 0.0469MN); \quad (4)$$

this equation was used in the present study for the estimation of parameter λ.

2) SATURATED VOLUMETRIC WATER CONTENT θ_s

The soil porosity ε (-) can be determined from the soil dry bulk density ρ_d (g cm⁻³) and soil particle density ρ_s (ρ_s = 2.65 g cm⁻³) through

$$\varepsilon = 1 - \frac{\rho_d}{\rho_s}. \quad (5)$$

The saturated soil volumetric water content θ_s (m³ m⁻³) can be related to the porosity using the following relationship proposed by Haverkamp et al. (1998c):

$$\theta_s = \varepsilon 2^{(\lambda/(2+\lambda)-M)}. \quad (6)$$

3) SHAPE PARAMETER h_bc OF THE WATER RETENTION CURVE

Zammit (1999) showed that predictive methods based on textural analysis failed to predict this parameter. The latter is more related to soil structure and therefore to local effects (Haverkamp et al. 1998a). Consequently, the most reliable procedure to determine h_bc seemed to be through direct measurement of matric potential and water content or indirectly from simple infiltration tests (BRAU). Infiltration is very sensitive to this parameter because h_bc dominates the retention behavior in the vicinity of saturation. It seems therefore advisable to make these measurements near saturation if only one field campaign can be conducted. If several campaigns can be realized different parts of the retention curve must be covered. Least square errors optimization techniques
on the data pairs (matric potential, volumetric water content) available can subsequently be used to derive $h_{bc}$, assuming that the other parameters are known.

### 4) SHAPE PARAMETER $\eta$ OF THE HYDRAULIC CONDUCTIVITY CURVE

The derivation of the shape parameter for the hydraulic conductivity curve $\eta$ can be performed by introducing a tortuosity factor $p$, given by

$$M = \frac{\lambda}{2 + \lambda}(1 + p), \qquad (7)$$

and related to $\eta$ by (Haverkamp et al. 1998a)

$$\eta = \frac{2}{\lambda} + 2 + p. \qquad (8)$$

### 5) SATURATED HYDRAULIC CONDUCTIVITY $K_s$

Like the scale parameter $h_{bc}$ models based on texture fail to predict the structure-related scale parameter $K_s$ (Zammit 1999). The most reliable way to get an estimation of this parameter was thus to perform in situ experiments. One of the simple ways, which does not require too much time, consists of artificially wetting the soil under a constant positive pressure head. Then the Green and Ampt (1911) approach can be used to characterize the infiltration process (Hillel 1980). Assuming uniform initial volumetric soil water content, $\theta_0$ (m$^3$m$^{-3}$), the Green and Ampt model leads to the following relationship between the vertical cumulative infiltration $I_{ID}$ (m) and the corresponding infiltration time $t_{inf}$ (s):

$$I_{ID} = K_s t_{inf} + (\theta_s - \theta_0)(h_{surf} - h_f) \times \log\left[1 + \frac{I_{ID}}{(\theta_s - \theta_0)(h_{surf} - h_f)}\right], \quad (9)$$

where $h_f$ (m) is the wetting front suction, and $h_{surf}$ (m) is the constant head at the soil surface. The saturated hydraulic conductivity can be estimated from this expression. The suction at the wetting front can be calculated by (10). Indeed, under positive head infiltration, the sorptivity is an intrinsic property of the soil and its value should not depend on the model chosen to represent the soil water retention and hydraulic conductivity curves. Therefore, expressions obtained using the Brooks and Corey model and the Green and Ampt models must be equal, leading to Eq. (10) (Haverkamp et al. 1998b; BRAU):

$$h_f = \frac{2\lambda\eta(\lambda\eta - 1) + \lambda(2\lambda\eta - 1)h_{bc}}{(\lambda\eta - 1)(\lambda\eta - 1 + \lambda)} \qquad (10)$$

### c. Experimental work and derivation of the water retention and hydraulic conductivity curves

At each location over the sampling grid, measurements were carried out at two depths: 2–12-cm depth

for layer 1 and 17–27-cm depth for layer 2, not situated on the same profile but at two adjacent spots. In June 1994, when the experiment was conducted, no rain had fallen for several months. Therefore, the surface layer was very dry and an artificial wetting was required. A constant head water supply device (“bucket method”) consisting of a single-ring infiltrometer (a bottom-free plastic bucket of 18-cm diameter) and a 2-L plastic bottle acting as a Mariotte flask were used. Appropriate wetting (pressure range and wetting front depth) was obtained after infiltration of 1–2.25 L of water.

After drainage for a known time, the volumetric water content $\theta_{TDR}$ was determined through time domain reflectrometry (TDR; Topp et al. 1980) using a Tektronix 1052 C cable tester. The soil matric potential $h$ was measured simultaneously using a handheld micropressure transducer tensiometer developed by the Division of Soils, Commonwealth Scientific and Industrial Research Organization (CSIRO), Townsville, Australia. Besides these in situ measurements of volumetric water content and water pressure, soil samples were collected at the same depths for dry bulk density and gravimetric water content $W_{ap}$ determination. An auger with a volume of about 500 cm$^3$ was used to sample approximately the same soil volume as that sensed by the TDR probe.

A team of five people accomplished the fieldwork in two weeks (seven sites per day). Average time taken per soil site was 67 min (20 min for location of site and transportation, 24 min for water infiltration, and 23 min for actual soil measurements and sampling). This makes this type of soil characterization campaign accessible even for large areas, because the human and time investment is relatively low.

The collected soil samples were also used for particle size distribution analysis, combining sieving and sedimentation techniques for the determination of the following particle size fractions: very coarse sand (2000–1000 µm), coarse sand (1000–500 µm), medium sand (500–200 µm), fine sand (200–100 µm), very fine sand (100–50 µm), coarse silt (50–20 µm), medium-fine silt (20–2 µm), and clay (<2 µm) (Gee and Bauder 1986). The percentage of stones and crust fragments in the coarse soil fraction (>2000 µm) was determined before the particle size analyses.

It has been mentioned that a large proportion of stones and crust fragments were present in the soil samples. Furthermore, it was observed that the crust fragments were also holding water. Therefore, the dry bulk density derived from $\theta_{TDR}/W_{ap}$ was only an apparent value and was different from the fine soil bulk density $\rho_d$, needed in Eq. (5) used to derive the porosity. This bias was obvious given the substantial number of sites having large values of apparent dry bulk density (higher than 2 g cm$^{-3}$ in some cases). A simple correction (R. Haverkamp 1996, unpublished manuscript) was thus developed to obtain first the fine soil volumetric water content $\theta$, obviously lower than $\theta_{TDR}$, and second the
fine soil dry bulk density $\rho_d$. Data discussed below will be the corrected values.

Parameters of Eq. (3) were fitted on the experimental fine soil particle size distribution functions. Then the shape parameter $\lambda$ was deduced from Eq. (4) and introduced in Eq. (6) to determine the saturated water content $\theta_s$. Once $\lambda$ and $\theta_s$ were known for each layer, the scale parameter $h_{bc}$ was calculated by introducing the only data pair $(h, \theta)$ available in Eq. (1). It was not possible to use the two layers to optimize the scale parameter, because the saturated water content of both layers were found to be statistically different (Haverkamp et al. 1996); and each layer had to be treated separately, which of course, reduces the robustness of the procedure.

For the hydraulic conductivity curve at each location, the shape parameter $\eta$ was evaluated using Eqs. (7) and (8). The saturated hydraulic conductivity was derived from (9), assuming a value of the surface suction $h_{surf}$ of 1 cm.

The initial soil moisture content $\theta_0$ was not available. It was assumed that, given the long time since the last rainfall, steady state over the whole area had been reached. A constant arbitrary value $h_0 = 10^7$ cm, representative of dry surface conditions, was thus assumed at each point and the corresponding water content calculated using (1). The $\theta_0$ values ranged between 2% and 10%. It was checked that the results on $K_s$ were not very sensitive to the choice of this value for $h_0$. Furthermore, the measured infiltrated height was representative of the 3D infiltration, whereas the 1D value must be used in (9). The following correction was applied (Haverkamp et al. 1994):

$$I_{1D} = I_{3D} - \frac{\gamma}{r_d} 2(h_{surf} - h_f) K_s t_{inf}, \quad (11)$$

where $\gamma = 0.7$ and $r_d$ is the radius of the infiltrometer. This leads to the following equation, solved iteratively to obtain $K_s$, once $h_f$ had been estimated using (10):

$$K_s t_{inf} = \frac{I_{3D} - (\theta_s - \theta_0)(h_{surf} - h_f) \log \left[ 1 + \frac{I_{3D}}{(\theta_s - \theta_0)(h_{surf} - h_f)} - \frac{\gamma}{r_d} \frac{2K_s t_{inf}}{(\theta_s - \theta_0)} \right]}{1 + \frac{2\gamma(h_{surf} - h_f)}{r_d}}. \quad (12)$$

The experimental protocol used during the EFEDA experiment and reported in this paper was a first attempt to develop such a large-scale strategy for soil hydraulic properties derivation. Since this early work, the method has been refined and theoretically justified. An update of the procedure, known as the Beerkan method, is provided in BRAU and could be used in practice. Figure 2 shows the retention curves and hydraulic conductivity curves drawn using (1) and (2) and the values as calculated from the method described above for all the measurement points. A large scatter can be observed showing the variability of those curves at this scale. Table 1 summarizes the statistical properties of the different parameters. The spatial correlation of the data was also examined through variogram analysis. No spatial correlation was detected at the scale of the measurements (1 km), justifying therefore the assumption of statistical independence of the soil profiles used in the subsequent modeling work. At the field scale (4 ha), Vandervaere (1995) found correlation lengths of about 150 m for soil dry bulk density and hydraulic conductivity. The coefficients of variation (8%–10% for porosity and 20%–30% for hydraulic conductivity) were much smaller than those reported in Table 1 calculated at the regional scale.

### d. The SiSPAT model

An extensive description of the Simple Soil Plant Atmosphere Transfer (SiSPAT) model can be found in

Braud et al. (1995c) and Braud (2000) and only a brief summary is presented here. SiSPAT is a vertical 1D model. The driving forces are climatic time series of air temperature and humidity, wind speed, incoming solar and longwave radiation, and rainfall. In the soil, coupled heat and mass transfer equations, derived from the Richards (1931) equation, are solved for temperature $T$ and matric potential $h$. They include both liquid and vapor transfers as formulated by Philip and De Vries (1957) or Milly (1982). The model deals with vertically heterogeneous soils. The upper boundary conditions are obtained by the solution of the soil–plant–atmosphere interface, modeled as a two-source model of heat and vapor, where bare soil and vegetation are considered separately (Deardorff 1978), as formulated by Shuttleworth and Wallace (1985) and Taconet et al. (1986). Five equations can be written: energy budget over bare soil and vegetation, continuity of the sensible and latent heat fluxes through the canopy, and continuity of the surface flux at the soil surface. Leaf temperature $T_v$, canopy temperature $T_{av}$, canopy specific humidity $q_{av}$, soil surface temperature $T_1$, and surface matric potential $h_1$ can thus be calculated. The resolution of this module provides the surface soil heat and mass fluxes and the surface matric potential $h_1$ and temperature $T_1$. In the case study, surface matric potential and temperature were used as the upper boundary condition of the soil module for all the simulations because it proved to be numerically more stable. When rainfall exceeds the in-
<figure><img src="image_3.png" /><figcaption>FIG. 2. (left) Retention curves of the 78 soil samples; (right) hydraulic conductivity curves of the 78 soil profiles for layer 1 (2–12 cm). Points are the measurements made at each location.</figcaption></figure>

filtration capacity of the soil, saturation of the surface occurs. The matric potential at the surface is set to 0 and the runoff is calculated from the mass budget equation.

The incoming energy is partitioned between bare soil and vegetation through a shielding factor $σ_f$ (Deardorff 1978; Taconet et al. 1986). In the soil, a root extraction model, based on Federer (1979) is included. For each soil layer, a soil-root and a root-leaf resistance are put in series. The moisture extraction in layer *j* is proportional to the water potential difference between the leaf $h_j$ and the soil $h_s$. The leaf water potential is calculated by assuming steady state at each time step and that total moisture extraction is equal to the transpiration calculated from the atmospheric conditions. The leaf water potential controls the water stress function of the stomatal resistance, which also depends on the incoming radiation and vapor pressure deficit. Iterative procedures are used to solve the various modules of the model. They are described in detail in Braud (2000).

### e. Method used to assess surface fluxes variability and effective parameters soundness, in response to surface soil parameters variability

At the end of the data processing, 78 sample points were available with estimation of both the retention and hydraulic conductivity curves for layers 2–12 and 17–27-cm depth. Some 22 sample points had to be removed from the analysis because the infiltrated volume was not recorded at the beginning of the field work and the

saturated hydraulic conductivity could not be determined for these points. These data were used to investigate the influence of soil surface hydraulic properties spatial variability on surface fluxes, using the SiSPAT SVAT model. The methodology was as follows.

1. One year of half-hourly values of climate forcing (air temperature, humidity, and wind speed at 2 m; solar and longwave radiation, and rainfall) typical of semiarid Spain had been generated by the 1D version of the U.K. Meteorological Office general circulation model (Lean 1992) in the framework of the Spatial Variability of Land Surface Processes II (SLAPS II; Dooge et al. 1994). This dataset was used as the forcing of the SiSPAT SVAT model. Daily values are presented in Fig. 3. The rainfall and the other climate forcing variables were assumed to be homogeneous over the whole area in order to focus the study on soil surface properties influence.

2. For the same reason, one vegetation type (vineyard) typical of the area of Tomelloso was chosen. On the 100 km² studied area, vineyards only covered 43% of the surface but Sene (1996) reported values of 70%-80% for the whole EFEDA zone. Plant parameters necessary to run the SiSPAT model had been calibrated by Braud et al. (1995a). The use of these parameters, including plant resistance, maximum root density, and parameters describing the stomatal conductance, allowed for the reproduction of one month of observations of surface fluxes measured in 1991 during the first EFEDA experiment. The root-
<table><caption>TABLE 1. Basic statistics of secondary soil hydraulic parameters for layer 1 (2–12 cm) and layer 2 (17–27 cm): min, max, mean, std dev, coefficient of variation CV, number of observations n<sub>obs</sub>, type of distribution tested with the χ² test (N = normal, LN = lognormal), and value of the probability p test of the χ² test. Values with * mean rejection at the 10% level; ** means rejection at the 5% level. Values in parentheses are the corresponding values for the logarithm of the variable.</caption><thead><tr><th></th><th>Layer</th><th>Min</th><th>Max</th><th>Mean</th><th>Std dev</th><th>CV (%)</th><th>n<sub>obs</sub></th><th>Type of distribution</th><th>p test value of the χ² test</th></tr></thead><tbody><tr><td>Dry bulk density ρ<sub>d</sub> (g cm<sup>-3</sup>)</td><td>1</td><td>0.94</td><td>1.56</td><td>1.27</td><td>0.16</td><td>12.5</td><td>100</td><td>N</td><td>0.10</td></tr><tr><td></td><td>2</td><td>0.86</td><td>1.69</td><td>1.32</td><td>0.20</td><td>14.9</td><td>98</td><td>N</td><td>0.05*</td></tr><tr><td>Soil porosity ε</td><td>1</td><td>0.41</td><td>0.65</td><td>0.52</td><td>0.06</td><td>11.5</td><td>100</td><td>N</td><td>0.10</td></tr><tr><td></td><td>2</td><td>0.36</td><td>0.68</td><td>0.50</td><td>0.07</td><td>14.8</td><td>98</td><td>N</td><td>0.05*</td></tr><tr><td>Saturated water content θ<sub>s</sub> (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>1</td><td>0.39</td><td>0.63</td><td>0.50</td><td>0.06</td><td>11.6</td><td>100</td><td>N</td><td>0.19</td></tr><tr><td></td><td>2</td><td>0.34</td><td>0.66</td><td>0.48</td><td>0.07</td><td>15.1</td><td>98</td><td>N</td><td>0.08*</td></tr><tr><td>Shape parameter of the particle size distribution M</td><td>1</td><td>0.102</td><td>0.196</td><td>0.138</td><td>0.018</td><td>13.0</td><td>100</td><td>N</td><td>0.47</td></tr><tr><td></td><td>2</td><td>0.100</td><td>0.194</td><td>0.136</td><td>0.018</td><td>13.2</td><td>98</td><td>N</td><td>0.44</td></tr><tr><td>Shape parameter of the Brooks and Corey retention model λ</td><td>1</td><td>0.119</td><td>0.243</td><td>0.165</td><td>0.026</td><td>15.7</td><td>100</td><td>N</td><td>0.51</td></tr><tr><td></td><td>2</td><td>0.113</td><td>0.269</td><td>0.165</td><td>0.027</td><td>16.5</td><td>98</td><td>N</td><td>0.42</td></tr><tr><td>Shape parameter of the Brooks and Corey model η</td><td>1</td><td>10.15</td><td>15.60</td><td>12.88</td><td>1.19</td><td>9.2</td><td>100</td><td>N</td><td>0.85</td></tr><tr><td></td><td>2</td><td>9.61</td><td>16.13</td><td>12.90</td><td>1.27</td><td>9.8</td><td>98</td><td>N</td><td>0.06*</td></tr><tr><td>Scale parameter of the particle size distribution d<sub>s</sub> (mm)</td><td>1</td><td>177.6</td><td>1774.0</td><td>629.9</td><td>272.3</td><td>43.2</td><td>100</td><td>LN</td><td>0.03**</td></tr><tr><td></td><td>2</td><td>160.8</td><td>1526.0</td><td>660.7</td><td>287.8</td><td>43.5</td><td>98</td><td>LN</td><td>0.28</td></tr><tr><td>Scale parameter of the Brooks and Corey retention model -h<sub>k</sub> (cm)</td><td>1</td><td>0.52</td><td>52.8</td><td>12.4</td><td>10.7</td><td>86.5</td><td>99</td><td>LN</td><td>0.36</td></tr><tr><td></td><td>2</td><td>(-0.66)</td><td>(3.97)</td><td>(2.12)</td><td>(0.96)</td><td>(45.1)</td><td>(99)</td><td>LN</td><td></td></tr><tr><td></td><td></td><td>0.19</td><td>76.8</td><td>20.0</td><td>17.7</td><td>88.8</td><td>93</td><td>LN</td><td>0.02**</td></tr><tr><td></td><td></td><td>(-1.67)</td><td>(4.34)</td><td>(2.46)</td><td>(1.26)</td><td>(51.3)</td><td>(93)</td><td>LN</td><td></td></tr><tr><td>Saturated hydraulic conductivity, K<sub>s</sub> (cm s<sup>-1</sup>)</td><td>1</td><td>0.00014</td><td>0.0075</td><td>0.00151</td><td>0.00129</td><td>85.5</td><td>83</td><td>LN</td><td>0.51</td></tr><tr><td></td><td></td><td>(-8.84)</td><td>(-4.89)</td><td>(-6.84)</td><td>(0.883)</td><td>(12.9)</td><td>(83)</td><td>LN</td><td></td></tr><tr><td></td><td>2</td><td>0.000056</td><td>0.0083</td><td>0.00122</td><td>0.00149</td><td>122.3</td><td>79</td><td>LN</td><td>0.23</td></tr><tr><td></td><td></td><td>(-9.79)</td><td>(-4.78)</td><td>(-7.28)</td><td>(1.102)</td><td>(15.1)</td><td>(79)</td><td>LN</td><td></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Daily value of the atmospheric forcing: (top) air temperature (°C), specific humidity (10⁻³ kg kg⁻¹), wind speed (m s⁻¹); (bottom) incoming solar and longwave radiation (W m⁻²), rainfall (mm⁻¹), and LAI (10⁻³ LAI) (—).</figcaption></figure>

mean-square error of latent and sensible heat fluxes was 37 and 25 W m⁻², respectively. These figures lie in the range of errors reported in an intercomparison study of various surface schemes (Linder et al. 1995) performed using the same dataset. A drill showed that the vineyard roots could reach a depth of 3 m (Santa Olalla Manas 1994); a soil profile of 4 m was chosen with a rooting depth of 3 m. A plausible 1-yr time series of leaf area index (LAI) was defined for a vineyard cultivation, based on information given in Sene (1996). It also appears in Fig. 3.

3) It was assumed that the whole area could be represented by independent soil profiles, with vertical heat and water transfers. This hypothesis was reasonable as the water table was located below 10-m depth (Sene 1996) and no spatial correlation between parameters was found at the scale of the study. Such an hypothesis was also retained in the work of Peck et al. (1977). Some 78 soil profiles were defined, and the water retention and hydraulic conductivity curves

measured at the grid points were assigned to one of them for layers 0–20 and 20–50 cm. No data about the variability of soil properties below the 17–27-cm layer were available. Therefore, the properties for the 50–400-cm layers were kept identical for all soil profiles and assigned a value measured for the calcite crust encountered at that depth (Haverkamp et al. 1996). The saturated hydraulic conductivity of that layer was measured using suction infiltrometers leading to $K_s = 2.78 \times 10^{-6} \text{ m s}^{-1}$ (J. Vandervaere 1996, personal communication). This value was consistent with a value of $K_v = 4.63 \times 10^{-6} \text{ m s}^{-1}$ used by Sene (1996) in a water balance study performed in the same region. Another trial was also performed with a value of the saturated hydraulic conductivity divided by 10.

4) For the 78 soil profiles, the same atmospheric forcing was applied during two or three years in order to obtain equilibrium of the solution. In analogy with climate modeling, the word “equilibrium” means that the soil water storage at the beginning of the
<table><caption>TABLE 2. Components of the mean water budget (in mm) for the two/three years of simulation and the three lower boundary condition choices. Deep drainage was positive for percolation and negative for capillary rises. Total rainfall was 352 mm in all the cases and runoff was 0. Evaporation of intercepted rainfall was not given in this table but was of the order of 7 mm in all cases.</caption><thead><tr><th></th><th>Total evaporation</th><th>Deep drainage</th><th>Change in water storage</th><th>Bare soil evaporation</th><th>Transpiration</th></tr></thead><tbody><tr><td>First-year constant h</td><td>574</td><td>-45</td><td>-177</td><td>336</td><td>231</td></tr><tr><td>Second-year constant h</td><td>549</td><td>-192</td><td>-5</td><td>319</td><td>223</td></tr><tr><td>First-year gravitational K<sub>s</sub> = 2.78 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>549</td><td>146</td><td>-343</td><td>334</td><td>208</td></tr><tr><td>Second-year gravitational K<sub>s</sub> = 2.78 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>410</td><td>6</td><td>-64*</td><td>303</td><td>97</td></tr><tr><td>First-year gravitational K<sub>s</sub> = 2.78 × 10<sup>-7</sup> m s<sup>-1</sup></td><td>564</td><td>35</td><td>-246</td><td>332</td><td>224</td></tr><tr><td>Second-year gravitational K<sub>s</sub> = 2.78 × 10<sup>-7</sup> m s<sup>-1</sup></td><td>424</td><td>12</td><td>-82</td><td>310</td><td>107</td></tr><tr><td>Third-year gravitational K<sub>s</sub> = 2.78 × 10<sup>-7</sup> m s<sup>-1</sup></td><td>378</td><td>4</td><td>-29</td><td>306</td><td>64</td></tr></tbody></table>

simulation was equal to the soil water storage at the end of the simulation (i.e., after 1 yr). This ensured that model results did not depend on the initial moisture conditions and that various scenarios could be compared with each other (Dooge et al. 1994). The same vegetation characteristics and the same lower boundary condition (sinusoidal annual cycle for temperature, and either a constant value of the matric potential or gravitational drainage for the moisture) were retained for all 78 soil profiles.

5) Then, the mean water budget was calculated at the annual and monthly timescales by averaging the results of the 78 profiles, each profile representing 1/78 of the whole area. The variability of those components of the water budget was assessed. A comparison with a single simulation using lumped parameters derived using various methods was conducted. Simulations using lumped parameter values associated with the average texture, as given by the Clapp and Hornberger (1978) classification, were also performed. The differences in the water budget obtained with the various aggregation scenarios and the average of the 78 profiles outputs could therefore be quantified.

Milly and Eagleson (1987), followed by Kim et al. (1997), proposed an analytical framework for the derivation of effective parameters, based on reasonable assumptions about the probability density functions of the soil parameters. Such a generality was not sought in the present paper, which aimed only to test the relevance of possible parameter aggregations readily available from measured data.

## 3. Results

### a. Analysis of the annual water balance

The assumptions retained in the modeling approach constrained greatly the modeled annual water balance. Due to equilibrium, annual change in water storage was 0. Rainfall was the same for all the soil profiles, because it was chosen to focus the analysis on soil surface properties influence. Calculated runoff was 0 in all the cases,

because the soil saturated hydraulic conductivity was very high (more than 120 mm day<sup>-1</sup>) as compared to rainfall intensity (less than 20 mm day<sup>-1</sup>). Therefore, the variability of the water balance was reduced to a balance between total evaporation and deep drainage. Deep drainage was proportional to the hydraulic conductivity at the bottom of the soil column. Soil characteristics of the lower layer were the same for all the soil profiles. Therefore, the lower boundary condition choice and the value of the saturated hydraulic conductivity of this layer mainly governed deep drainage. Table 2 shows that if a constant matric potential was assumed at the bottom of the soil profile, capillary rises were generated (upward flux), whereas the flux was directed downward if gravitational flux was assumed. In the latter case, drainage at the end of the first year decreased from 146 to 35 mm when the bottom saturated hydraulic conductivity was divided by 10. When equilibrium was reached, variability of drainage was low (Tables 4–6). As a consequence, a change in surface soil hydraulic properties mainly affected the partition of total evaporation between bare soil evaporation and transpiration by the plants.

#### 1) INFLUENCE OF LOWER BOUNDARY CONDITION CHOICE

Table 2 provides the terms of the mean annual water budget (average over the 78 profiles) for the three choices of lower boundary condition [case 1, constant matric potential (with $K_s = 2.78 \times 10^{-6} \text{ m s}^{-1}$); case 2, gravitational flow with $K_s = 2.78 \times 10^{-6} \text{ m s}^{-1}$; case 3, gravitational flow with $K_s = 2.78 \times 10^{-7} \text{ m s}^{-1}$] and the two/three years of the simulation.

Table 2 shows that quite different results were obtained according to lower boundary condition choice for the water flow equation. In the case of constant matric potential, large capillary rises were simulated (almost 200 mm yr<sup>-1</sup>). This figure was very high and one might suspect that a lower boundary condition defined with a value of the matric potential evolving through the season would be more realistic. Unfortunately, no data were available to define such an annual

<footnote>* Equilibrium was not reached at the end of the second year. However, runs on a third year were not possible because of a divergence of the model due to very small water content within the soil profile.</footnote>
course. When the gravitational flow was considered, the choice was not very satisfactory because only downward fluxes were allowed and capillary rises were systematically excluded. Furthermore, the value of the saturated hydraulic conductivity, $K_s = 2.78 \times 10^{-6} \text{ m s}^{-1}$ (chosen consistently with earlier published work on this area), at the bottom led to large values of drainage, especially for the first year, and then little water remained in the soil profiles for the second year. At the end of the second year equilibrium was not reached. However, it was not possible to run the model for a third year, because of model divergence associated with a very low water content in the soil profiles. Therefore, another trial was done by dividing the value of the saturated hydraulic conductivity of the lower layer by 10 in order to avoid too much drainage.

Boundary condition choice was crucial for such a long-term study and results in terms of water balance were very different according to this choice. When matric potential was constant at the bottom of the soil profile, total evaporation was larger than rainfall, due to almost 200 mm of capillary rise. With a gravitational flow, total evaporation was almost equal to rainfall, because deep percolation was very small once equilibrium was reached. It was also interesting to see that lower boundary condition choice mainly affected plant transpiration, which was greatly reduced in the gravitational flow case, whereas bare soil evaporation was almost the same. The change in mean water storage implied that less water was available for transpiration in deeper layers, whereas bare soil evaporation was more linked to surface soil moisture, and therefore to rainfall time evolution (see also discussion in section 3b).

Finally, it seems that for such a study, the ideal choice for the lower boundary condition would be an imposed value of the matric potential evolving with time, provided that such a time evolution could be defined. This is in fact very difficult in practice, because such an evolution depends on rainfall time evolution and redistribution of water within the soil profiles, that is, on model results themselves. Another choice could be to use a very deep soil profile, where the lower boundary condition could be defined by the water table. However, in this case, the determination of soil hydrodynamic properties for deeper layers becomes a problem because these layers are poorly known and a coupling with a 2D groundwater model could be required.

The problem posed by the specification of the soil lower boundary condition is more crucial for models based on the Richards (1931) equation than for reservoir models, not resolving explicitly the diffusion equation within the soil [e.g., the Interface Soil Biosphere Atmosphere (ISBA) model of Noilhan and Planton 1989]. The predictive potential of Richards equation-based models is therefore limited. Boone and Wetzel (1996) and Lee and Abriola (1999) also reported sensitivity of such models, used within general circulation models, to the lower boundary conditions representation and to the

<table><caption>TABLE 3. Values of the parameters for the three 1D cases. Here “aaaa” means arithmetic mean for all the four parameters, and “aagg” means arithmetic mean for $\lambda$ and $\theta_s$, and geometric mean for $h_{oc}$ and $K_s$, and corresponds to the median value of the four parameters. “Mean CH” represents the Clapp and Hornberger derived parameters for the average soil texture over the whole area (class 3 for both layers).</caption><thead><tr><th></th><th>Layer</th><th>λ</th><th>θ<sub>s</sub><br>(cm<sup>3</sup>)</th><th>-h<sub>we</sub><br>(cm)</th><th>K<sub>s</sub><br>(cm s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="2">aaaa</td><td>1</td><td>0.163</td><td>0.504</td><td>11.8</td><td>0.00151</td></tr><tr><td>2</td><td>0.162</td><td>0.486</td><td>11.7</td><td>0.00111</td></tr><tr><td rowspan="2">aagg</td><td>1</td><td>0.163</td><td>0.504</td><td>7.86</td><td>0.00105</td></tr><tr><td>2</td><td>0.162</td><td>0.486</td><td>9.97</td><td>0.00066</td></tr><tr><td rowspan="2">Mean CH</td><td>1</td><td>0.204</td><td>0.435</td><td>21.8</td><td>0.00341</td></tr><tr><td>2</td><td>0.204</td><td>0.435</td><td>21.8</td><td>0.00341</td></tr></tbody></table>

soil discretization. Nevertheless, results obtained using the SiSPAT SVAT model are worth discussing because some common features can be extracted, independent of the lower boundary condition choice.

## 2) EVALUATION OF SEVERAL AGGREGATION RULES FOR SURFACE SOIL HYDRAULIC PROPERTIES

The average water budget was compared with one 1D run conducted with aggregated parameters. Two choices were tested. Here “aagg” means arithmetic mean for $\lambda$ and $\theta_s$ and geometric mean for $h_{oc}$ and $K_s$, and corresponds to the median of the four parameters, according to their fitted probability density function (see Table 1). The designation “aaaa” means arithmetic mean for all the four parameters. Water budget values derived with parameters calculated from the Clapp and Hornberger (1978) classification for the average soil texture of the two surface layers are also given. Corresponding parameter values are summarized in Table 3. Tables 4–6 provide the statistical analysis of water budget annual components at the end of the second or third year.

For the first lower boundary condition choice (imposed constant matric potential value), the variability in the water balance components was less than 10% and was maximum for bare soil evaporation and deep drainage. Transpiration was not affected very much because, due to capillary rise, vegetation was always well supplied with water. On the other hand, in the second and third cases (gravitational flow), a large amount of water left the soil profiles through percolation and the average water content could become lower than the wilting point. Vegetation water stress appeared for some soil profiles and transpiration variability was the highest (31%–37%) whereas bare soil evaporation showed a lower variability than in the first case (7%). Contrasts in variability between cases 1 and 2 or 3 resulted from the steady state being reached. At the end of the first year (when steady state was not reached), the water balance components’ variability was similar for the three lower boundary condition choices, except for deep percolation (not shown). Note also that although tran-
<table><caption>TABLE 4. Case 1: constant matric potential at the bottom of the soil profile for the second year; statistics of the annual water balance calculated by using measured surface hydraulic properties in the 78 soil profiles. Deep drainage was positive for percolation and negative for capillary rises. Mean values obtained through one 1D run using either the aagg or aaaa averaging or the soil surface hydraulic properties derived from the Clapp and Hornberger classification for the average soil texture are also given.</caption><thead><tr><th></th><th>Total evaporation</th><th>Deep drainage</th><th>Change in water storage</th><th>Bare soil evaporation</th><th>Transpiration</th></tr></thead><tbody><tr><td>Mean</td><td>549</td><td>-192</td><td>-4</td><td>319</td><td>223</td></tr><tr><td>Median</td><td>553</td><td>-195</td><td>-4</td><td>322</td><td>224</td></tr><tr><td>Min</td><td>497</td><td>-139</td><td>-8</td><td>244</td><td>204</td></tr><tr><td>Max</td><td>571</td><td>-216</td><td>-1</td><td>359</td><td>246</td></tr><tr><td>Std dev</td><td>17</td><td>17</td><td>1</td><td>27</td><td>11</td></tr><tr><td>Coefficient of variation (%)</td><td>3</td><td>9</td><td>-36</td><td>9</td><td>5</td></tr><tr><td>1D run aagg</td><td>550</td><td>-193</td><td>-4</td><td>319</td><td>224</td></tr><tr><td>1D run aaaa</td><td>562</td><td>-204</td><td>-5</td><td>335</td><td>220</td></tr><tr><td>1D run mean CH</td><td>569</td><td>-216</td><td>-1</td><td>358</td><td>204</td></tr></tbody></table>

spiration and bare soil evaporation variability might be large, total evaporation variability was always very small, due to compensation effects.

When looking at the results in terms of aggregated parameters, median values for the four soil parameters (aagg) led to the closest agreement between the 1D run with the aggregated parameters and the 78 soil profiles simulations average in the three test cases. The arithmetic mean led to a larger bias in terms of partition between bare soil evaporation and transpiration (average scale parameter $h_{bc}$ and saturated hydraulic conductivity $K_s$ were larger than median values used in the previous case because the probability density functions were log-normal). In general, transpiration was underestimated with the aggregated parameters and bare soil evaporation overestimated. Due to compensation effects, the bias on total evapotranspiration was, however, small (less than 5%, except for run CH in case 2 where it reached 15%). The bias on evaporation and transpiration was small using median values of the parameters (less than 5%), larger using the arithmetic mean (between 5% and 30%), and very large using mean parameters derived from the Clapp and Hornberger classification, especially in the gravitational flow case (case 3) where it reached 55%. In this case, predicted transpiration was half the true value.

Therefore, in the case study where no runoff was

simulated, median values used as aggregated soil surface parameters led to a satisfactory simulation of the water balance with one “equivalent” 1D simulation. This result was consistent with studies reported by Milly and Eagleson (1987), Kim and Stricker (1996), Kim et al. (1997), Braud (1998), and Boulet et al. (1999). They showed that runoff occurrence was triggering the effect of the spatial variability of soil properties on surface fluxes. Estimation of the median value for the parameters remains, however, an unsolved problem given the large variability of some of them (see Fig. 4) requiring a large sample for achieving a robust estimation.

### b. Analysis of the monthly water balance

Results are summarized in Fig. 4 for the second year of case 1 (imposed matric potential at the bottom), Fig. 5 for the second year of case 2 (gravitational flow with saturated hydraulic conductivity of $2.78 \times 10^{-6} \text{ m s}^{-1}$), and Fig. 6 for the third year of case 3 (gravitational flow with saturated hydraulic conductivity of $2.78 \times 10^{-7} \text{ m s}^{-1}$). The mean monthly water balance of the 78 soil profiles ± one standard deviation are shown, as well as the monthly water balance obtained with the three choices of aggregated parameters. Rainfall time evolution is also shown. Bare soil evaporation appeared mainly related to rainfall, whereas transpiration time

<table><caption>TABLE 5. Case 2: gravitational flow at the bottom of the soil profile with $K_s = 2.8 \times 10^{-6} \text{ m s}^{-1}$ for the second year statistics of the annual water balance calculated by using measured surface hydraulic properties in the 78 soil profiles. Deep drainage was positive for percolation and negative for capillary rises. Mean values obtained through one 1D run using either the aagg or aaaa averaging or the soil surface hydraulic properties derived from the Clapp and Hornberger classification are also given.</caption><thead><tr><th></th><th>Total evaporation</th><th>Deep drainage</th><th>Change in water storage</th><th>Bare soil evaporation</th><th>Transpiration</th></tr></thead><tbody><tr><td>Mean</td><td>410</td><td>6</td><td>-63</td><td>303</td><td>97</td></tr><tr><td>Median</td><td>410</td><td>6</td><td>-63</td><td>305</td><td>94</td></tr><tr><td>Min</td><td>393</td><td>4</td><td>-80</td><td>240</td><td>46</td></tr><tr><td>Max</td><td>423</td><td>9</td><td>-46</td><td>339</td><td>175</td></tr><tr><td>Std dev</td><td>8</td><td>1</td><td>8</td><td>22</td><td>30</td></tr><tr><td>Coefficient of variation (%)</td><td>2</td><td>14</td><td>13</td><td>7</td><td>31</td></tr><tr><td>1D run aagg</td><td>403</td><td>6</td><td>-56</td><td>303</td><td>92</td></tr><tr><td>1D run aaaa</td><td>402</td><td>6</td><td>-55</td><td>323</td><td>72</td></tr><tr><td>1D run mean CH</td><td>352</td><td>6</td><td>-59</td><td>302</td><td>43</td></tr></tbody></table>
<table><caption>TABLE 6. Case 3: gravitational flow at the bottom of the soil profile with $K_s = 2.78 \times 10^{-7} \text{ m s}^{-1}$ for the third year. Statistics of the annual water balance calculated by using measured surface hydraulic properties in the 78 soil profiles. Deep drainage was positive for percolation and negative for capillary rises. Mean values obtained through one 1D run using either the aagg or aaaa averaging or the soil surface hydraulic properties derived from the Clapp and Hornberger classification for the average soil texture are also given.</caption><thead><tr><th></th><th>Total evaporation</th><th>Deep drainage</th><th>Change in water storage</th><th>Bare soil evaporation</th><th>Transpiration</th></tr></thead><tbody><tr><td>Mean</td><td>378</td><td>4</td><td>-29</td><td>306</td><td>64</td></tr><tr><td>Median</td><td>378</td><td>4</td><td>-29</td><td>309</td><td>61</td></tr><tr><td>Min</td><td>374</td><td>4</td><td>-39</td><td>245</td><td>29</td></tr><tr><td>Max</td><td>388</td><td>6</td><td>-25</td><td>338</td><td>136</td></tr><tr><td>Std dev</td><td>3</td><td>0.4</td><td>3</td><td>21</td><td>24</td></tr><tr><td>Coefficient of variation (%)</td><td>1</td><td>8</td><td>10</td><td>7</td><td>37</td></tr><tr><td>1D run aagg</td><td>378</td><td>4</td><td>-29</td><td>307</td><td>63</td></tr><tr><td>1D run aaaa</td><td>376</td><td>4</td><td>-28</td><td>322</td><td>47</td></tr><tr><td>Mean Clapp and Hornberger</td><td>375</td><td>4</td><td>-25</td><td>339</td><td>28</td></tr></tbody></table>

evolution was more related to LAI evolution (significantly different from 0 between Apr and Sep). Tables 7–9 provide the monthly coefficients of variation for total evaporation, bare soil evaporation, and transpiration for the last simulation year and the three lower boundary conditions choices.

Bare soil evaporation variability was maximum in May (32%), a month with no rainfall following a month with a large rainfall. This was especially true for case

1 (constant matric potential at the bottom of the soil profile), where transpiration variability was very small, because the vegetation was well supplied with water due to the capillary rise. For cases 2 and 3, transpiration variability was maximum (up to 66%) during the dry period because of plant water stress. Note that in case 3, where the lack of water was larger, the coefficients of variation were also the highest. For total evaporation, the monthly coefficient of variation did not exceed 30%,

<figure><img src="image_5.png" /><figcaption>FIG. 4. Monthly evolution of the mean value (*) ± 1 std dev (error bars) of the 78 soil profiles for total evaporation, bare soil evaporation, transpiration, and deep drainage (positive for percolation and negative for capillary rises). Monthly values of rainfall are also shown together with the monthly evolution of the components of the water balance for the case aagg (open triangle) or aaaa (full square) averaging and the case where the soil hydraulic properties of the surface were derived from the Clapp and Hornberger classification (open diamond). Case 1: imposed matric potential at the bottom of the soil profiles (second year).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Monthly evolution of the mean value (*) ± 1 std dev (error bars) of the 78 soil profiles for total evapotranspiration, bare soil evaporation, transpiration, and deep drainage (positive for percolation and negative for capillary rises). Monthly values of rainfall are also shown together with the monthly evolution of the components of the water balance for the case aagg (open triangle) or aaaa (full square) averaging and the case where the soil hydraulic properties of the surface were derived from the Clapp and Hornberger classification (open diamond). Case 2: gravitational flow with saturated hydraulic conductivity of $2.78 \times 10^{-6} \text{ m s}^{-1}$ (second year).</figcaption></figure>

showing once again that compensation effects occurred between bare soil evaporation and transpiration. In the gravitational flow case, deep drainage variability was very small and was much larger in the imposed matric potential case (not shown).

When looking at the results obtained with the aggregated parameters, aagg and aaaa choices almost always fell in the ± one standard deviation interval, whereas the Clapp and Hornberger results were often outside this interval, especially for transpiration in the dry period. At the monthly scale, aggregated runs performed using parameters' median values led to the best agreement obtained with the 78 soil profiles average.

## 4. Discussion and conclusions

A methodology to derive soil hydraulic properties on a large area was presented. It was shown that soil water retention and hydraulic conductivity curves at a large number of sample points could be derived using simple measurements and the variability could be characterized at the regional scale. Since the early EFEDA campaign in 1994, the methodology has been improved and is fully described in BRAU.

The numerical study conducted using this dataset focused on the influence of variability of surface hydraulic properties on the annual and monthly water budget. Rainfall and vegetation cover were the same for all the soil profiles. Results were discussed at equilibrium when the annual change in water storage was 0 and no runoff was generated, due to high values of the hydraulic conductivity. Consequently, total evaporation and deep drainage variability was less than 10%, regardless of the lower boundary condition, whereas the choice of the lower boundary condition resulted in quite different scenarios for bare soil evaporation and transpiration. When the combination of lower boundary condition and soil hydraulic properties generated water stress for the vegetation on some soil columns, transpiration variability over the whole area was large. When the lower boundary condition ensured that no water stress occurred for the vegetation, transpiration variability over the whole area was small, whereas that of bare soil evaporation variability was higher. At the monthly timescale, differences observed between the various cases were enhanced, leading to a very large variability of transpiration in case of water stress. A limitation of the study is the lack of interaction between surface condition, the atmo-
<figure><img src="image_7.png" /><figcaption>FIG. 6. Monthly evolution of the mean value (*) ± 1 std dev (error bars) of the 78 soil profiles for total evapotranspiration, bare soil evaporation, transpiration, and deep drainage (positive for percolation and negative for capillary rises). Monthly values of rainfall are also shown together with the monthly evolution of the components of the water balance for the case aagg (open triangle) or aaaa (full square) averaging and the case where the soil hydraulic properties of the surface were derived from the Clapp and Hornberger classification (open diamond). Case 3: gravitational flow with saturated hydraulic conductivity of 2.78 × 10⁻⁷ m s⁻¹ (third year).</figcaption></figure>

sphere, and vegetation growth. The use of a vegetation growth model coupled with an atmospheric model would be the next step to investigate further interaction between soil surface properties variability and surface fluxes.

The study has shown that using median parameter values in a 1D run leads to a good agreement with the 78 profiles associated with the measured surface properties average. The use of arithmetic means or worse, of the Clapp and Hornberger classification, leads to a serious bias, especially in case of water stress for the vegetation. Peck et al. (1977) found that for a forest cover an averaging procedure of soil parameters, based on the scaling theory of Miller and Miller (1956), provides a good agreement between average fluxes and fluxes calculated for the equivalent medium. However,

their period of study was much shorter than in our case (only a few months) and the effect of initial conditions on the results was not considered.

When examining the consequences in terms of GCM modeling, it must be stressed that the use of texture-derived parameters (such as the Clapp and Hornberger classification) can lead to significant bias in simulated partition between soil evaporation and transpiration, although total evaporation might be correctly simulated. These results were obtained using an assumption of homogeneous rainfall and vegetation characteristics and no runoff was simulated. The bias reported in this paper was therefore certainly a lower bound. An intercomparison of SVAT models used by climate modelers and hydrologists showed that the runoff term was responsible for the largest differences between models and that

<table><caption>TABLE 7. Case 1: constant matric potential at the bottom of the soil profile for the second year. Monthly coefficients of variation (%) of total evaporation (EVT), bare soil evaporation (BSE), and plant transpiration (TR).</caption><thead><tr><th></th><th>Jan</th><th>Feb</th><th>Mar</th><th>Apr</th><th>May</th><th>Jun</th><th>Jul</th><th>Aug</th><th>Sep</th><th>Oct</th><th>Nov</th><th>Dec</th></tr></thead><tbody><tr><th>EVT</th><td>12.7</td><td>17.3</td><td>7.7</td><td>10.2</td><td>17.1</td><td>3.2</td><td>3.3</td><td>6.6</td><td>4.1</td><td>1.0</td><td>9.9</td><td>8.2</td></tr><tr><th>BSE</th><td>12.9</td><td>17.8</td><td>8.1</td><td>10.8</td><td>32.6</td><td>8.9</td><td>7.5</td><td>9.9</td><td>3.4</td><td>2.8</td><td>12.2</td><td>8.6</td></tr><tr><th>TR</th><td>5.8</td><td>5.6</td><td>1.3</td><td>2.2</td><td>0.9</td><td>1.6</td><td>5.2</td><td>8.9</td><td>7.4</td><td>3.7</td><td>4.3</td><td>1.3</td></tr></tbody></table>
<table><caption>TABLE 8. Case 2: gravitational flow at the bottom of the soil profile with $K_s = 2.78 \times 10^{-6} \text{ m s}^{-1}$ for the second year. Monthly coefficients of variation (%) of total evaporation (EVT), bare soil evaporation (BSE), and plant transpiration (TR).</caption><thead><tr><th></th><th>Jan</th><th>Feb</th><th>Mar</th><th>Apr</th><th>May</th><th>Jun</th><th>Jul</th><th>Aug</th><th>Sep</th><th>Oct</th><th>Nov</th><th>Dec</th></tr></thead><tbody><tr><th>EVT</th><td>12.8</td><td>17.6</td><td>5.9</td><td>10.7</td><td>11.1</td><td>12.5</td><td>32.8</td><td>27.8</td><td>6.9</td><td>2.3</td><td>6.5</td><td>3.7</td></tr><tr><th>BSE</th><td>13.0</td><td>18.2</td><td>6.2</td><td>11.3</td><td>22.2</td><td>4.6</td><td>4.3</td><td>6.4</td><td>2.4</td><td>1.8</td><td>11.9</td><td>5.8</td></tr><tr><th>TR</th><td>5.8</td><td>5.5</td><td>1.1</td><td>2.6</td><td>1.7</td><td>24.5</td><td>54.0</td><td>52.3</td><td>34.4</td><td>25.0</td><td>44.1</td><td>40.6</td></tr></tbody></table>

when rainfall variability was taken into account, simulated runoff was considerably modified (Dooge et al. 1994).

Given the large soil properties spatial variability, and therefore the large sample needed to obtain a representative value, the experimental effort needed to get the median value is the major obstacle to the derivation of such values for large areas. The methodology presented in the first part of this paper represents the smallest effort achievable at the present time, but is still difficult to apply in a large number of regions. Future research should probably focus on the derivation of a representative set of soil parameters at the scale at which the modeling is conducted, provided that these units represent homogeneous areas. As an illustration, Soria et al. (2002) showed that the fluxes calculated for an equivalent medium of two soils exhibiting the same structural parameters and different textural parameters could not match the reference flux deduced from a full 2D model of water transport, whereas the average flux was well represented by the weighted average of the individual vertical fluxes. These results plead for parameterization of subgrid processes as weighed averages of individual components (combination of soil type and vegetation cover, for instance) referenced as the mosaic approach (Koster and Suarez 1992). This seems more promising than increasing model complexity, assuming homogeneous surfaces. Information on soil hydraulic properties will probably still remain necessary. A research line that requires further development (viz., in order to restrain the number of independent parameters) is certainly the inversion of SVAT models in order to match fluxes estimated at the scale of interest, which could be provided by remote sensing. The availability of evaporation fluxes at the scale of the model could be of great help in reaching this goal and for GCMs' validation, but such remote sensing estimations are still not reliable enough. Assimilation of remote sensing within SVAT or hydrological models data could also be a promising way of achieving this goal, provided that the number of parameters

to be estimated remains sufficiently small to obtain robust estimations.

*Acknowledgments.* This study was funded by the European Union in the context of the ECHIVAL Experiment in a Threatened Desertification Area (EFEDA) project (EPOCH-CT 90-0030 and EV5V-CT93-0272). F. N. Dalton, R. Garcia, and P. J. Ross are thanked for their contribution to the field work. Three anonymous reviewers helped to improve the quality of the paper.

## REFERENCES

Avissar, R., 1992: Conceptual aspects of a statistical-dynamical approach to represent landscape subgrid-scale heterogeneities in atmospheric models. *J. Geophys. Res.*, **97** (D3), 2729–2742.

—, 1995: Scaling of the land-atmosphere interactions: An atmospheric modeling perspective. *Hydrol. Processes*, **9**, 679–695.

Bell, K. R., B. J. Blanchard, T. J. Schmugge, and M. W. Witczake, 1980: Analysis of surface moisture variations within large field sites. *Water Resour. Res.*, **16**, 796–810.

Bolle, H. J., and Coauthors, 1993: EFEDA: European Field Experiment in a Desertification Threatened Area. *Ann. Geophys.*, **11**, 173–189.

Boone, A., and P. J. Wetzel, 1996: Issues related to low resolution modeling of soil moisture: Experience with the PLACE model. *Global Planet. Change*, **13**, 161–181.

Boulet, G., J. D. Kalma, I. Braud, and M. Vauclin, 1999: An assessment of effective parameterization of soil physical and land surface properties in regional-scale water balance studies. *J. Hydrol.*, **217**, 225–238.

Braud, I., 1998: Spatial variability of surface fluxes versus spatial variability of surface properties: Application to a fallow savannah of the Hapex-Sahel experiment using the SiSPAT SVAT model. *Agric. For. Meteor.*, **89**, 15–44.

—, 2000: SiSPAT version 3.0, user's manual. LTHE, 83 pp. [Available from LTHE, BP 53, 38041 Grenoble Cedex 9, France, and online at http://www.lthe.hmg.inpg.fr.]

—, R. Angulo-Jaramillo, R. Haverkamp, J. P. Laurent, J. Noilhan, and J. P. Vandervaere, 1995a: Modélisation locale 1-D des transferts de masse et d'énergie d'une vigne à Tomelloso (EFEDA) incluant une croûte profonde de calcite. *Proc. Atelier de Modélisation de l'Atmosphère*, Toulouse, France, 173–182.

—, A. C. Dantas-Antonino, and M. Vauclin, 1995b: A stochastic approach to studying the influence of the spatial variability of

<table><caption>TABLE 9. Case 3: gravitational flow at the bottom of the soil profile with $K_s = 2.78 \times 10^{-7} \text{ m s}^{-1}$, third year. Monthly coefficients of variation (%) of total evaporation (EVT), bare soil evaporation (BSE), and plant transpiration (TR).</caption><thead><tr><th></th><th>Jan</th><th>Feb</th><th>Mar</th><th>Apr</th><th>May</th><th>Jun</th><th>Jul</th><th>Aug</th><th>Sep</th><th>Oct</th><th>Nov</th><th>Dec</th></tr></thead><tbody><tr><th>EVT</th><td>12.0</td><td>16.7</td><td>5.3</td><td>8.4</td><td>13.4</td><td>21.4</td><td>29.1</td><td>15.8</td><td>4.3</td><td>1.8</td><td>9.9</td><td>7.1</td></tr><tr><th>BSE</th><td>12.2</td><td>17.5</td><td>5.7</td><td>9.0</td><td>20.8</td><td>3.5</td><td>3.6</td><td>6.0</td><td>2.2</td><td>1.9</td><td>12.8</td><td>7.9</td></tr><tr><th>TR</th><td>6.2</td><td>6.1</td><td>1.7</td><td>3.1</td><td>16.1</td><td>50.9</td><td>66.5</td><td>38.5</td><td>26.3</td><td>20.6</td><td>32.7</td><td>14.1</td></tr></tbody></table>
soil hydraulic properties on surface fluxes, temperature and humidity. *J. Hydrol.*, **165**, 283-310.

—, —, —, J. L. Thony, and P. Ruelle, 1995c: A Simple Soil-Plant-Atmosphere Transfer model (SiSPAT): Development and field verification. *J. Hydrol.*, **166**, 231-250.

Brooks, R. H., and A. T. Corey, 1964: Hydraulic properties of porous media. Hydrology Paper 3, Colorado State University, Fort Collins, CO, 27 pp.

Clapp, R. B., and G. M. Hornberger, 1978: Empirical equations for some soil hydraulic properties. *Water Resour. Res.*, **14**, 601–604.

Deardorff, J. W., 1978: Efficient prediction of ground surface temperature and moisture with inclusion of a layer of vegetation. *J. Geophys. Res.*, **83**, 1889–1903.

Dooge, J. C. I., P. R. Rowntree, M. Vauclin, E. Todini, K. Dümenil, J. C. André, and H. Stricker, 1994: Spatial Variability of Land Surface Processes (SLAPS II). *Proc. Symp. on Global Change: Climate Change and Climate Change Impacts*, Copenhagen, Denmark, EEC.

Entekhabi, D., and P. S. Eagleson, 1989: Land surface hydrology parameterization for atmospheric general circulation models including subgrid scale spatial variability. *J. Climate*, **2**, 816–831.

Famiglietti, J. S., and E. F. Wood, 1991: Evapotranspiration and runoff from large land areas: Land surface hydrology for atmospheric general circulation models. *Land-Surface-Atmosphere Interactions for Climate Modelling: Observations, Model and Analysis*, E. F. Wood, Ed., Kluwer Academic, 179–204.

—, and —, 1992: Effects of spatial variability and scale on areally averaged evapotranspiration. *Water Resour. Res.*, **31**, 699–712.

Federer, C. A., 1979: A soil-plant-atmosphere model for transpiration and availability of soil water. *Water Resour. Res.*, **15**, 555–562.

Gee, G. W., and J. W. Bauder, 1986: Particle-size analysis. *Methods of Soil Analysis: Physical and Mineralogical Methods*, 2d ed., A. Klute, Ed., American Society of Agronomy, 383–411.

Green, W. H., and G. A. Ampt, 1911: Studies on soil physics. *J. Agric. Sci.*, **4**, 1–24.

Haverkamp, R., P. J. Ross, K. R. J. Smetten, and J. Y. Parlange, 1994: Three-dimensional analysis of infiltration from the disc infiltrometer. II. Physically based infiltration equation. *Water Resour. Res.*, **30**, 2931–2935.

—, and Coauthors, 1996: Hydrological and thermal behaviour of the vadose zone in the area of Barrax and Tomelloso (Spain): Experimental study, analysis and modelling. Final integrated Rep. EFEDA II Spain, Project UE No. EV5C-CT 92 00 90, 144–185.

—, F. Bouraoui, R. Angulo-Jaramillo, C. Zammit, and J. W. Delleur, 1998a: Soil properties and moisture movement in the unsaturated zone. *The Handbook of Groundwater Engineering*, J. W. Delleur, Ed., CRC Press, chap. V, 1–50.

—, J.-Y. Parlange, R. Cuenca, P. J. Ross, and T. S. Steenhuis, 1998b: Scaling of the Richards equation and its application to watershed modeling. *Scale Dependence and Scale Invariance in Hydrology*, G. Sposito, Ed., Cambridge University Press, 190–223.

—, C. Zammit, F. Bouraoui, K. Rajkai, and J. L. Arrué, 1998c: Survey of soil field data and description of particle size, soil water retention and hydraulic conductivity functions. GRIZZLY: Grenoble Soil Catalogue of Soils, LTHE, 117 pp. [Available from LTHE, BP 53, 38041 Grenoble Cédex 9, France.]

Hillel, D., 1980: *Application of Soil Physics*. Academic Press, 385 pp.

Kabat, P., R. W. A. Hutjes, and R. A. Feddes, 1997: The scaling characteristics of soil parameters. *J. Hydrol.*, **190**, 363–396.

Kim, C. P., and J. N. M. Stricker, 1996: Influence of spatially variable soil hydraulic properties and rainfall intensity on the water budget. *Water Resour. Res.*, **32**, 1699–1712.

—, —, and R. A. Feddes, 1997: Impact of soil heterogeneity on the water budget of the unsaturated zone. *Water Resour. Res.*, **53**, 991–999.

Koster, R. D., and M. J. Suarez, 1992: A comparative analysis of two land-surface heterogeneity representations. *J. Climate*, **5**, 1379–1390.

Lean, J., 1992: A guide to the UK Meteorological Office Single Column Model. Hadley Centre for Climate Prediction and Research, Bracknell, United Kingdom.

Lee, D. H., and L. M. Abriola, 1999: Use of the Richards equation in land surface parameterization. *J. Geophys. Res.*, **104** (D22), 27 519–27 526.

Lewan, E., 1996: Evaporation, discharge and nitrogen leaching from a sandy soil in Sweden: Simulations and measurements at different scales in space and time. Department of Soil Science, Reports and Dissertations 27, University of Uppsala, Uppsala, Sweden, 125 pp.

Linder, W., and Coauthors, 1995: Intercomparison of surface schemes using EFEDA flux data. Note du Groupe de Modélisation à Moyenne Echelle 34, Météo-France/CNRM, 105 pp.

Mallants, D., B. P. Mohanty, D. Jacques, and J. Feyen, 1996: Spatial variability of hydraulic properties in a multi-layered soil profile. *Soil Sci.*, **161**, 167–181.

Miller, E. E., and R. D. Miller, 1956: Physical theory for capillary flow phenomena. *J. Appl. Phys.*, **27**, 324–332.

Milly, P. C. D., 1982: Moisture and heat transport in hysteretic inhomogeneous porous media: A matric head-based formulation and a numerical model. *Water Resour. Res.*, **18**, 489–498.

—, and P. S. Eagleson, 1987: Effects of spatial variability on annual average water balance. *Water Resour. Res.*, **23**, 2135–2143.

Nielsen, D. R., J. W. Biggar, and K. T. Erh, 1973: Spatial variability of field-measured soil-water profiles. *Hilgardia*, **42**, 215–260.

Noilhan, J., and S. Planton, 1989: A simple parameterization of land surface processes for meteorological models. *Mon. Wea. Rev.*, **117**, 536–549.

—, and P. Lacarrère, 1995: GCM gridscale evaporation from mesoscale modelling. *J. Climate*, **8**, 206–223.

Peck, A. J., R. J. Luxmoore, and J. L. Stolzy, 1977: Effects of spatial variability of soil hydraulic properties in water budget modelling. *Water Resour. Res.*, **13**, 348–354.

Philip, J. R., and D. A. De Vries, 1957: Moisture movement in porous materials under temperature gradients. *Eos, Trans. Amer. Geophys. Union*, **38**, 222–232.

Raupach, M. R., and J. J. Finnigan, 1995: Scale issues in boundary-layer meteorology: Surface energy balances in heterogeneous terrain. *Hydrol. Processes*, **9**, 589–612.

Richards, L. A., 1931: Capillary conduction of liquids through porous mediums. *J. Phys.*, **1**, 318–333.

Russo, D., and E. Bresler, 1981: Soil hydraulic properties as stochastic processes. I. An analysis of field spatial variability. *Soil Sci. Soc. Amer. J.*, **45**, 682–687.

Santa Olalla Manas, J. F., 1994: EFEDA II Project: Vegetation, soil physics inventory and impact. Technical annual report (August 1 1993–July 31 1994), Albacete, Spain, 242 pp.

Sene, K. J., 1996: Meteorological estimates for the water balance of a sparse vine crop growing in semiarid conditions. *J. Hydrol.*, **179**, 259–280.

Shao, Y., and A. Henderson-Sellers, 1996: Modelling of soil moisture: A project for intercomparison of land surface parameterisation schemes Phase 2(b). *J. Geophys. Res.*, **101** (D3), 7227–7250.

Shuttleworth, W. J., and J. S. Wallace, 1985: Evaporation from sparse crops—An energy combination theory. *Quart. J. Roy. Meteor. Soc.*, **111**, 839–855.

Soria, J. M., P. Reggiani, R. Angulo-Jaramillo, and R. Haverkamp, 2002: Aggregation of soils in presence of spatial changes in texture for representation of large scale watershed hydrological fluxes. *Proc. Int. Association for Hydraulic Research Conf.*, Berkeley, CA, IAHR, 573–578.

Taconet, O., R. Bernard, and D. Vidal-Madjar, 1986: Evapotranspiration over an agricultural region using a surface flux/temperature model based on NOAA-AVHRR data. *J. Climate Appl. Meteor.*, **25**, 284–307.

Topp, G. C., J. L. Davis, and A. P. Annan, 1980: Electromagnetic
determination of soil water content: Measurements in coaxial transmission lines. *Water Resour. Res.*, **16**, 574–582.

Vandervaere, J. P., 1995: Caractérisation hydrodynamique du sol in situ par infiltrométrie à disques. Analyse critique des régimes pseudo-permanents. Méthodes transitoires et cas des sols en-croûtés. Ph.D. thesis, Université Joseph Fourier, Grenoble I, Grenoble, France, 329 pp.

Vauclin, M., D. E. Elrick, J. L. Thony, G. Vachaud, P. Revol, and P. Ruelle, 1994: Hydraulic conductivity measurements of the spatial variability of a loamy soil. *Soil Technol.*, **7**, 181–195.

Zammit, C., 1999: Analyse et évaluation des paramètres hydrodynamiques des sols: Prédiction par un modèle analytique à base physique à partir de données texturales. Ph.D. thesis, Université Joseph Fourier, Grenoble I, Grenoble, France, 200 pp.