## Abstract

Spatially averaged soil moisture dynamics are studied under seasonally fixed conditions. We consider rainfall as a marked Poisson process, uniformly covering a spatial domain consisting of multiple plant types. Each plant type is considered to have different characteristics in terms of evapotranspiration functions, root-zone depth and rainfall interception. Equations for the evolution of joint probability density functions for individual soil moistures associated with different plant types are developed, and the non-uniqueness of the spatially averaged evapotranspiration function as a function of the average soil moisture is demonstrated and quantified in an example.

## 1. Introduction

Large-scale global models attempt to represent the main features of the climate dynamics taking place in the earth–atmosphere system. These models use grid cells of a relative large area (of the size order of 10^{4} km^{2}), where conditions are represented by their average spatial values. An important variable in the above dynamics is the evapotranspiration taking place over the grid cell. Its correct characterization has important implications both in the mass transfer of water between earth and atmosphere as well as for the latent heat-sensible heat partition, which is a key factor in the energy balance over the region.

Evapotranspiration is a function whose structure depends on the characteristics of climate, soil and vegetation. As the area of the region under consideration increases, those controlling factors become more and more heterogeneous in space and larger fluctuations take place in the evapotranspiration inside the region. At a point level, the soil moisture interacting with the vegetation at the site will control the evapotranspiration fluxes. At larger spatial scales, even if the average soil moisture is well known, its interaction with spatially varying vegetation can lead to a wide range of responses. Thus, the same average soil moisture over the region can lead to different total evapotranspiration from the region. Therefore, over some range of length-scales, unique relationships between average soil saturation and evapotranspiration cannot be obtained when the averaging domain includes plants with different evapotranspiration properties.

The main objective of this paper is to provide a probabilistic framework for the characterization of total evapotranspiration over a region with heterogeneous vegetation. We will consider the stochastic aspects of the temporal occurrence of precipitation, but the region will be assumed to be uniformly affected by rainfall in space. While this assumption of spatial uniformity of rainfall limits the length-scale over which averages may be taken, the overall procedure illustrates the underlying non-uniqueness while providing a general framework within which additional spatial variations, including rainfall, can be included.

## 2. Soil moisture balance equation

In a region with heterogeneous vegetation, the relationship between spatially averaged soil moisture and spatially averaged evapotranspiration depends on the relationships between evapotranspiration and saturation for the individual plants within the system. For the analysis presented herein, we assume that the plant root zones are non-overlapping in space, and that the horizontal cross-flow is negligible. Then for each root zone we may write the soil moisture balance equation (e.g. Rodriguez-Iturbe & Porporato 2004),(2.1)where *n* is porosity, *Z* is depth of active soil or root depth, is relative soil moisture content, *t* is time, is rainfall rate, is rate of losses due to canopy interception, is runoff rate and and are rates of evapotranspiration and leakage, respectively. Subscript *i* denotes plant species. Note that we allow the plants to prefer different soil types through the species dependency of the porosity *n*.

Following previous authors (e.g. Rodriguez-Iturbe *et al*. 1999*b*; Guswa *et al*. 2002), we will model the rainfall as a marked Poisson process in time, where the depth of rainfall is exponentially distributed. The canopy interception will be modelled as a shift in the mean rainfall depth. Further, the equations for runoff, leakage and evapotranspiration are all nonlinear and non-decreasing. For each plant type equation (2.1) therefore takes the form of a nonlinear, stochastic ordinary differential equation in the variable . This one-plant equation has been extensively studied in Cox & Isham (1986), Rodriguez-Iturbe *et al*. (1999) and Laio *et al*. (2001*a*). While these equations have been treated in a probabilistic framework individually, we are interested in the relationship between the probability distribution functions they generate, so we can study the influence of multiple plant species coexisting within a spatial averaging area of interest.

### (a) Infiltration from rainfall

We will model rainfall as point processes in time arriving as a Poisson process with rate *λ*, occurring uniformly on all plants. This can be interpreted as an idealization of a climate where the rainfall duration is short compared to the time-scale of evapotranspiration, and where the spatial scale of variability of rainfall is large compared to the area of interest. The rain events are assumed to have varying depths *R*, according to the probability density function(2.2)where is the mean depth of rain events.

Canopy interception will typically reduce observed rainfall for a root zone by a fixed constant dependent on the structure of the plant. While most rain events will thus be reduced by this factor, some rain events might have sufficiently small depths that they would not be observed at all. For the purposes of the treatment herein, we will model this behaviour simply as a shift in the mean rainfall depth, neglecting the fact that some rain events might be below the canopy interception.

In what follows, we will work with the normalized version of equation (2.2):(2.3)where(2.4)is the normalized rainfall, and is the normalized mean rainfall depth, both corrected for canopy interception. Note that since , , and (and thus also ) are constants, the functions will be fully correlated through their relationship to .

### (b) Dimensionless system of equations

With the above description of rainfall, we will proceed by collecting the runoff, leakage and evapotranspiration functions in equation (2.1) into a single term. Dividing through by , we obtain(2.5)where is a known positive monotonic dimensionless function with , and is a dimensionless stochastic forcing term with depth given by . The stationary probability distribution function for a saturation *i* is (eqn. (29) of Laio *et al*. (2001*a*))(2.6)where the constant of integration is determined by requiring the probability density function to be normalized. When the decay function has a simple structure, the integration in equation (2.6) can be performed analytically, as was done for piecewise smooth loss functions in Cox & Isham (1986), Rodriguez-Iturbe *et al*. (1999*a*,*b*) and Laio *et al*. (2001*a*).

## 3. Joint probability distribution function

We will in this section obtain an equation which defines the joint probability density function of the vector of *m* saturations occurring in the root zone of *m* different plant species at a time *t* subject to identical rainfall. The development will be similar to the development in Rodriguez-Iturbe *et al*. (1999*b*), although we will herein be more general as we account for multiple plants and keep the time term throughout the analysis.

The probability density *p* of finding in a given state at time *t* can be related to the probability density at time as(3.1)where *k* is the index indicating a soil saturation satisfying at time *t*, and is therefore independent of . The above expression arises from the weighted sum of the change in saturation when no rain occurs (first term), and the change in saturation due to a point-in-time rainfall (integral term) (Cox & Miller 1965). We have disregarded the atom of probability at the origin, since we see from equation (2.6) that the marginals are finite. If we divide through by , subtract on both sides, divide by and let , we get(3.2)To eliminate the integral, we manipulate the equation using the rainfall depth probability distribution explicitly. If we insert equation (2.3) into (3.2), we obtainFinally, by taking the directional derivative in saturation space of this equation with respect to the vector and adding the equation itself to the result, we arrive at(3.3)Note that taking the directional derivative in the direction leaves the index *k* unchanged.

Equation (3.3) is a linear partial differential equation, which defines the dynamic probability density function for the states .

Equation (3.3) is defined for saturations that are non-negative. However, at steady state, the saturations with non-zero probability density will define a region which is a subset of all possible saturations. This region of non-negative saturations in saturation space will be a function of the evapotranspiration functions and the infiltration from rainfall. To define this region, we will start by looking at the case of two plants with equal root-zone depth and interception (implying ), but different evapotranspiration functions (e.g. figure 1).

First, we partition the saturation space into two regions, dependent on the sign of . These regions are separated by the line of equal evapotranspiration :(3.4)Combining equation (2.5) for each plant, we have that when ,(3.5)We see from this equation that rainfall events change both saturations equally (shifting us along the vector in saturation space). Accordingly, if and only if the saturations are on the line , saturation changes due to evapotranspiration are in the exact opposite direction to the rainfall shifts. However, within each of the regions separated by , the change in the saturation vector with time will have a normal component to the vector , which has uniform sign within the region. This is illustrated in figure 2. A typical path in saturation space is shown in figure 3.

Motivated by this simple discussion of the movement in saturation space, we see that saturations have a tendency to go towards the line . Any saturation on can therefore occur after a (possibly infinite) sequence of rainfall and drying events. Therefore, any drying curve or rainfall shift from this curve will also lead to physical saturation combinations. Finally, revisiting equation (3.5) shows us that all saturations consistent with the equations, which are not artefacts of the initial conditions, must be spanned by drying, and rainfall curves originating from . The resulting domain of physically possible saturation states is illustrated in figure 4 for the evapotranspiration functions in figure 1.

The extension of the region of accessible saturations, *Ω*^{m}, to allow for *m* plants and different follows the same reasoning as above; however, the expressions become slightly more technical. Extend to be the subset of , such that(3.6)where ** s**\

*s*

_{k}is the

*m*−1 dimensional saturation vector obtained from removing the

*k*th dimension of

*, and*

**s***Ω*

^{m−1}is the corresponding

*m*−1 dimensional space of accessible saturations. Here we have modified equation (3.4) to take into account that the saturation shifts due to rainfall are, in general, along the vector . The sets defined by (3.6) will be analogous to those from equation (3.4) in the sense that they partition the saturation space into domains where the net effect of rainfall and evapotranspiration sequences have different directions.

For any point in saturation space, let the evapotranspiration curve be defined by(3.7)and the precipitation line as(3.8)We will denote the set of points on these curves for as and , respectively. Then the accessible saturation states of the *m*-plant system are(3.9)From the definition of the rainfall curves, we see that extends above unity to . This can be interpreted as pooling of water above the top of the soil column (possibly to an infinite height). Alternatively, saturations limited to the interval can be enforced by letting the leakage term of the soil balance equation go to infinity.

At this point, we again emphasize the meaning of the region . This region bounds possible saturations in the sense that for any saturation vector in , no chain of rainfall events can ever result in a saturation outside of . However, initial conditions may contain saturations outside of with positive probability densities. We may infer from the above discussion, and in particular from equation (3.5), that the probability of saturations outside of must be a monotonically decreasing function with time, approaching zero as time goes to infinity.(3.10)

As boundary conditions, equation (3.3) must satisfy the Chapman–Kolmogorov equations(3.11)Further, if we assume that the potential evapotranspiration is greater than the average rainfall, we have that(3.12)

When we only consider a single plant, *m*=1, then equation (3.3) reduces to an extension of the one found in, for example, Rodriguez-Iturbe *et al*. (1999*a*,*b*),(3.13)The stationary form obtained by setting the time derivatives equal to zero can be integrated up to yield the solution reproduced in equation (2.6).

In the case of multiple plants, analytical solutions to equation (3.3) are not, in general, obtainable. For the nonlinear loss functions we are interested in, we will therefore have to use numerical methods to obtain the probability density function *p*.

## 4. Average evapotranspiration function

Because of the complexity of a heterogeneous system, the term ‘average soil moisture’ can have multiple meanings. For some applications one may wish to take into account that the root zones of different plants have different depths, while this may be undesirable in other settings. We will leave the exact definition to be determined by the problem at hand, assuming only that the average saturation can be expressed as a weighted average of the saturations ,(4.1)where the vector contains weights, such that and for all *i*. From this definition we see that fixing an average saturation defines an -dimensional plane in saturation space with normal vector . Let us denote the points on such planes as .

In a similar fashion, we assume that average evapotranspiration can be expressed as a linear combination of ,(4.2)such that and for all *i*. Fixing an average evapotranspiration then defines non-planar surfaces or volumes in or *m* dimensions, respectively. We will denote the set of points covered by these surfaces and volumes as .

We can now obtain the probability density function of average evapotranspiration for a given average saturation,(4.3)The numerator of this equation is simply the integral of the probability density function over all saturations with the desired average saturation and evapotranspiration, while the denominator is the appropriate normalization.

## 5. Application to the Nylsveley savanna

The Nylsveley savanna in South Africa is a warm savanna with relatively homogeneous soil texture, yet with a presence of both woody and herbaceous plants (Rodriguez-Iturbe & Porporato 2004). Both the savanna itself and the physiological characteristics of the plants are well documented elsewhere (e.g. Scholes & Walker 1993; Rodriguez-Iturbe *et al*. 1999*a* and Laio *et al*. 2001*b*). The extensive literature on this savanna, as well as its varied species makes this a good example for study.

While the theory presented above is general in the sense that it allows for time-varying probability distributions, we wish to emphasize herein the non-uniqueness of the average evapotranspiration as a function of average saturation, rather than the time evolution of this relationship. Therefore, we will limit these examples to equilibrium distributions.

This example will serve two main purposes. The first is to illustrate the magnitude of uncertainty in the average relationship. Second, we wish to investigate if the solution to equation (3.3) can be approximated by the product of the marginals, i.e. we will investigate if(5.1)where and solve equation (2.6).

We have chosen to use the plants *Ochna pulchra* and *Eragostris pallens* (Laio *et al*. 2001*b*), which are common tree and grass types, respectively. Simplified forms of the evapotranspiration functions for these plants are shown in figure 5, and all the data used in this example are reproduced in table 1. We see that the evapotranspiration is approximated by a piecewise linear function within the low saturation range typical for the climate in the Nylsveley savanna. Thus, the primary cause of non-uniqueness in the average saturation evapotranspiration relationship will be due to effects around the transition zone between full evapotranspiration and stress-limited evapotranspiration.

We found the steady-state probability density function satisfying equations (3.3) and (3.10)–(3.12) numerically. For the ecosystem in question, extensive numerical simulations have shown that goes to zero at (M. Puma, 2005, personal communication). Therefore, we simplified the solution by solving on a finite domain and replacing boundary condition (3.12) with the assumption that *p* is continuous at , implying that on . The resulting distribution is shown in figure 6. For comparison, the marginal distributions and are plotted in figure 7. These marginal distributions look slightly different than in, for example, Laio *et al*. (2001*b*), since those authors included a different treatment of evapotranspiration for low saturations.

The root-zone depth of both *O. pulchra* and *E. pallens* is 100 cm, thus the average soil saturation will be the average of the soil saturations for the two plants weighted by the proportion of plants. For this illustration we set the proportion of plants equal, so that .

We may now apply equation (4.3) to obtain probability density functions for average evapotranspiration, such as the one shown in figure 8. We see that while the mean of the average evapotranspiration distribution is around 0.4 cm d^{−1}, there is a significant uncertainty associated with this number. In fact, average evapotranspiration values may occur that are lower than the evapotranspiration of either plant at this saturation. The probability density function obtained by using the product of the marginals allows more than twice as wide range of evapotranspiration values, and is clearly not a good approximation in this case.

When we consider the mean and the standard deviations of the average relationship as shown in figure 9, we see that the mean average evapotranspiration deviates up to 20% from the arithmetic average. As expected, the largest standard deviation occurs close to the point where the individual plant evapotranspiration curves have discontinuous derivatives. In fact, we see that at high average saturations, evapotranspiration can be predicted with virtually deterministic certainty. This contrasts with the observations drawn from figure 10, where we again see that the product of the marginals is not a satisfactory approximation to the joint probability. This approximation consistently leads to a gross underestimation of evapotranspiration and overestimation of standard deviation.

## 6. Concluding remarks

A general framework has been presented that allows for the computation of the time-varying probability distributions of average soil saturation and evapotranspiration over a region with heterogeneous vegetation coverage and spatially uniform rainfall. The theory allows for the calculation of the probability density function of average evapotranspiration over the region for a given average soil saturation. This is particularly interesting for studies where remote-sensing measurements provide the average saturation and dynamic models require the characterization of evapotranspiration.

It was observed that the variance of the above distribution may be quite large, implying that there is considerable uncertainty associated with the estimation of an average spatial evapotranspiration. It was also observed that the mean average evapotranspiration for a given saturation may deviate substantially from the arithmetic average calculated from the evapotranspiration values associated with the different plants in the region. The magnitude of the discrepancy will depend on the vegetation and climate characteristics, as well as in the saturation at which the average evapotranspiration is calculated.

The approach makes clear the need to use a joint probability density to characterize the soil saturation. In case the mean average evapotranspiration is calculated using the marginal probability density functions of soil saturation for the different plant types, the result grossly underestimates the average evapotranspiration for a given soil saturation over an extended range of saturation values.

## Acknowledgments

I.R.-I. acknowledges the support of the USA National Science Foundation under grants of the National Centre for Earth Surface Dynamics (EAR-0120914) and Biocomplexity (DEB-0083566).

## Footnotes

- Received January 12, 2005.
- Accepted December 13, 2005.

- © 2006 The Royal Society