Implications

The use of viability theory provides a solid framework to deal with non-equilibrium rangelands where the rancher has different management options (herd dimension and composition, burning/not burning). Viability theory is particularly suited to defining resilience, robustness and management adaptability. In our application, we showed that a pure grazing strategy is viable only for high-rainfall regimes and that management is more adaptable with mixed grazing–browsing herds in states with intermediate values of woody and grass biomass.

Dominant rangeland model paradigms and the concept of resilience

The first model of rangeland dynamics was the range succession model. On the basis of the Clementsian theory (Clements, Reference Clements1916), this model predicted vegetation changes to be continuous and reversible along a gradient (Dyksterhuis, Reference Dyksterhuis1949; Briske et al., Reference Briske, Fuhlendorf and Smeins2003). Following this view, the rancher should find the optimal stocking rate that maintains the rangeland in the climax where grass and livestock production are maximized in the long term (Westoby et al., Reference Westoby, Walker and Noy-Meir1989). However, the reality of rangeland management showed that vegetation changes were not necessarily reversible (Westoby et al., Reference Westoby, Walker and Noy-Meir1989; Laycock, Reference Laycock1991). This means that, if increasing the stocking rate could lead to woody plant encroachment, a successive decrement would not necessarily lead to a symmetrical return to a non-encroached condition. Westoby et al. (Reference Westoby, Walker and Noy-Meir1989) proposed the state-and-transition model, according to which the manager/modeller was expected to define a catalogue of possible states of the system and of possible transitions between states, with the respective causes (Bestelmeyer et al., Reference Bestelmeyer, Brown, Havstad, Alexander, Chavez and Herrick2003; Briske et al., Reference Briske, Fuhlendorf and Smeins2003). The problem of this model was its limited predictive capacity (Ward, Reference Ward2005b; Bashari et al., Reference Bashari, Smith and Bosch2008). In fact, the model could be designed to explore states, transitions, and causes of transitions from experience and observation. However, because not all possible configurations of the system were necessarily considered, the model could not predict the behaviour of the system out of the range of the already observed situations (Ward, Reference Ward2005b).

This problem could be overcome with models based on difference and differential equations (Tietjen and Jeltsch, Reference Tietjen and Jeltsch2007). Describing system dynamics with mathematical laws (i.e. difference or differential equations) rather than with a priori defined states and transitions allowed the modeller to explore all the possible states under all possible values of the system’s drivers. In particular, introducing non-linearity in these models could demonstrate the existence of alternative steady states (May, Reference May1977; Anderies et al., Reference Anderies, Janssen and Walker2002; Beisner et al., Reference Beisner, Haydon and Cuddington2003, Accatino et al., Reference Accatino, De Michele, Vezzoli, Donzelli and Scholes2010), each one characterized by its own attraction basin, that is, the set of initial configurations of the system from which a trajectory would lead to that steady state. It was generally possible to move from one steady state to another, but such a shift required important interventions from the manager. From a non-encroached state, increasing stocking density made it possible to shift to the encroached state. However, starting from an encroached situation, reducing the stocking density might not be sufficient for the system to move back to the non-encroached situation (Anderies et al., Reference Anderies, Janssen and Walker2002). To address the existence of multiple stable states, the concept of threshold (Scheffer et al., Reference Scheffer, Carpenter, Foley, Folke and Walker2001) was particularly appropriate. Thresholds referred to boundaries that, once passed, caused the transition of the system to another steady state, whereas the reverse transition was difficult to accomplish (May, Reference May1977). Because these transitions could be abrupt and irreversible, they were often defined as catastrophic shifts (Scheffer et al., Reference Scheffer, Carpenter, Foley, Folke and Walker2001). In particular, in rangelands, we highlight the concept of grazing catastrophe, according to which plant abundance could vary discontinuously in response to a change in the stocking rate (May, Reference May1977). In the literature, several spatially implicit differential equation models demonstrated grazing catastrophes (May, Reference May1977; Rietkerk and van de Koppel, Reference Rietkerk and Van de Koppel1997; De Michele et al., Reference De Michele, Accatino, Vezzoli and Scholes2011).

The formalization of rangeland dynamics with multiple steady states naturally led to the concept of resilience. Holling (Reference Holling1973) gave two definitions of resilience. Engineering resilience is the rate at which the system returns to a desirable steady state (inverse of return time) after a disturbance. Ecological resilience is the largest magnitude of disturbance that the system can absorb without leaving the attraction basin of a desired steady state. The concept of ecological resilience implies the existence of attractors with their attraction basins (Gunderson, Reference Gunderson2000; Martin et al., Reference Martin, Deffuant and Calabrese2011) and was evaluated as a distance to bifurcation points (Ludwig et al., Reference Ludwig, Walker and Holling1997), or was linked to the width of the attraction basin of the desired attractor (van Coller, Reference Van Coller1997; Anderies et al., Reference Anderies, Janssen and Walker2002; Janssen et al., Reference Janssen, Anderies and Walker2004).

In the literature, there are mainly two concepts that broadened the concept of resilience and contrasted with its link to the concept of attractor. The first was the concept of non-equilibrium (Ellis and Swift, Reference Ellis and Swift1988): rangeland vegetation dynamics are driven by fire events and inter-annually highly variable rainfall (Ellis and Swift, Reference Ellis and Swift1988; Behnke and Scoones, Reference Behnke and Scoones1993; Illius and O’Connor, Reference Illius and O’Connor1999). Thus, variables describing rangeland dynamics do not reach equilibrium; they rather fluctuate in response to environmental inputs (Westoby et al., Reference Westoby, Walker and Noy-Meir1989; Schröder et al., Reference Schröder, Persson and De Roos2005). The second concept is that rangelands are managed agroecosystems that ranchers may be interested in maintaining in certain conditions, which guarantee sustainable exploitation of resources (Janssen et al., Reference Janssen, Anderies and Walker2004); these conditions may be a region of the state space that does not include a steady state (Martin et al., Reference Martin, Deffuant and Calabrese2011).

To address this situation where the desired states are far from the steady states, Martin (Reference Martin2004) adapted the concept of viability to define resilience. In the framework of viability theory (Aubin, Reference Aubin1991), the system has dynamics whose trajectories can be modified and driven by controls. The challenge is to find the controls that maintain through time the trajectories within a so-called domain of constraints defined by a set of viability constraints. Applied to the management of renewable resources, the viability constraints are generally defined to allow good levels of production in the long term while ensuring the preservation of the resource (see e.g. Tichit et al., Reference Tichit, Hubert, Doyen and Genin2004; Martinet and Doyen, Reference Martinet and Doyen2007; Sabatier et al., Reference Sabatier, Doyen and Tichit2012). The viability kernel is the set of initial states (a subset of the domain of constraints) for which there exists at least a sequence of controls maintaining the trajectory in the domain of constraints. The capture basin (called resilience basin in this context) of the viability kernel is the set of initial states for which there is at least one trajectory leading to the viability kernel. Following Martin (Reference Martin2004) and Martin et al. (Reference Martin, Deffuant and Calabrese2011), out of the resilience basin, the resilience is null, because it is impossible to drive the system to the viability kernel. Within the resilience basin, the resilience is positive and may be measured as the inverse of the cost necessary to bring the system to the constraint set (Martin, Reference Martin2004). According to this view, emphasis is put on the sets of states and dynamics rather than on steady states.

The use of viability theory for defining resilience opens new perspectives for an operational definition of robustness and adaptive capacity. Robustness measures how likely the system, in a certain state with a certain control, is to remain in desirable conditions in an uncertain environment (Anderies et al., Reference Anderies, Janssen and Ostrom2004). Alvarez and Martin (Reference Alvarez and Martin2011) formalized this idea in the case where uncertainty applies directly to the state determination and define robustness of a state as its distance to the boundary of the viability kernel (or resilience basin). In this study, we measure robustness as the probability of a control to keep the system in a viable state. This operational definition is close to the one used by Rougé et al. (Reference Rougé, Mathias and Deffuant2013).

Adaptability is the capacity to adjust management choices in the face of changing conditions (Galoppín, Reference Galoppín2006; Smit and Wandel, Reference Smit and Wandel2006). In this study, we measure adaptive capacity in each state with the proportion of viable management controls, accounting for their robustness. The more robust controls the rancher has, the more he/she can cope with changing conditions.

Model

The model describes the interactions between grass, woody plants and accounts for rainfall, grazing and browsing. It is spatially implicit and refers to 1-ha area (variables are expressed in terms of density – biomass per unit area), it is time-discrete with seasonal time step, alternating between wet and dry seasons (Figure 1). Similarly to a study by Accatino and De Michele (Reference Accatino and De Michele2013), in the wet season, vegetation growth depends on rainfall, and is limited by competition, grazing and browsing. In the dry season, vegetation does not grow but is removed owing to senescence, herbivory and possible fire occurrence. If X is a state variable, X _t ^S indicates the value of the variable in the year t at the end of the season S (which can be WET or DRY). Each year conventionally begins with the wet season, thus the time succession of the variable is: X _t ^WET, X _t ^DRY, X _t+1 ^WET, X _t+1 ^DRY … (Accatino and De Michele, Reference Accatino and De Michele2013).

Figure 1

Scheme of the interactions between the model components. State variables are represented by rectangles, the driver is represented by a circle and the controls are represented by hexagons.

The components of the model are the state variables, the driver and the controls (Figure 1). The state variables are T, the density of woody plants, and G, the density of grass. Both variables are expressed in kg/ha. The driver is the rainfall R in the wet season that is expressed in mm/year. The controls are the harvested biomass and the grazing fraction in the wet and the dry season, and the fire decision I in the dry season. The harvested biomass (wet season: H ^WET, dry season: H ^DRY) is defined as the sum of browsing b and grazing g, that is, the amount of woody plant and grass biomass [kg/ha] removed by herbivores in a season, respectively. The grazing fraction (wet season: γ ^WET, dry season: γ ^DRY) is defined as the ratio γ=g/H and represents the fraction (ranging between 0 and 1) of grass in the biomass harvested by the herd. This way, we implicitly consider that the rancher can adjust the grazing/browsing fraction of the herd by adjusting the composition of the herd (mixed grazing) or by other management techniques (i.e. grazing routes) that influence the feeding behaviour of the herd (Meuret, Reference Meuret1996). The harvested biomass and the grazing fraction determine the amount of biomass grazed g and browsed b in a season. Concerning fire decision, the rancher can decide whether to burn or not at the end of each dry season, and thus fire is represented by a binary variable (I=0, the rancher does not burn, I=1, the rancher decides to burn).

Dynamics in the wet season

In the wet season, vegetation growth depends on the rainfall amount as well as on competition between woody plants and grasses, and competition within and among grass species.

$$\openup 9pt\left\{ {\matrix{ {T_{{\rm t}}^{{{\rm WET}}} } \cr {G_{{\rm t}}^{{{\rm WET}}} } \cr } \matrix{ = \cr = \cr } \matrix{ {T_{{{\rm t - 1}}}^{{{\rm DRY}}} {\plus}r_{0} f_{{{\rm RT}}} \left( {R_{{\rm t}} } \right)\varphi _{{{\rm GT}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)\left( {k_{{\rm T}} {\minus}T_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right){\minus}b_{{\rm t}}^{{{\rm WET}}} } \cr {\varphi _{{{\rm GG}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)\varphi _{{\rm TG}} \left( {T_{{\rm t}}^{{{\rm WET}}} } \right)f_{{{\rm RG}}} \left( {R_{{\rm t}} } \right){\minus}g_{{\rm t}}^{{{\rm WET}}} } \hfill\cr } } \right.$$

Woody plant growth (equation (1a)) is modelled as a logistic curve where r ₀ (dimensionless) is a growth term and k _T (kg/ha) is the carrying capacity of woody plants. The growth term r ₀ is reduced by the dimensionless factors f _RT and φ _GT, both ranging between 0 and 1 and representing, respectively, the rainfall limitation on woody plant growth, and the negative influence of the grass stratum from the previous dry season on woody plant growth and establishment of new woody plant seedlings. These two factors are expressed by:

$$f_{{{\rm RT}}} \left( {R_{{\rm t}} } \right)=R_{{\rm t}} /(a{\plus}R_{{\rm t}} )$$

$$\varphi _{{{\rm GT}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)=1{\minus}\delta _{{{\rm GT}}} G_{{{\rm t - 1}}}^{{{\rm DRY}}} /k_{{\rm G}} $$

In equation (2), parameter a (mm) is a half-saturation constant, that is, it represents the amount of rainfall for which the rainfall limitation factor f _RT is equal to 0.5. In equation (3), δ _GT (dimensionless) is a parameter ranging from 0 to 1, describing the negative influence of a unit of grass biomass on the woody plant population, whereas parameter k _G (kg/ha) is the carrying capacity of grass.

In equation (1b), f _RG represents the maximum potential grass biomass made possible, given the amount of rainfall, and ranges between 0 and the carrying capacity of grass k _G. Factors φ _GG and φ _TG are dimensionless and range between 0 and 1 and represent the negative influence on grass production of the previous dry season and woody plants, respectively. These quantities are determined as follows:

$$f_{{{\rm RG}}} \left( {R_{{\rm t}} } \right)=\min \left\{ {\lambda R_{{\rm t}} ,k_{{\rm G}} } \right\}$$

$$\varphi _{{{\rm GG}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)=1\,{\minus}\,\delta _{{{\rm GG}}} G_{{{\rm t - 1}}}^{{{\rm DRY}}} /k_{{\rm G}} $$

$$\varphi _{{{\rm TG}}} \left( {T_{{\rm t}}^{{{\rm WET}}} } \right)=1\,{\minus}\,\delta _{{{\rm TG}}} T_{{\rm t}}^{{{\rm WET}}} /k_{{\rm T}} $$

In equation (4), we use the same linear relation between grass production and rainfall as the one given in Higgins et al. (Reference Higgins, Bond and Trollope2000), where the coefficient λ is equal to 3.369 kg/(mm ha), and we give as the upper limit the carrying capacity of grass k _G. In equation (5), factor δ _GG (dimensionless) ranges from 0 to 1 and represents the negative influence of a unit of grass biomass from the previous dry season on the growth of new grass. In equation (6), δ _TG (dimensionless) ranges from 0 to 1 and represents the negative influence of a woody plant biomass unit on grass production.

Dynamics in the dry season

In the dry season, both the woody plants and grass densities decrease because of senescence. The dynamics with or without fire are different. If the rancher decides not to burn (I=0), the equations are:

$$\openup 9pt\left\{ {\matrix{ {T_{{\rm t}}^{{{\rm DRY}}} =} \cr {G_{{\rm t}}^{{{\rm DRY}}} =} \cr } } \right.\matrix{ {\left( {1{\minus}\mu _{{\rm T}} } \right)T_{{\rm t}}^{{{\rm WET}}} {\minus}b_{{\rm t}}^{{{\rm DRY}}} } \cr {\left( {1{\minus}\mu _{{\rm G}} } \right)G_{{\rm t}}^{{{\rm WET}}} {\minus}g_{{\rm t}}^{{{\rm DRY}}}} \cr } $$

where parameters μ _T and μ _G (dimensionless) represent, respectively, the fraction of woody plant and grass biomass removed because of senescence within the entire dry season. If the rancher decides to burn (I=1), the equations are:

$$\openup 9pt\left\{ {\matrix{ {T_{{\rm t}}^{{{\rm DRY}}} =} \cr {G_{{\rm t}}^{{{\rm DRY}}} =} \cr } } \right.\matrix{ {\varphi _{{\rm F}} (G_{{\rm t}}^{{{\rm WET}}} )\left( {\left( {1{\minus}\mu _{{\rm T}} } \right)T_{{\rm t}}^{{{\rm WET}}} \,{\minus}\,b_{{\rm t}}^{{{\rm DRY}}} } \right)} \cr 0 \hfill \cr } $$

where the factor φ _F (dimensionless) ranges between 0 and 1, represents the damage of fire on trees and depends on the amount of grass in the ecosystem. This factor is given by the following equation:

$$ \varphi _{{\rm F}} \left( {G_{{\rm t}}^{{{\rm WET}}} } \right)=1{\minus}{\varepsilon}\left( {\left( {1{\minus}\mu _{{\rm G}} } \right)G_{{\rm t}}^{{{\rm WET}}} {\minus}g_{{\rm t}}^{{{\rm DRY}}} } \right)/k_{{\rm G}} $$

We make the simplifying assumption that the rancher burns at the end of the dry season, and thus the severity of fire is determined by the amount of grass at the end of the dry season before burning. Parameter ε (dimensionless) represents the fire damage provoked by a unit of grass biomass (normalized by grass-carrying capacity) on woody plants.

Viability constraints for the model

The constraints refer to the controls and do not impose any a priori conditions on the states the system could visit. In this way, our analysis differs from the analysis by Calabrese et al. (Reference Calabrese, Deffuant and Grimm2011).

A production constraint ensures that the harvested biomass in both seasons is always greater than or equal to a quantity H ^min

$$H^{{{\rm WET}}} \ges H^{{{\rm min}}} \,\,\,{\rm and}\,\,H^{{{\rm DRY}}} \ges H^{{{\rm min}}} $$

Two biological constraints ensure that the needs of the animals are satisfied. In both seasons, the grazing fraction should be between two limits γ ^inf and γ ^sup, reflecting the dietary ability of the herd.

$$\gamma ^{{{\rm inf}}} \les \,\gamma ^{{{\rm WET}}} \les \,\gamma ^{{{\rm sup}}} \,{\rm and}\,\gamma ^{{{\rm inf}}} \les \,\gamma ^{{{\rm DRY}}} \les \,\gamma ^{{{\rm sup}}} $$

The biomass needed by the herd should not exceed the available biomass, that is, in the wet season:

$$\openup 9pt\left\{ {\matrix{ {b_{{\rm t}}^{{{\rm WET}}} \les } \cr {g_{{\rm t}}^{{{\rm WET}}} \les } \cr } } \right.\matrix{ {T_{{{\rm t - 1}}}^{{{\rm DRY}}} {\plus}r_{0} f_{{{\rm RT}}} \left( {R_{{\rm t}} } \right)\varphi _{{{\rm GT}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)\left( {k_{{\rm T}} {\minus}T_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)} \cr {\varphi _{{{\rm GG}}} \left( {G_{{{\rm t - 1}}}^{{{\rm DRY}}} } \right)\varphi _{{{\rm TG}}} \left( {T_{{\rm t}}^{{{\rm WET}}} } \right)f_{{{\rm RG}}} \left( {R_{{\rm t}} } \right)}\hfill \cr } $$

and in the dry season:

$$\openup 9pt\left\{ {\matrix{ {b_{{\rm t}}^{{{\rm DRY}}} } \cr {g_{{\rm t}}^{{{\rm DRY}}} } \cr } } \right.\matrix{ {\les \left( {1{\minus}\mu _{{\rm T}} } \right)T_{{\rm t}}^{{{\rm WET}}} } \cr {\les \left( {1{\minus}\mu _{{\rm G}} } \right)G_{{\rm t}}^{{{\rm WET}}} } \cr } $$

For computing the viability kernels for different levels of rainfall, we discretize the values of the state variables and of the controls. We define a state s as a couple of discretized values of T and G. For the wet season, we define as a control c a couple of discretized values of H and γ, and for the wet season a triplet of discretized values of H and γ and a value of I. The controls of the dry season can be divided into the controls with I=0 (without fire), and the controls with I=1 (with fire). We use the algorithm by De Lara and Doyen (Reference De Lara and Doyen2008) to compute viability kernels. Starting from the last time step and going backwards, keeping the rainfall level constant, the viability kernel quickly becomes periodic. It is different for the wet and the dry season, and for the dry season it is possible to find a viability kernel allowing only controls without fire and a viability kernel allowing only controls with fire.

Robustness and adaptability indices

Wet season

A first index ROB(c,s), accounting for robustness against rainfall stochasticity, is defined for the wet season. We define n different rainfall levels (R ₁, R ₂, …R _n ) each one respectively characterized by a probability of occurrence p(R ₁), p(R ₂), …p(R _n ). For each control–state pair (c,s), the robustness index ROB(c,s) is defined as:

$$ROB({c,s})=\Sigma _{i} Y^{{\rm w}} ({c,s,}R_{i} ) \cdot p(R_{i} )$$

where Y ^W (c,s,R _i ) takes value 1 if the wet season dynamics, starting from state s, applying control c, with the rainfall value R _i (using equation 1a,b) leads to a state that is viable for all rainfall values, and takes value 0 in the opposite case. If no control is viable for the state s, the robustness index ROB(c,s) is null. This index represents the probability that being in the state s and applying the control c, the system stays in a viable state (Rougé et al., Reference Rougé, Mathias and Deffuant2013). This index is an adaptation of the robustness index developed by Alvarez and Martin (Reference Alvarez and Martin2011) to the situation where stochasticity occurs as a perturbation of the dynamics rather than a direct perturbation of the state.

We then define an index of adaptability accounting for the proportion of viable decisions in a context of rainfall stochasticity AIW(s). This index is obtained for each state, by summing ROB(c,s) for all possible controls, and dividing by the number of possible controls.

$$AIW({s})={1 \over {{\rm Card(\it C)}}}\Sigma _{{{c}\in {\it C}}} ROB({c,s})$$

If AIW(s) is 0, it means that in state s there are no available controls leading to a state that is viable, regardless of the rainfall.

Dry season

Two adaptability indices are defined for the dry season: AID _NF, if no fire is applied, and AID _F, if fire is applied. There is no stochasticity during the dry season, and these indices only account for the number of viable decisions. They are defined as follows:

$$AID_{{x}} ({s})={1 \over {{\rm Card}({C})}}\Sigma _{{{c}\in {C}}} Y^{{\rm D}} ({c,s})$$

where x stands for F or NF, Y ^D(c,s) takes value 1 if the dynamics in the dry season is the state s, and applying the control c (using equations 7a,b or 8a,b) leads to a state whose index AIW is positive. Indices AID _NF and AID _F are obtained using equation (16), respectively, considering only controls without fire, and only controls with fire.

For readability, we normalized the robustness indices to compare the value they assume in different scenarios. The values of AIW(s) are normalized by dividing them by the maximum value of AIW(s) obtained in different scenarios. The values of AID _NF(s) and AID _F(s) are normalized by dividing them by the maximum value of AID _NF(s) or AID _F(s) obtained in the different scenarios. In this way, the value of the most robust state is 1, and the other values range from 0 to 1 proportionally.

Scenarios explored

We explore the viability kernels and the robustness indices in four different scenarios: low-rainfall cattle-only, low-rainfall mixed-herd, high-rainfall cattle-only and high-rainfall mixed-herd. We consider the range 200 to 700 mm/year as ‘low rainfall’ and the range 700 to 1200 mm/year as ‘high rainfall’. Within each range, we consider a discretization step of 50 mm/year, and we assume that all the rainfall levels have the same probability of occurrence equal to 1/11. The grazing fraction constraint in the cattle-only scenarios is defined to ensure a high grazing fraction (0.8 to 1.0) to reflect the low browsing ability of cattle. Although the herd is only composed of cattle, we consider that some management practices can drive (Meuret, Reference Meuret1996; Provenza and Balph, Reference Provenza and Balph1988) their feeding behaviour towards a little share of browsed biomass (Meuret et al., Reference Meuret, Débit, Agreil and Osty2006). In mixed-herd scenarios, we consider herds composed of cattle and goats. With mixed herds, the farmer has the ability to adjust the grazing fraction of his herd through both management and herd composition. Goats have a higher browsing ability than cattle but are not pure browsers. The browsing fraction in the diet can reach 60% (Pérez-Barberìa et al., Reference Pérez-Barberìa, Elston, Gordon and Illius2004). We therefore defined the constraints on the grazing fraction as 0.4 to 1.0. The values of the other parameters are kept the same for all scenarios (see Table 1). For the state variables T and G we use a discretization step of 200 kg/ha ranging from 0 to, respectively, k _T and k _G. For the control H, we use a discretization step of 200 kg/ha ranging from 0 to k _T+k _G. For the control g, we use a discretization step of 0.05 ranging from γ _i ^nf to γ ^sup.

Table 1

List of parameters with their values kept constant in all the situations explored

The symbol ‘–’ indicates that the parameter is dimensionless.