## Abstract

Feedbacks between water use, biomass and infiltration capacity in semiarid ecosystems have been shown to lead to the spontaneous formation of vegetation patterns in a simple model. The formation of patterns permits the maintenance of larger overall biomass at low rainfall rates compared with homogeneous vegetation. This results in a bias of models run at larger scales neglecting subgrid-scale variability. In the present study, we investigate the question whether subgrid-scale heterogeneity can be parameterized as the outcome of optimal partitioning between bare soil and vegetated area. We find that a two-box model reproduces the time-averaged biomass of the patterns emerging in a 100 × 100 grid model if the vegetated fraction is optimized for maximum entropy production (MEP). This suggests that the proposed optimality-based representation of subgrid-scale heterogeneity may be generally applicable to different systems and at different scales. The implications for our understanding of self-organized behaviour and its modelling are discussed.

## 1. Introduction

Heterogeneities such as periodic patterns can be observed at different scales in a range of natural and experimental systems. Examples include waves on water, sand or snow surfaces, convection cells in turbulent flow experiments, wave-like clouds behind mountain ridges, hummocks and hollows in peat bogs and vegetation patterns in semiarid vegetation, the latter of which are the subject of the present paper. For more examples, see Rietkerk & Van de Koppel (2008).

Whereas the shape and dynamics of the different patterns is an interesting research area in itself, this paper focuses on the general problem of distributed numerical models to account for the heterogeneity associated with such patterns, often synonymously referred to as subgrid-scale heterogeneity or subgrid-scale variability. Subgrid-scale variability is related to heterogeneity that is not resolved in a numerical model because the grid cells are larger than the elements of the underlying pattern. The neglect of subgrid-scale variability can lead to substantial model bias (Larson *et al.* 2001). On the other hand, setting the resolution of a model high enough to resolve all underlying patterns is usually not feasible especially in large-scale models.

Different mechanisms that could lead to the formation of periodic patterns have been explored in the literature and a range of numerical models that generate patterns have been developed. Klausmeier (1999) developed a simple model capable of simulating the formation of vegetation patterns in response to homogeneous forcing. The model (henceforth referred to as the ‘Klausmeier model’) is based on two coupled differential equations, one for water and one for biomass. It is based on the assumption that the existence of biomass permits increased infiltration and hence increased water use, while increased water use in return permits increased biomass growth. This leads to a positive feedback between biomass and transpiration. A description of the Klausmeier model is given below in §2. Downhill flow of water and diffusion-like dispersal of biomass provide the lateral feedbacks necessary for spatial organization in the distributed model. This can lead to the spontaneous formation of vegetation stripes that follow the contour lines of the hillslope and slowly move uphill. Klausmeier (1999) demonstrated that a realistic choice of a few model constants leads to a reproduction of the observed wave length of vegetated stripes in semiarid tiger bush and the observation that the wave length increases with decreasing rainfall or increasing mortality (e.g. due to grazing). Depending on rainfall, mortality and the initial conditions, the system can take on a homogeneously vegetated steady state, a patterned steady state or a bare soil steady state. Using a similar model, Lejeune *et al.* (2002) pointed out that the average biomass of the patterned state can be higher than the biomass of the alternative homogeneous state. This suggests that the neglect of pattern formation can lead to a model bias towards systematical under-prediction of the average biomass.

More recently, other models such as cellular automata (Rietkerk *et al.* 2004*a*) or additional processes such as the reduction of bare soil evaporation by wilted biomass (Zeng *et al.* 2004) or nutrient transport (Rietkerk *et al.* 2004*b*) were explored to simulate vegetation patterns other than stripes. It appears that many different mechanisms can lead to the emergence of patterns and many plausible answers can be found to the question *how* the patterns form. Moreover, the quantitative effects of the proposed mechanisms are largely unknown and hence the results of the mechanistic models are uncertain and possibly biased. However, given that the proposed mechanisms themselves are adaptive properties of the system (e.g. seed dispersal strategies, facilitation and inhibition between plants), the question arises *why* the ecosystem produces regular patterns under certain conditions. The answer to the latter question could permit the formulation of models that simulate the effects of the patterns without resolving the exact mechanism causing the patterns, i.e. large-scale models that account for the effect of subgrid-scale variability, independent of initial conditions.

One possible explanation for self-organization and the formation of patterns is optimal adaptation to fulfil a certain objective function of the system (Schymanski 2008). The finding that the patterned vegetation state has a higher biomass than a homogenous distribution for a given rainfall (Lejeune *et al.* 2002) suggests that pattern formation may be an adaptation to increase the biomass in the system. However, the question remains whether the maximization of biomass does indeed explain the formation of vegetation patterns and how universal such an objective function is. A range of different objective functions have been proposed for biological systems, mostly relating to the carbon cycle or to plant water stress (Schymanski 2009*a*,*b*). On the other hand, several authors suggested that the maximum entropy production (MEP) principle could be a universal objective function for self-organized systems with sufficient degrees of freedom, such as atmospheric, biological or hydrological systems (Lorenz *et al.* 2001; Dewar 2003; Ozawa *et al.* 2003; Kleidon 2004; Kleidon & Schymanski 2008). The MEP principle can be briefly explained as follows.

The second law of thermodynamics tells us that all spontaneous flows and reactions lead to an increase in entropy (d*S*) in an isolated system. This results in the predictability of the direction of chemical reactions or flow of matter (e.g. water flows downhill). In more general terms, all spontaneous reactions and flows can be formulated as functions of gradients, also called ‘thermodynamic forces’ (Kondepudi & Prigogine 1998). Examples are temperature gradients driving heat flow, pressure gradients driving the flow of momentum and chemical potential gradients driving flow of matter and chemical reactions. In these cases, entropy increases owing to the depletion of the gradient driving the flow or reaction.

In open systems, such as biological or hydrological systems, changes of entropy are a result of exchange of energy and matter with the exterior and of irreversible processes in the interior of the system:
1.1
where d_{e}*S* and d_{i}*S* denote the entropy exchange across the system boundary and the entropy production due to irreversible processes inside the system, respectively (Kondepudi & Prigogine 1998). Open systems can thus maintain their steady state (d*S* = 0) by balancing their internal entropy production by their entropy exchange with their surroundings (d_{e}*S* = –d_{i}*S*). This allows them to maintain a highly organized state of low entropy as long as they are able to export entropy to their surroundings. Isolated systems, in contrast, only reach steady state at thermodynamic equilibrium when they have maximized their own entropy.

The entropy production inside the system is the sum of all internal fluxes (d*X*_{i}) that are driven by thermodynamic forces (*F*_{i}) (Kondepudi & Prigogine 1998):
1.2
Neglecting temperature and pressure gradients and considering only gradients in the chemical potential of water, equation (1.2) can be written as (Kleidon & Schymanski 2008):
1.3
where *μ*_{1} and *μ*_{2} are the chemical potentials of water at two different locations within the system, *T* is the temperature of the system and d*W* is the flux of water between these locations. The chemical potential of soil water is a nonlinear function of the soil water content, which is usually derived from water retention curves.

Figure 1 illustrates the entropy exchange between two grid boxes in the Klausmeier model due to the exchange of water and the entropy production within each grid box due to the dissipation of gradients created by the fluxes. Fluxes of water dissipate gradients between boxes but create gradients within the receiving boxes because they add water with higher chemical potential to the existing water. This results in a chemical potential gradient (thermodynamic force) within the receiving box, which is dissipated by the mixing of the incoming water with the water already present. The MEP principle suggests that the system would adapt its degrees of freedom (e.g. the conductivities to flow by means of its biomass distribution) in such a way as to maximize its entropy production. According to figure 1, two processes lead to the entropy production within the system: the mixing of rain water with the water in each box and the mixing of water from different boxes. It is obvious that vegetation patterns would increase the entropy production due to mixing of water between the boxes, as patterns imply heterogeneity and heterogeneity implies the existence of gradients within the system. However, the effect of patterns on the entropy production due to infiltration of rainfall is not as obvious and depends on the water content in the boxes.

Although vegetation patterns have been identified as dissipative structures (Lejeune *et al.* 2004) before, they have not been analysed from a thermodynamic perspective to date. The aim of the current study is to test whether the formation of patterns in the Klausmeier model is consistent with the MEP principle and whether this principle can be used to predict the average biomass in the system without the need to know how the patterns form. We wish to stress that the aim of the study is not to evaluate the Klausmeier model and how it relates to natural processes. The Klausmeier model is merely used as an example for a simple dynamic model that reproduces some characteristics of patterns found in semiarid ecosystems. The principles described in this paper could equally be applied to more complex models, as long as their formulation of the dynamics can be translated into thermodynamic terms.

## 2. Material and methods

### (a) Dynamic, distributed model

The Klausmeier (1999) model is based on the premise that plant water use is a function of soil water (*W*), biomass (*B*) and infiltration rate, while the infiltration rate is assumed to be a linear function of biomass, leading to a quadratic dependence of the transpiration rate on biomass. The mass balance for water was thus given as
2.1
The terms on the right-hand side denote precipitation, soil evaporation, transpiration and advection, where *t* is time (years), *W* is the water content (kg m^{−2}), *P* is precipitation (kg m^{−2} yr^{−1}), *B* is dry biomass (kg m^{−2}), *V* is the horizontal flow velocity of water (m yr^{−1}), *X* is the upslope coordinate and *c*_{1} and *c*_{2} are empirical constants. The mass balance for biomass was given as
2.2
The terms on the right-hand side denote growth (constant *c*_{3} multiplied by the transpiration rate), mortality and dispersal of biomass, where *Y* is the coordinate along a contour line of the hillslope, while *c*_{3}, *c*_{4} and *D* are constants representing the water use efficiency, the mortality coefficient and the dispersal coefficient, respectively. For the implementation as a distributed model, the model is usually run with periodic boundary conditions, i.e. water and biomass leaving the domain at the bottom of the hillslope is added at the top of the hillslope and vice versa.

The model represented by equations (2.1) and (2.2) was formulated for semiarid ecosystems (Klausmeier 1999), hence no energy limit was imposed on evapotranspiration. However, due to the concentration of water in the vegetated stripes, evapotranspiration would tend towards infinity as the vegetated surface area tends towards 0 (e.g. for decreasing values of *D* and *P*). To prevent such unrealistic behaviour, transpiration (*E*_{t}, kg m^{−2} yr^{−1}) was set to
2.3
where *a* was set to a value of 4 in order to represent semiarid conditions where the energy limitation of evapotranspiration (*a**P*) is four times larger than annual precipitation. In general, the assumption of semiarid conditions would be equivalent to any value of *a* between 2 and 5 (Arora 2002). Accordingly, soil evaporation (*E*_{s}, kg m^{−2} yr^{−1}) was set to the minimum of the water-limited rate (*c*_{1}*W*) and the energy-limited rate minus transpiration (*aP* − *E*_{t}):
2.4
Equations (2.3) and (2.4) were then substituted for the terms *c*_{2}*W**B*^{2} and *c*_{1}*W* respectively in equations (2.1) and (2.2). With this modification, the distributed model was run on a 100 × 100 grid using different values of constant and homogenous *P*. The model parameters were taken from Klausmeier (1999) and set to the following values: *c*_{1} = 4 yr^{−1}, *c*_{2} = 100 yr^{−1} (kg dry biomass m^{−2})^{−2}, *c*_{3} = 0.003 kg dry biomass (kg H_{2}O)^{−1}, *c*_{4} = 1.8 yr^{−1}, *V* = 365 m yr^{−1} and *D* = 0.5 m^{2} yr^{−1}. To allow the formation of patterns, the model was initialized with a random biomass distribution that varied spatially between 0 and 0.8 kg m^{−2} and a homogeneous water content of 75 kg m^{−2}. For each value of *P* between 150 and 1000 kg m^{2} yr^{−1}, the model was run for 1000 years and the entropy production was calculated for the last year as described below in §2*c*. The time steps of the model were varied such that the state variables could not change by more than 10 per cent in any grid cell in a single time step. Generally, an apparent steady state in biomass was reached after a couple of years, but in some cases, parallel stripes only developed after several hundreds of years. The choice of 1000 years for all model runs ensured that the results stem from fully developed patterns.

### (b) Two-box and one-box models

To simulate the heterogeneity between vegetated and bare surface area in the simplest possible way, we formulated a two-box version of the Klausmeier model, consisting of a vegetated box (denoted by the subscript v) covering an area fraction *A*_{v} and a bare soil box (denoted by the subscript b) covering an area fraction *A*_{b} = 1 − *A*_{v} (figure 2). The same equations were used for the fluxes and the mass balance as in the 100 × 100 grid model, with the following modifications: biomass dispersal was set to 0 in equation (2.2) and *V* (∂*W*/∂*X*) was replaced by *V* (*W*_{b} − *W*_{v})/*Δ**X* in equation (2.1), where we set *Δ**X* = 1 m. This yielded the following mass balance equations for the two-box model:
2.5
and
2.6
for water in the vegetated and bare soil box, respectively, and
2.7
for biomass in the vegetated box. *W*_{v} and *W*_{b} hereby refer to the water content in the vegetated and bare soil box, respectively (kg H_{2}O), and *B*_{v} refers to the biomass in the vegetated box (kg dry biomass). Note that the units are not expressed per unit area as was the case in equations (2.1) and (2.2).

Steady-state conditions for given *P* and *A*_{v} were calculated by setting equations (2.5)–(2.7) to 0 and solving the system of equations for *B*_{v}, *W*_{v} and *W*_{b} using the SAGE Mathematical Software, V. 3.2 (freely available at http://www.sagemath.org).

To simulate homogenous conditions, the two-box model was converted to a one-box model by setting *A*_{v} = 1.

### (c) Calculation of entropy production

In equation (1.2), entropy production was expressed as a function of a flux and the thermodynamic force driving this flux. Inspection of equation (2.1) reveals that all water fluxes in the present model are driven by the water content *W* while *c*_{1} and *c*_{2}*B* can be seen as resistances for the fluxes. To relate *W* to a thermodynamic force (e.g. a gradient in chemical potential), the chemical potential of water in the atmosphere (*μ*_{a}) was set to 0 and the chemical potential of water in the soil (*μ*_{s}) was assumed to increase linearly with *W*:
2.8
where, for simplicity, *c*_{μ} was set to *c*_{μ} = 1 J kg^{−2} m^{2}. This ensures that *W* represents both the chemical potential of water in the soil (*μ*_{s}) and the gradient *μ*_{s} − *μ*_{a}. Note that equation (2.8) is a direct consequence of the implicit linearization of the soil water retention curve in the Klausmeier model.

For every grid cell, we write the entropy balance as (see also figure 1):
2.9
and
2.10
where *Q*_{i} is the incoming flux rate from a neighbouring cell with a water content *W*_{i} while *Q*_{o} is the water flux out of the cell. The temperature was set to *T* = 300 K.

Note that inwards flux of water in figure 1 results in an outward flux of entropy and vice versa (Kondepudi & Prigogine 1998). The chemical potential of precipitation (*μ*_{p}) had to be set to a high enough value so that it was always higher than *μ*_{s} of any of the boxes or grid cells, otherwise the infiltration of rainfall would not always be thermodynamically feasible. For rainfall rates ranging between 0 and 2000 kg m^{−2} yr^{−1} we used a value of *μ*_{p} = 500 J kg^{−1}. Note that this setting is somewhat artificial as the maximum infiltration rate should be a function of the difference in chemical potentials between rain and soil water in a more physically based model, so that a violation of the second law of thermodynamics could not occur. In the present model, the exact value for *μ*_{p} does not affect the overall results as long as it is always greater than the highest *μ*_{s} in each grid cell.

## 3. Results and Discussion

### (a) Entropy production in the Klausmeier model

Entropy production in the Klausmeier model was positively related to the average biomass and negatively related to the average water in the domain (figure 3). Starting from random initial conditions, both biomass and entropy production increase while the patterns begin to form and reach a steady state when the patterns are fully developed (data not shown). This suggests that the self-organization leading to the spontaneous formation of patterns follows the MEP principle. It further illustrates that self-organization leads to an increase of biomass and a reduction of water in the investigated system as compared with the unorganized, homogeneous steady state. The link between pattern formation and entropy production can be understood in terms of the dissipation of gradients. The more organized the system is, the more gradients it contains and the larger is its internal entropy production owing to the dissipation of the gradients. In order to maintain the organized state, the system must be able to export the entropy produced to the surroundings. The system is at steady state only if the entropy exchange with the surroundings (d_{e}*S*) equals the entropy production within the system (d_{i}*S*). Furthermore, ∫(*E*_{t} + *E*_{s}) = ∫*P* at steady state, so that d_{e}*S* is strongest negative if *W* is low and negatively correlated with (*E*_{t} + *E*_{s}) (see equation (2.9)). This is achieved if the biomass is high and heterogeneous, as bare soil patches have a higher *W* but *E*_{t} = 0, while vegetated patches have a lower *W* due to *E*_{t}.

### (b) Optimised two-box model

The finding that the spontaneous formation of patterns in the distributed Klausmeier model follows the MEP principle does not enable us to directly predict the optimal biomass or water distribution in the system, as the possible states of the system are countless. However, the simplified two-box model allows investigation whether a certain value for the vegetated area fraction (*A*_{v}) would lead to MEP. Using the analytical steady-state solutions, we found that for any given rainfall, there is an optimal value of *A*_{v} that maximizes both entropy production and average biomass. Figure 4 illustrates the average biomass in the optimal two-box model as a function of rainfall compared with the average biomass of the distributed Klausmeier model and with the biomass of a homogeneous one-box model. The latter results were simulated by setting *A*_{v} = 1 in the two-box model. The comparison reveals that the optimal two-box model simulates very similar average biomass to the distributed model, whereas the homogeneous one-box model underestimates biomass for low values of rainfall. In fact, the neglect of heterogeneity would lead to a failure of simulating any biomass below a rainfall of 240 kg m^{−2} while the optimized two-box model predicts biomass down to a rainfall of 60 kg m^{−2}. Similar differences between homogeneous and distributed model results have been reported earlier in the literature (Lejeune *et al.* 2002; Dekker & Rietkerk 2005), but the consequence has always been to resort to fine-scale distributed models. The present study suggests that the use of optimality and the MEP principle allows for a simplified representation of self-organized heterogeneity and the simulation of the resulting sensitivity of biomass to rainfall in a simple two-box model.

It is to be noted that the optimal two-box solution leads to slightly higher biomass than the distributed model. This is because the dispersal rate of biomass was set to 0 in the two-box model. In the dynamic, distributed model, biomass dispersal is necessary in order to achieve the organized heterogeneity starting with random initial conditions, but the dispersal leads to a loss of biomass that is neglected in the two-box model. The two-box model also simulates a slightly higher overall entropy production compared with the distributed model, although the entropy production owing to the advection of water is slightly higher in the distributed model (data not shown). The increased overall entropy production is mainly due to the higher overall biomass and consequently decreased overall water content, which has a positive impact on the entropy production by infiltration of rain water.

## 4. Conclusions and outlook

The present study demonstrates the usefulness of the MEP principle for modelling the effects of self-organization in semiarid vegetation without the need to resolve the mechanism by which the self-organization takes place (e.g. biomass dispersal and fine-scale dynamics of water flow). The approach presented here has the potential to be widely applicable to many advection–diffusion systems similar to the Klausmeier model, since it is based on non-equilibrium thermodynamics. The application should also be independent of the scale of the system, as long as the system has sufficient degrees of freedom to self-organize.

## Footnotes

↵† Present address: Department of Bioengineering, Rice University, Houston, TX 77005, USA.

One contribution of 17 to a Theme Issue ‘Maximum entropy production in ecological and environmental systems: applications and implications’.

- © 2010 The Royal Society