## Abstract

In semi-arid ecosystems, successful use of the limited water resources is of central importance in determining the evolutionary trends of the vegetation. The competition between different species and individuals for this resource is driven by variations in physiology and metabolic regulation strategies, expressed by such parameters as rooting depth, wilting point or stomatal opening, among others. It is typically not practical to simulate the full evolutionary dynamics of every plant individual in the landscape because of the difficulties introduced by the spatial heterogeneity, as well as the many timescales involved, ranging from hourly up to intergenerational. Instead, the amount of biomass of a given species assimilated per unit area of the landscape may serve as a proxy for its competitiveness and evolutionary success. It is the behaviour of the biomass, which must be described probabilistically due to the stochasticity of the rainfall, which is the subject here. This paper develops a new analytical description of the stationary and transient joint behaviour of plant biomass and soil moisture. Additionally, the effects of climatic fluctuations are considered, including the important case of a bi-seasonal climate regime consisting of alternating wet and dry seasons, which is characteristic of many ecosystems of interest.

## 1. Introduction

Characterizing the vegetation in water-limited ecosystems, with regards to quantity, species composition and stability, is a long-standing problem in ecology, in which a variety of techniques and paradigms have been brought to bear. Some authors have analysed the problem in terms of general organizing principles such as minimizing water stress or maximizing average biomass [1,2], though it is increasingly argued (e.g. [3,4]) that semi-arid or savannah-type ecosystems can be better understood in terms of detailed dynamical behaviour, involving the influences of factors such as water availability, fire, grazing or competitive invasion. As discussed by Accatino *et al.* [5], among the various works which follow this course, there is a dichotomy between models which emphasize disturbances (especially due to fire) and recovery therefrom, and those which emphasize resource competition (especially over water). In particular, Accatino *et al.* [5] note that fire processes are often invoked to explain complex composition structures like tree–grass coexistence (e.g. [6]) when water limitation alone ostensibly cannot; this observation is generalized to the notion that savannahs reflect states of disequilibrium rather than equilibrium.

In this work, by contrast, the water-limited model is re-explored with an emphasis on developing new analytical results for the soil moisture and biomass probabilities induced by a stochastic rainfall process. The primary objective is to obtain and describe the distribution functions generated by equations (2.3) and (2.4), as given in the next section, for three different cases. The first (§3) is the long-time limit, in which the system evolves to an equilibrium (steady state) decoupled from the initial conditions and under the influence of constant parameters. The second case (§4) corresponds to the short-time or transient behaviour, for which the initial conditions play a role, and an approximation scheme is derived which exploits the time-scale separation between the soil moisture and biomass processes. The third case (§5) is an extension of the transient case which allows for time-varying physiological and climate parameters, and culminates with a treatment of seasonal variability, in which the long-term interplay between the wet and dry seasons once again generates an equilibrium, bringing the analysis full circle. These results not only simplify the task of making predictions for landscape behaviour, but also illuminate the sensitivities of the results to the physiological, climatic and soil parameters which govern the system. These sensitivities are essential for understanding questions of evolutionary stability and competition (e.g. with regards to resilience and invasibility), topics which will be explored in an upcoming paper, currently in preparation, using the framework developed here.

## 2. The dynamic equations: form and interpretation

The dynamics of the plant biomass and soil moisture at a point (or within an averaging area) in a landscape are coupled via the transpiration performed by vegetation, which withdraws water from the soil and converts it into new biomass. Meanwhile, soil moisture is replenished by precipitation, and biomass is lost as old plant tissues are destroyed and shed, and stored energy is consumed. The quantitative features of these system dynamics are captured by the following equations^{1} :
*B* (kg m^{−2}) is the biomass per unit area, and *S* (dimensionless) is the relative soil moisture in the root zone (i.e. the fraction of the pore space there occupied by water), described according to a ‘bucket’ model of soil moisture, vertically integrated over the root zone (e.g. [8,9]); the total volume of soil pore space per unit area in the root zone available to store water is the bucket depth *nZ*_{r}, the product of porosity *n* (dimensionless) and root zone depth *Z*_{r} (mm). The various terms appearing on the right-hand side above can be understood as follows. *I* (d^{−1}) is the infiltration rate (i.e. rainfall post interception) normalized by the bucket depth; it is taken to be a marked Poisson process with arrival rate λ (d^{−1}) and exponential depths with mean size 1/*θ* (dimensionless). The function *η*(*S*) gives the soil moisture dependence of evapotranspiration; it is increasing and convex, and goes to unity at the point where *S* is no longer limiting. The function *h*(*b*)(s^{−1}) captures the biomass dependence of evapotranspiration rate. The function *g*(*B*) (kg m^{−2}) gives the scaling of the net biomass gain rate with the existing biomass; it is usually chosen to be approximately linear for small values of *B*, but may be convex at larger values to account for carrying capacity effects [10]. In light of these definitions, the parameters *α* (s^{−1}) and *β* (s^{−1}) may be understood as, respectively, the non-moisture-limited assimilation rate and the biomass loss rate per unit of biomass in the landscape before *g*(*B*) begins to saturate.

This formulation of the dynamics broadly agrees with those given by other authors (e.g. [11–13]), though some simplifications have been introduced, such as the neglect of leakage or spatial diffusion. The exact form of equations (2.1) and (2.2) is chosen *a priori*, though it will prove to be the most general which permits a steady-state solution by the techniques here (§3). In most of what follows, the choices *η*(*S*)=[*S*−*S*_{w}]_{+}/(*k*+[*S*−*S*_{w}]_{+}), *h*(*B*)=*γg*(*B*) and *g*(*B*)=*B* are made, following Zea-Cabrera *et al.* [4] and Nordbotten *et al.* [14], so that *γ* (kg^{−1}s^{−1}) is the maximum rate of transpiration (normalized by bucket depth) per unit of biomass, and *k* (dimensionless) is the half-saturation point, at which transpiration and assimilation occur at half of their respective maxima, and *S*_{w} is the soil moisture level at which evapotranspiration stops. Introducing also *X*≡*S*−*S*_{w}, the dynamic equations are given by

Before proceeding, it is worth addressing a few interpretive points. Biomass losses due to destruction of woody tissue or plant mortality are not incorporated in this model, and so the biomass considered here should be understood as the non-woody biomass of the leaves and fine roots. The linearity with respect to *b* of the assimilation rate implies that the marginal productivity of each additional unit of biomass is constant, so structural, nutrient and light (i.e. energy) limitations are neglected; water is assumed to be the dominant limiting factor. Leakage and run-off are not explicitly resolved at the outset, but a post hoc correction is given later to account for their influence, and this serves to confirm their smallness for parameter values of interest, thereby justifying their neglect in the dynamic equations.

Finally, and perhaps most importantly, the ostensible neglect of evaporation is worth examination. From an atmospheric demand perspective, evaporation is not really absent; intercepted water, which fails to infiltrate the soil, is necessarily removed by evaporation. While the absence of an evaporation (or other non-transpirative loss) term in equation (2.4) is at first glance somewhat disconcerting, since for *B*=0 soil moisture gain becomes unbounded, the system is still well behaved; the biomass cannot completely decay in finite time, and when a low biomass state is reached, the resulting rise in soil moisture causes it to rebound. Moreover, vegetation can suppress evaporation from the soil by limiting radiation and reducing the surface windspeed [15,16], physically justifying the absence of the evaporation term.

However, it is important to note that if the biomass trajectory goes sufficiently low, these evaporation suppression mechanisms may fail, and the dynamics may undergo a paradigm shift from a state where precipitation is mostly balanced by transpiration to one in which it is mostly balanced by soil evaporation. Cueto-Felgueroso *et al.* [11] show that a basin of attraction can emerge around the non-vegetated state *B*=0, separated from the remaining state space by a region of lower probability. In light of this fact, appendix A uses the results of §3 for a brief example on how this scenario might arise from the dynamics given by equations (2.1) and (2.2), while the results of §§4 and 5 should be considered conditional on the rareness of such crossings into such a non-vegetated state.^{2} With these considerations in mind, attention returns to deriving the promised distributions.

## 3. Steady-state analysis

Under the more general dynamics of equations (2.1) and (2.2), provided there is initially some biomass in the system, and if a steady state exists, then either that steady state may be non-vegetated (*B*=0) or vegetated (*B*>0). In the former scenario, equation (2.1) becomes null, and equation (2.2) is precisely as considered by Laio *et al.* [8] for the soil moisture equation. In the latter scenario, *g*(*B*)>0 and may be used to divide equation (2.1). Since *η* given in terms of *X*=*S*−*S*_{w})
*η*(*S*)≤1, this can be satisfied only for *β*/*α*<1; the biomass growth rate must sometimes exceed the loss rate for a vegetated state to persist. Meanwhile, taking the expectation of equation (2.2) at steady state (written for *X*) gives
*B*.

### (a) The joint distribution

The soil moisture probability density *f*_{B,X}(*b*,*x*,*t*) of the soil moisture and biomass system will evolve according to the Kolmogorov equation, which gives the rate at which probability accumulates around a given point (*b*,*x*). The derivation of this equation corresponding to the soil moisture component only is given by Rodriguez-Iturbe & Porporato [9] and Rodriguez-Iturbe *et al.* [17]. The corresponding equation for the coupled biomass/soil moisture system is very similar; an advection term for *B* corresponding to the continuous evolution generated by equation (2.3) will appear in addition to the advection and jump terms for the soil moisture, yielding (see also [18])
*ρ*_{B}(*b*,*x*)=−*g*(*b*)(*αη*(*x*)−*β*) and *ρ*_{X}(*b*,*x*)=*h*(*b*)*η*(*x*) give the continuous loss rates of variables *B* and *X*, respectively, as defined by equations (2.1) and (2.2); these represent the advection of probability as the system evolves absent rainfall infiltrations. The two terms preceded by λ correspond to the infiltration jump process. In particular, the first of these gives the rate at which trajectories jump away from state (*b*,*x*), whereas the second gives the rate at which trajectories jump into it; the sum of the pre-infiltration soil moisture state and the random infiltration depth (with density *f*_{I}) results in a probability space convolution with respect to variable *x*. The soil moisture is nominally unbounded in the above equation, in that its trajectories are permitted to exceed full saturation (i.e. *S*=1 or *X*=1−*S*_{w}), but in the arid and semi-arid cases of interest, such trajectories, which represent saturated run-off, occur with small probability and neglecting the upper bound has minimal effect on the results, particularly for the biomass. The steady-state condition is realized by dropping the time dependence, so that the left-hand side vanishes. When the explicit definitions of the loss functions are substituted, and exponential distribution of the marks is used in the convolution term, the equation becomes

In general, it would be difficult or impossible to obtain exact solutions to this equation, but the chosen forms of the loss functions allow for a solution to be obtained in this case, and that solution turns out to be separable: *f*_{B,X}(*b*,*x*)=*f*_{B}(*b*) *f*_{X}(*x*), as derived below. Separability implies that *X* and *B* are statistically independent; knowledge about the state of one reveals nothing about the state of the other. This is not an obvious result, but it can be made sense of qualitatively. Suppose a below-average soil moisture value is observed, and consider the alternative system trajectories which could generate such points. Short dry spells occur relatively frequently which suppress the soil moisture significantly and the biomass slightly, but the soil moisture tends to rebound as soon as a few rains occur. A long wet spell occurs with low frequency, but the accumulated biomass persists for a long time, suppressing the soil moisture in an ongoing fashion. Thus, the low soil moisture observation corresponds in the first case to a slightly reduced biomass with high probability, or in the second case to a substantially increased biomass, but with low probability. The separability result indicates that these two effects tend to cancel each other out not only approximately, but exactly.

To prove that this holds, return to equation (3.4), which upon rearrangement becomes
*b*, this becomes
*g*(*b*)*f*_{B}(*B*) must vanish at the integration bounds to ensure normalizability. Equation (3.6) reduces to the Kolmogorov equation for the soil moisture process with 〈*h*〉 as a parameter, and its solution is shown by Laio *et al.* [8], Rodriguez-Iturbe *et al.* [17] to be equivalent to
*f*_{B}(*b*) and substituting it into equation (3.5) yields (dropping the functional dependencies for brevity)
_{x}[*ηf*_{X}] and regrouping terms gives
*b*, but the first is multiplied by *η*(*x*), so the above can only hold if each term vanishes independently
*hη*〉=〈*h*〉〈*η*〉. Thus, the separable solution holds.^{3} Using the expression for 〈*h*〉 above in equation (3.7) and either equation (3.10) or equation (3.11) gives the two differential equations which determine the distributions *f*_{X} and *f*_{B}, respectively, as well as the integrated form of those equations.
*g*(*b*)∼*b* when *b* is small, it must be the case that *f*_{B} normalizable. It follows then from equation (3.14) that this will only hold when *h*(0^{+})<〈*h*〉=λ*α*/*βθ*. This innocuous statement has deep implications. Since the assimilation rate vanishes for *b* near zero, *h*(*b*) in that region corresponds to other water losses (e.g. evaporation), and if *h*(*b*) is too big, even if only in a small neighbourhood near *b*=0, biomass trajectories will get trapped there, unable to take a big enough share of any future precipitation to escape with certainty. In this case, the steady state will have zero biomass.

Attention now turns to the specific form of the dynamics. The biomass distribution may be determined by taking *h*(*b*)=*γb* and *g*(*b*)=*b* according to equations (2.3) and (2.4), which are the forms used in the following sections. Since *h*(0)=0, the steady state will be vegetated, and it is seen to follow a gamma distribution with shape parameter *q*=λ/*β* and rate parameter ^{4} The density function and moments are given by
*κ*_{n} are also introduced in the above, as these prove more natural to work with later; the first and second cumulant are the mean and the variance, and the higher order cumulants can be related to the moments in a straightforward way (e.g. [19]). It should be noted that with respect to the first two moments or cumulants (i.e. the mean *μ*_{B}≡*κ*_{1} and variance *et al.* [14], which was obtained by approximating the biomass evolution as a diffusion process. However, the shape of the distribution can be quite different when the shape parameter *q* is small; the distribution looks approximately Gaussian only for *q*>>1. This distinction is illustrated in figure 1.

The soil moisture distribution is independent of *h*(*b*) and *g*(*b*); taking *η*(*X*) as given in equations (2.3) and (2.4), the density function and moments are (in terms of 〈*h*〉 and its explicit definition 〈*h*〉=λ*α*/*βθ*):
*et al.* [8] and Rodriguez-Iturbe *et al.* [17], corresponding to an effective loss function of 〈*h*〉*η*(*x*). Because *h*(*b*)=*γb*, it follows that 〈*h*〉=*γμ*_{B}, i.e. the transpiration rate (the only source of moisture loss in this system) associated with the mean steady-state biomass.

### (b) Steady-state biomass: an example

To provide some graphical intuition for the distributions described above and to illustrate the difference between the exact (gamma) and approximate (Gaussian) distributions, it is worth considering a concrete example. Following Nordbotten *et al.* [14], the values of the necessary parameters for a typical semi-arid system are given in table 1, and these are used to generate plots of the steady-state biomass and soil moisture distributions in figure 1.

The gamma shape parameter derived from these values is λ/*β*=21.3, which, while rather large compared with unity, still results in a significant discrepancy between the exact and approximate density functions, especially with respect to the extreme values, far from mean. In dry climates with a smaller λ, or under conditions with a high loss rate, these differences become even more pronounced.

The soil moisture distribution, presented for the translated variable *X*=*S*−*S*_{w}, is exact with respect to the system defined by equations (2.4) and (2.3), and even though leakage was neglected in that formulation, the relatively low values of the resulting soil moisture suggest that this will not be very important; the problem of estimating the error associated with neglecting leakage is discussed in the next section in the context of the time-dependent distributions.

## 4. Time-dependent analysis

The time-dependent Kolmogorov equation does not admit a simple solution as in the steady-state case, since the state variables cannot be separated from each other or from time. However, a useful approximation to the time-dependent distribution is derived in what follows. The linearity of the loss term in equation (2.3) allows it to be combined with the time derivative via the integrating factor e^{βt}, and substitution for the assimilation term from equation (2.4) gives
*t*≤*T* gives
*β*. The second term is the biomass which would be added due to all infiltrations up to time *T* provided all the water of each infiltration were immediately transpired. The third term accounts for the accumulation of water in the soil and tends to reduce or increase the final biomass depending on whether the soil moisture is generally increasing or decreasing, i.e. whether more water is stored from new infiltrations than is withdrawn from the initial supply.

As will be shown in what follows, the first two terms will dominate the third term for the typical parameters associated with the dry land systems of interest; when the third term is neglected entirely, the system is said to be in the no-storage limit, where all water (above the wilting point) is transpired as soon as it infiltrates. More generally, this term will give rise to a correction to the no-storage limit, whose magnitude and behaviour can be estimated from the baseline properties of the no-storage case for the biomass, combined with the steady-state analysis of the previous section for the soil moisture.

### (a) No-storage limit

Deferring temporarily the discussion of under what parameter conditions the no-storage limit might apply, it is worth describing its behaviour. This limit is governed by equation (4.1) or equation (4.2), with the term d*X*/d*t* set to 0, so that the biomass is driven directly by infiltration:
*B* can be obtained from the integral formulation above by using the accumulation properties of marked Poisson processes to characterize the infiltration term, but it is simpler to proceed from the Kolmogorov equation as previously.^{5} Defining *W*=e^{βt}*B* and *H*(*t*)=(*α*/*γ*) e^{βt}*I*(*t*), the density function *f*_{W}(*w*,*t*) obeys
*w*. Since *H* is a rescaling of *I*, the LTs are related by *I* (which has *B* (via the scaling relation *B* in LT space is due to the representation of *B* as the sum of independent random variables, as in the integral form of equation (4.3). The first term corresponds to the distribution of the initial value *B*(0), namely *B*(0) may be a delta function, if for instance the initial value is fixed by a definite observation, but *B*(0) may also be a non-trivial random variable with a distribution corresponding to either measurement uncertainty (i.e. in the case of a particular system), or due to stochastic influences on the biomass evolution prior to the initial time of the interval [0,*T*] under consideration (i.e. for an ensemble of systems). The logarithm of the transformed distribution is the generating function of the cumulants *κ*_{j}(*B*), each of which is a sum of the contribution from the initial state distribution and the biomass assimilated from the infiltrations up to time *T*:
*et al.* [14]:
*B*(*T*) can be obtained from the (e.g. numerical) inverse Laplace transform of equation (4.6) for any initial distribution. However, it is also possible to isolate the distribution for the integral term in equation (4.3) which is proportional to the variable *Y* (*T*) just defined. The result is
*T*], and this term can usually be neglected for timescales of interest, which are on the order of a season. The biomass *B*(*T*) is related to *Y* (*T*) by a translation and a scaling, so that upon computation of the density in equation (4.10), the density *f*_{B(T)} may be recovered according to
*T* goes to infinity, the distribution of *Y* (*T*) takes a simple form, namely that of a gamma distribution:
*X*/d*t*) is added to it, the correlation structure between them is precisely such that the distribution is preserved.

### (b) Storage and leakage effects

It was previously suggested that the soil moisture storage embodied by the last term in equation (4.2) would typically result in only a small correction to the biomass distribution. Proving this requires a way of estimating the soil moisture distribution, and it is argued here that this can be done by exploiting the time-scale separation (observed by Nordbotten *et al*. [14]) between the soil moisture and the biomass; in particular, the soil moisture adjusts (i.e. tends towards equilibrium) on a timescale much faster than the biomass.

Equation (4.2) shows that the effect of the initial biomass decays according to e^{−βT}, so that the natural correlation timescale of the biomass process is
*X*=0. In particular, the fractional change in soil moisture due to the loss function sets the timescale on which this occurs:
*a priori* that the typical value of *X* is no more than *k*, and in fact it will be shown subsequently to be significantly less. If the time-scale ratio is much less than one, the soil moisture will achieve a quasi-steady state around the local biomass process. The steady-state distribution of the soil moisture was shown to be independent of the biomass value and depend only on its steady-state mean. The transient distribution of the soil moisture may be approximated by the same distribution with the time-dependent mean biomass *μ*_{B}(*T*), as determined by the no storage limit in equation (4.8), replacing the steady-state value *μ*_{B} as a parameter in equation (3.15).
*μ*_{B}(*T*), is close to the true value, 〈*B*(*T*)〉. Confirmation may be obtained by taking the expectation of equation (4.2) after integrating its last term by parts
*ε*_{B}(*T*). The assumption that the no-storage limit dominates the biomass distribution and that accordingly equations (4.17) and (4.18) are valid, is established by confirming that the error term is relatively small, |*ε*_{B}(*T*)|≪*μ*_{B}(*T*).

The interpretation here is that the steady-state soil moisture results, e.g. of Laio *et al.* [8] or Rodriguez-Iturbe *et al.* [17], can still apply, even when the soil moisture is coupled to a varying biomass which has not yet reached steady state. This is because the biomass varies on a much longer timescale, and so looks approximately constant from the view of the soil moisture process. Consequently, the same procedure used in this section can be used to estimate the sensitivity of the results to other forms of soil moisture loss. For instance, if a leakage term *L*(*X*) were added to equation (2.4), the soil moisture distribution at time *T* would be approximately the one induced by the steady-state properties of

### (c) Time-dependent biomass: an example

As a demonstration of the accuracy of the preceding approximations, the time-dependent analogue of the example in §3b may be considered. Again following Nordbotten *et al.* [14], the soil, vegetation and precipitation parameters are furnished by table 1. For these parameters, the typical time-scale ratio can be evaluated using the steady-state mean *μ*_{B} from equation (3.15) in place of the biomass, so
^{6}

The small deviations which do exist between the simulated and analytical results have physical interpretations. Consider the initial condition of the first plot of figure 2, in which the system starts with steady-state soil moisture but less than steady-state biomass. The initial rate of transpirative loss due to the reduced biomass will be accordingly reduced, so that the soil moisture will tend to increase initially, and the subtracted soil moisture term of equation (4.2) will be positive. As a result, the no-storage limit slightly overpredicts the biomass at early times, as seen in the plot for *T*=70 d, and this may be understood as the failure to immediately realize all of the potential new biomass of the incoming precipitation, some of which is instead converted to stored water. Similar reasoning explains the slight underprediction of no-storage limit when the system starts with inflated rather than reduced biomass, as in the third plot of figure 2. The size of these discrepancies can be described quantitatively using the mean correction term *ε*_{B}(*T*) given in equation (4.19), and this analysis is carried out in the next section for the case of climatic fluctuations.

## 5. Climatic fluctuations

In §4a, the climate and plant physiological parameters were taken to be constant over the integration time *T*. This assumption may be dropped, whereby the integrating factor in equation (4.3) becomes *T*) for time-varying parameters is then
*J* epochs, with time *T*_{j} spent in each so that the total evolution time is *B*(*T*). If each epoch *j* has its own constant parameters *α*_{j},*β*_{j},λ_{j},*γ*_{j},*θ*_{j}, then carrying out the integral in equation (5.1) and exponentiating both sides gives the result
*X*_{j}(*t*_{j}) is the soil moisture state at time *t*_{j} after the start of epoch *j*. The time *T*_{j} of each epoch causes the existing biomass to decay by e^{−βjTJ}, while new assimilation, given by the parenthetical term in the sum, replenishes it. The longer ago a unit of biomass was assimilated relative to the observation time *T*, the more it has decayed. If the integral soil moisture storage term is neglected, all that remains are the mutually independent no-storage-limit assimilations *Y*):

### (a) Seasonal fluctuations

An important case of distinct climate epochs arises when each epoch corresponds to an annual season, with each year comprising a wet and dry season. This scenario can be modelled by making the rainfall parameters seasonally dependent, so that λ∈{λ_{w},λ_{d}} and *θ*∈{*θ*_{w},*θ*_{d}}, where the wet and dry parameters apply, respectively, to seasons of length *T*_{w} and *T*_{d}, and where the season durations must sum to the annual duration *T*_{w}+*T*_{d}=*T*_{a}=365 days. Suppose for concreteness that each year *j* consists of a dry season followed by a wet season, the initial biomass at the start of year *j*=1 is distributed as an arbitrary variable *B*_{0}, and that the biomass evolves under these dry/wet fluctuations into year *j*=*J*. Then the quantity of interest is the biomass after time *t* has elapsed in the dry season, *B*_{J,d}(*t*), or the wet season, *B*_{J,w}(*t*) of this final year *J*.

The LTs of the density functions follow from equation (5.2), though it is convenient to work with the scaled biomass *p*=*p*(*t*)=e^{−βt}, and *p*_{(⋅)}=*p*(*T*_{(⋅)}) for (⋅)∈{*a*,*w*,*d*}. For either the dry or wet LT above, the corresponding biomass density function can be obtained through the inverse LT, which may be implemented numerically, e.g. via Talbot's method. Since all of the singularities lie along the negative real axis, there are very efficient techniques for performing this inversion [21]. The inversion may be performed with respect to the variable *B*′, and the distribution in the physical biomass units can be recovered by rescaling
*θ*_{w}=*θ* and λ_{0,w}=λ_{0}. A dry season can be defined with a mean rainfall depth and rainfall arrival rate which are half of the wet season values, so that *θ*_{d}=2*θ*_{w} and λ_{0,d}=λ_{0,w}/2. The plant parameters governing transpiration, assimilation and loss are taken to be the same for both seasons. A indicative measure of the system performance is the biomass distribution at the end of each wet or dry season, which is presented in figure 3. For both the wet and dry cases, the initial value *B*_{0} is taken to follow the corresponding wet or dry steady-state distribution (equation (3.15)), so that the initial distributions (in scaled units) correspond to

With the given parameter choices, figure 3 shows that the biomass distribution induced by seasonal fluctuations at the end of a wet or dry season is significantly different from the distributions of their respective steady states, which would tend to prevail if the seasons were made very long. However, the effects of the alternating regime are well established after 1 year, regardless of whether the system was initially prepared in the dry or wet state, and continued exposure to the seasonal climate has only a small impact. If a seasonal steady state is defined as corresponding to an infinite number of years spent alternating between the wet and dry seasons (i.e. *J*. Consequently, the seasonal steady-state mean and variance are useful descriptors of the system and can be obtained from the general cumulant expressions in appendix C, and the subscript *J* has been dropped since every year becomes statistically identical in the limit.
*p*_{(⋅)}=e^{−β(⋅)T(⋅)}, (⋅)∈{d,*w*,*a*}, which are the decay factors, respectively, per dry season, per wet season and per annum. These means and variances correspond to the end of their respective seasons and to the case where only the climate parameters vary, as in the example that follows.

### (b) Comparison with simulation results

The analytic results in figure 3 were obtained for the no-storage limit. As in the previous section, these results are checked against the distributions obtained by direct simulation of the dynamic equations, and simulated distributions are also plotted there. The simulations were performed both with and without leakage losses from the soil, but the results were virtually indistinguishable; the results presented here are for the no-leakage case. Once again, the no-storage limit provides very good agreement with the simulated results, though close inspection of the figure shows that while the shape of the analytic distribution is extremely accurate, there is slight tendency of the no-storage limit to underpredict the mean value, even after several years. The error in the end-of-season mean, defined as the exact, simulated average *μ*_{{d,w}}, are given below for the wet and dry season in the third year, at which point the seasonal climate is well established:
*J* increases in equations (C1) and (C2). This suggests that system has reached a seasonal steady state, whose properties are given to high accuracy by the long-time limit *ε*_{B}. When this is evaluated at the end of the dry or wet season, the result is
^{7} As previously, the soil moisture expectations are specified from equation (4.18) using the no-storage estimates of the mean biomass, which now refer to the steady bi-seasonal regime. The expression for *μ*_{{w,d}}, the end-of-season mean biomass, is given by equation (5.9). The mean biomass at time *t* within a given season follows most easily from equation (4.8), evolving from the value at the end of the preceding season:

## 6. Discussion

The model presented here describes the dynamics of non-woody biomass and soil moisture. The steady-state joint distribution of the biomass and soil moisture was derived exactly for the case of a system in which leakage and run-off are negligible, and this distribution was shown to be separable, so that the biomass and soil moisture are statistically independent in the long-time limit. Although the primary distribution of interest here was unimodal, corresponding to no soil evaporation, it was shown that the steady-state results are sufficiently general to generate more complex distributions.

When the dynamics were taken to be those of a landscape limited only by water, so that the assimilation, transpiration and biomass loss rates were linear with respect to the biomass, approximate transient distributions of the biomass and soil moisture were also derived. These solutions exploited the time-scale separation between the slowly varying biomass process and the quickly varying soil moisture process, and demonstrated that the biomass dynamics in the limiting case of no soil moisture storage (i.e. instantaneous transpiration of rainfall) are a good predictor of the full system behaviour. The no-storage biomass distribution was then generalized to allow for time-varying parameters, in particular to account for the seasonal behaviour in the precipitation.

Several challenges remain to be explored in this area, in particular those related to the scale and interpretation of the dynamics. The model here involves the area-averaged biomass, and the appropriate averaging length scale must be determined carefully. For example, the averaging area must be small relative to the typical storm cell in order for the point-process description of rainfall used here to be correct. Once the averaging scale is chosen, another problem arises in upscaling from a landscape description which resolves the transpiration features of individual plants (e.g. [22]) or even components of intra-plant hydrology, such as in the whole-tree conductance model given in Manoli *et al.* [23]. Such upscalings can be tricky, as the spatially averaged dynamics will not always be exactly defined in terms of the spatially averaged variables (for instance, as noted by Nordbotten *et al.* [24], the average transpiration function will not in general be a unique function of the average soil moisture), but they are necessary to maintain a solid physical grounding for the sort of analysis given here.

It is not intended to suggest that the analytical framework developed in this work will capture the nuance of every landscape. Rather, the objective has been to describe the baseline behaviour of water-limited ecosystems in a way which is simple enough to be tractable, but nevertheless possesses the explanatory richness and transparency to tackle some important and difficult topics in ecology and ecohydrology. Chief among these are vegetation stability, especially in the face of variable climates, competitive invasion by different species and long-term evolutionary behaviour. It is hoped that the results derived here will help to characterize and address some of these problems.

## Data accessibility

There are no data associated with the results, which are theoretical in nature.

## Authors' contributions

The authors all conceived and developed the results herein collaboratively, and jointly prepared the manuscript.

## Competing interests

The authors declare that there are no competing interests.

## Funding

B. E. Schaffer's research at Princeton University is supported by The Mary and Randall Hack '69 Graduate Award (through Princeton Environmental Institute) and NSF Training Grant in the mathematics of water (through Princeton Program in Applied and Computational Mathematics).

## Appendix A. A bimodal distribution

As previously remarked, the steady-state biomass distribution for the specific dynamics of equations (2.3) and (2.4) is always unimodal. The recent work of Cueto-Felgueroso *et al.* [11] calls attention to the important possibility of an attractor at *B*=0 which makes the distribution bimodal, and accordingly it is worth giving an example in which equations (2.1) and (2.2) do generate this feature.^{8} Keeping *g*(*b*)=*b* as previously, define *h*(*b*)=*γg*(*b*)+*E*(*b*), so that it encapsulates the original transpiration part *γg*(*b*), but also includes a dedicated evaporation term. Equation (3.14) can be written for this case in differential form as
_{b}*f*_{B}(*b*). If *E*(*b*) starts sufficiently large and decays sufficiently quickly, this sign can change twice and a new mode can emerge. As an example, suppose *E*(*b*)=*E*_{0}e^{−cb}, following a Beer's law dependence [25]. The condition
*B*=0, in the sense that ∂_{b}*f*_{B}(0^{+}) is negative, and is consistent with the vegetation existence requirement (§3) that *E*_{0}<〈*h*〉=*α*λ/*θβ* for suitable choice of parameters. For sufficiently large *b*, ∂_{b}*f*_{B}(*b*) is also negative, so if it switches sign in the interim, there is another distinct mode. This is guaranteed if the function on the right-hand side of equation (A1) has a maximum on the interior of its range at which it achieves a positive value, which is equivalent to the two conditions
*c* sufficiently large provided λ/*β*>1. Thus, the absence of a bimodal distribution from the examples in the preceding sections is due the particular choices for *g*(*b*) and *h*(*b*), which were appropriate under the given assumptions, and necessary for the developments in §§4 and 5. However, as shown above, the form of equations (2.1) and (2.2) is sufficiently general to generate bimodality for the steady-state case, as demonstrated for a common evaporation model. A more detailed analysis of the correspondence between functional forms of *h*(*b*) and *g*(*b*), parameter values and behaviour of the modes is deferred to a future work.

## Appendix B. Time-dependent distribution: inverse transform

The LT of the distribution of newly assimilated biomass in the no-storage limit is given by the second term in equation (4.3). In dimensionless units, this is the quantity *Y* contains a kernel of probability (i.e. a delta function) at *y*=0, corresponding to the realization of the infiltration process where there are no infiltrating rainfall arrivals in the interval [0,*T*], which happens with probability e^{−λT}, and so the density function may be broken up into its continuous and kernel parts:
*q*=λ/*β* and *p*=e^{−βT}. The inverse Laplace transform above can be carried out numerically using, e.g. Talbot's method, which along with the delta function term provides a way of computing with the time-dependent distribution. However, it is also possible to obtain an analytical result for this inverse transform; this result is not found in the standard references, and so is derived here. To carry out the inverse Laplace transform in equation (B1) analytically, first translate the variable *z*→*z*−1 so that the integral of interest is
*p*)/*p*,0) and the denominator has branch point at (0,0). The desired quantity is

Since the contour contains no poles, Cauchy's theorem gives
*C*_{2} and *C*_{6} will vanish due to the exponential term for any *y*>0; the delta function at *y*=0 is the same as previously discussed. The integrals along the negative real axis have to be carefully parametrized. In particular, because of the branch cut, the integrand will be different on each segment and they will not cancel. On *C*_{3}, the various components of the integrand can be parametrized according to the distance *l* along the axis:
*C*_{5} gives
*l*∈[(1−*p*)/*p*,*R*], leaving
*C*_{4} is *z*=*ϵ*e^{iθ} for *θ*∈(−*π*,*π*), so
*ϵ*, then upon dividing by the term e^{q−1} in the denominator, there will be *n*≡⌊*q*⌋ terms which have *ϵ* to a negative power. The rest of the terms will vanish as *ϵ*→0. The *n* non-vanishing terms would tend to blow up, but they are in fact precisely cancelled by the *n* boundary terms when the first integral is integrated by parts *n* times. The only part of the whole expression which remains is the (convergent) integral left by the integration by parts (introducing also the remainder of *q* as *r*≡*q*−*n* and *w*≡*l*/*k*):
*w*=0) separately by either subtracting and adding *w*^{−r} or e^{−ykw}*w*^{−r}, as shown in the first and second expressions below, respectively. The first expression is ideal for the computations of interest, when 0<*T*< ∼1/*β*, while the second is included to make the convergence to the steady-state gamma distribution more evident when *p*→0 and *j*=0 term survives, the integral term vanishes and the lower incomplete gamma function *γ*(1−*r*,*ky*)→*Γ*(1−*r*). This expression facilitates estimation of error when the limiting distribution is used instead of the full time-dependent distribution.

## Appendix C. Cumulants of the inter-seasonal distribution

Following from §5a, the no-storage biomass distributions after time *t* has elapsed in the dry or wet season of year *J* have cumulants derived from equation (5.4). The sum gives rise to two geometric series, one for the dry and wet season, respectively, each with a geometric multiplier *p*_{a}. Performing the summation gives
*Y* _{{d,w}} is the (non-dimensional) biomass accumulation variable for a (dry or wet) season. Its cumulants at time *t* in the season are given by equation (4.7), with subscripts now assigned to distinguish the seasonal parameters:
*p*_{(⋅)}=e^{−β(⋅)T(⋅)}, (⋅)∈{d,*w*,*a*} as given in §5*a*. When the seasonal steady state applies (*n*=1,2, respectively.

## Footnotes

↵1 These dynamics are written in continuous time, but should be interpreted at the daily timescale, for which they provide a very good representation [7].

↵2 The results of Cueto-Felgueroso

*et al.*[11] appeared while this paper was under revision, and so it seems natural to include a complementary analytic introduction to the idea of bimodality, with more extensive analysis deferred to a future work.↵3 This result is due to the precise forms of equations (2.1) and (2.2), the key features being the soil moisture loss function is a product of a biomass function

*h*(*b*) and a soil moisture function*η*(*X*); the biomass net gain rate also has a product form, which allows the logarithmic equivalent used in equation (3.1); the resulting logarithmic loss rate is constant; the evapotranspiration and assimilation soil moisture dependencies are the same, namely*η*(*X*), which implies that the water-use efficiency is constant with respect to soil moisture. The solution also requires the exponential distribution for the infiltrations.↵4 The emergence of a second mode due to evaporation is considered in appendix A by using different forms of

*h*(*b*) and*g*(*b*).↵5 Viola

*et al.*[20] give a derivation for the transient soil moisture dynamics similar to what follows, although the final solution is presented in different form.↵6 Simulations incorporating leakage losses were also performed using

*L*(*x*)=(*K*_{S}/*nZ*_{r})(*x*+*S*_{w})^{c}with*c*=10 and*K*_{S}=2000 mm d^{−1}[4], though the results were indistinguishable from the no-leakage case for the given parameter values.↵7 In the seasonal steady-state limit, each wet or dry season is statistically indistinguishable from the one occurring in the previous year.

↵8 Here, bimodal is taken to mean simply that the distribution has two peaks; see Cueto-Felgueroso

*et al.*[11] for the connection to the associated attractors of the deterministic system.

- Received March 11, 2015.
- Accepted October 19, 2015.

- © 2015 The Author(s)