Non-uniqueness of evapotranspiration due to spatial heterogeneity of plant species

Jan Martin Nordbotten, Ignacio Rodriguez-Iturbe, Michael A Celia


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 104 km2), 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),Embedded Image(2.1)where n is porosity, Z is depth of active soil or root depth, Embedded Image is relative soil moisture content, t is time, Embedded Image is rainfall rate, Embedded Image is rate of losses due to canopy interception, Embedded Image is runoff rate and Embedded Image and Embedded Image 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. 1999b; 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 Embedded Image. This one-plant equation has been extensively studied in Cox & Isham (1986), Rodriguez-Iturbe et al. (1999) and Laio et al. (2001a). 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 functionEmbedded Image(2.2)where Embedded Image is the mean depth of rain events.

Canopy interception will typically reduce observed rainfall for a root zone by a fixed constant Embedded Image 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 Embedded Image 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):Embedded Image(2.3)whereEmbedded Image(2.4)is the normalized rainfall, and Embedded Image is the normalized mean rainfall depth, both corrected for canopy interception. Note that since Embedded Image, Embedded Image, Embedded Image and Embedded Image (and thus also Embedded Image) are constants, the functions Embedded Image will be fully correlated through their relationship to Embedded Image.

(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 Embedded Image, we obtainEmbedded Image(2.5)where Embedded Image is a known positive monotonic dimensionless function with Embedded Image, and Embedded Image is a dimensionless stochastic forcing term with depth given by Embedded Image. The stationary probability distribution function for a saturation i is (eqn. (29) of Laio et al. (2001a))Embedded Image(2.6)where the constant of integration is determined by requiring the probability density function to be normalized. When the decay function Embedded Image 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. (1999a,b) and Laio et al. (2001a).

3. Joint probability distribution function

We will in this section obtain an equation which defines the joint probability density function Embedded Image of the vector Embedded Image 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. (1999b), 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 Embedded Image in a given state at time t can be related to the probability density at time Embedded Image asEmbedded Image(3.1)where k is the index indicating a soil saturation satisfying Embedded Image at time t, and is therefore independent of Embedded Image. 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 Embedded Image are finite. If we divide through by Embedded Image, subtract Embedded Image on both sides, divide by Embedded Image and let Embedded Image, we getEmbedded Image(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 obtainEmbedded ImageFinally, by taking the directional derivative in saturation space of this equation with respect to the vector Embedded Image and adding the equation itself to the result, we arrive atEmbedded Image(3.3)Note that taking the directional derivative in the direction Embedded Image leaves the index k unchanged.

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

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 Embedded Image), but different evapotranspiration functions (e.g. figure 1).

Figure 1

The dimensionless loss functions for two plants.

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

Figure 2

The solid line from the origin indicates the line Embedded Image, where Embedded Image. This separates saturation space into a region where Embedded Image (right of Embedded Image), and opposite (left of Embedded Image). The straight arrows pointing parallel to Embedded Image indicate the shift in saturation space induced by rainfall, while the curved arrows show movement in saturation space when no rain is falling. The net effect of a rainfall and evapotranspiration event is indicated, and has opposing direction in the two regions.

Figure 3

This figure illustrates a possible path in saturation space for starting at saturation A. We see that during the three rain-events AB, CD and EF, the saturation changes instantaneously along parallel lines, while between rainfalls the rainfall moves on curves with slopes dependent on position in saturation space. Note how the slope is parallel to the rainfall events at the intersection with Embedded Image (point G).

Motivated by this simple discussion of the movement in saturation space, we see that saturations have a tendency to go towards the line Embedded Image. Any saturation on Embedded Image 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 Embedded Image. The resulting domain of physically possible saturation states is illustrated in figure 4 for the evapotranspiration functions in figure 1.

Figure 4

Illustration of possible saturations. Dotted and dashed lines are the evapotranspiration and rainfall paths originating from Embedded Image, respectively. The domain of all possible saturations is the shaded region.

The extension of the region of accessible saturations, Ωm, to allow for m plants and different Embedded Image follows the same reasoning as above; however, the expressions become slightly more technical. Extend Embedded Image to be the subset of Embedded Image, such thatEmbedded Image(3.6)where s\sk is the m−1 dimensional saturation vector obtained from removing the kth dimension of s, and Ω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 Embedded Image. The sets Embedded Image 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 Embedded Image in saturation space, let the evapotranspiration curve Embedded Image be defined byEmbedded Image(3.7)and the precipitation line Embedded Image asEmbedded Image(3.8)We will denote the set of points on these curves for Embedded Image as Embedded Image and Embedded Image, respectively. Then the accessible saturation states Embedded Image of the m-plant system areEmbedded Image(3.9)From the definition of the rainfall curves, we see that Embedded Image extends above unity to Embedded Image. 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 Embedded Image 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 Embedded Image. This region bounds possible saturations in the sense that for any saturation vector in Embedded Image, no chain of rainfall events can ever result in a saturation outside of Embedded Image. However, initial conditions may contain saturations outside of Embedded Image 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 Embedded Image must be a monotonically decreasing function with time, approaching zero as time goes to infinity.Embedded Image(3.10)

As boundary conditions, equation (3.3) must satisfy the Chapman–Kolmogorov equationsEmbedded Image(3.11)Further, if we assume that the potential evapotranspiration is greater than the average rainfall, we have thatEmbedded Image(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. (1999a,b),Embedded Image(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 Embedded Image can be expressed as a weighted average of the saturations Embedded Image,Embedded Image(4.1)where the vector Embedded Image contains weights, such that Embedded Image and Embedded Image for all i. From this definition we see that fixing an average saturation defines an Embedded Image-dimensional plane in saturation space with normal vector Embedded Image. Let us denote the points on such planes as Embedded Image.

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

We can now obtain the probability density function of average evapotranspiration for a given average saturation,Embedded Image(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. 1999a and Laio et al. 2001b). 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 ifEmbedded Image(5.1)where Embedded Image and Embedded Image solve equation (2.6).

We have chosen to use the plants Ochna pulchra and Eragostris pallens (Laio et al. 2001b), 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.

Figure 5

Evapotranspiration functions for Ochna pulchra (plant 1) and Eragostris pallens (plant 2).

View this table:
Table 1

Plant parameters for Ochna pulchra (1) and Eragostris pallens (2). (The frequency of rainfall is 0.167 d−1, with mean depth 1.5 cm.)

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 Embedded Image goes to zero at Embedded Image (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 Embedded Image, implying that Embedded Image on Embedded Image. The resulting distribution Embedded Image is shown in figure 6. For comparison, the marginal distributions Embedded Image and Embedded Image are plotted in figure 7. These marginal distributions look slightly different than in, for example, Laio et al. (2001b), since those authors included a different treatment of evapotranspiration for low saturations.

Figure 6

Joint probability density function of saturation.

Figure 7

Individual probability density functions (marginal distributions) of saturation.

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 Embedded Image.

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 Embedded Image 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.

Figure 8

Probability density function of average evapotranspiration given an average soil moisture of Embedded Image, based on the joint probability distribution (3.3) and the product of the marginal distributions (5.1).

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.

Figure 9

Mean average evapotranspiration obtained using the probability density function shown in figure 6. The arithmetic mean of the plant evapotranspirations given in figure 5 is included as a dashed line for comparison.

Figure 10

Mean average evapotranspiration obtained using the product of the marginal probability density functions shown in figure 7. The arithmetic mean of the plant evapotranspirations given in figure 5 is included as a dashed line for comparison.

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.


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).


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


View Abstract