Introduction

In arid and semi-arid rangelands worldwide, livestock utilization gradients around watering points are commonly used to investigate the ecological effects of grazing on the biotic and abiotic environment (e.g. Adler and Hall 2005; Peper et al. 2011; Wesuls et al. 2013). These so-called ‘piospheres’ (Lange 1969) are ideal model systems to understand the interplay between grazing, vegetation, and the abiotic environment and can be used to guide range management with respect to landscape functionality and rangeland productivity (Thrash and Derry 1999). Thus, piospheres should be incorporated into more detailed rangeland simulation models that aim at understanding vegetation dynamics under different land-use and climate scenarios.

A piosphere can be understood as a spatially concentrated form of land degradation. It derives from a radial symmetry of ecological impacts around watering points (or similar animal concentration points like stables or kraals). The livestock density per unit area—and thus effects of grazing, trampling, and feces accumulation—is highest in direct vicinity to the watering point and decreases with distance from it until effects start to vanish in the rangeland matrix (Andrew 1988), which may appear in a few hundred meters or several kilometers distance. A characteristic vegetation zoning develops, which commonly includes a largely vegetation-free ‘sacrifice zone’ at the watering point, a zone of denser woody vegetation with reduced grass cover (‘bush thickening’), and beyond this a transition into the typical vegetation composition and structure with usually increasing occurrences of perennial palatable forage grasses (Tolsma et al. 1987; Jeltsch et al. 1997; James et al. 1999; Thrash and Derry 1999). Along this gradient, the spatial distribution of plant species and plant functional types show specific response patterns in their frequency and cover according to the plants’ niche breadths, i.e. the ability to cope with different intensities of livestock disturbance and related changes in biotic interactions and the abiotic environment (Van Rooyen et al. 1991; Heshmatti et al. 2002; Todd 2006; Wesuls et al. 2013).

Though the degradation of vegetation and soil around watering points may represent a permanent state with a low inherent capacity to regenerate (e.g. Jeltsch et al. 1997; Croft et al. 2007), piospheres are not static in time and space. Their overall shape, spatial extent, and the zone-specific vegetation formation may vary considerably with factors such as age of watering point and land-use history, precipitation pattern and season, grazing system and stocking rate, type of livestock, size of management units (camps or paddocks), or due to environmental gradients (Perkins and Thomas 1993; Jeltsch et al. 1997; James et al. 1999; Derry 2004; Washington-Allen et al. 2004; Smet and Ward 2005). Hence, piospheres can show considerable spatial and temporal variation in both the quantity and quality of forage resources. This can have ecological and economic consequences that are highly relevant for animal production systems (Tainton 1999; Todd 2006). Therefore, models of rangeland vegetation dynamics at local (e.g. camp) to landscape (e.g. farm) scales should preferably include the ecology and management of livestock in relation to watering points (Duraiappah and Perkins 1999). While it is possible to initialize spatial simulation models with field data or remotely sensed maps, the ability to generate a range of plausible scenarios is a powerful means for systematically understanding ecological processes (Pe’er et al. 2013). In this regard, complementary to simulation models of grazing gradients, landscape generators can be effective tools.

Landscape generators tend to be less complex than ecological simulation models concerned with the spatiotemporal dynamics of natural systems. Landscape generators are widely used in ecology to systematically create landscapes ranging from near-natural landscapes to artificial extremes for a better understanding of scale-dependent interactions and processes that shape landscape patterns (Gardner et al. 1987; With and King 1997; van Strien et al. 2016). There are basically two approaches to create such patterns by means of landscape generators: process-based and pattern-based (Pe’er et al. 2013). Process-based landscape generators take into account interactions between the different agents in a landscape and their interaction with the environment both in space and time. Agents are collective or individual objects in the observed system (such as livestock or vegetation) that underlie biological and physical processes. Pattern-based models generate realistic landscapes focusing on spatial characteristics of the landscape irrespective of the underlying processes acting in nature (Pe’er et al. 2013). While there exist many different approaches to model piospheres (e.g. Derry 2004; Adler and Hall 2005; Frank et al. 2012; Wesuls et al. 2013; see also Thrash and Derry (1999) for strengths and weaknesses of some), these focus primarily on the emergence of patterns without providing contextualized landscapes as input for simulation models.

In the present study, we report on the development of a piosphere landscape generator (PioLaG) based on a hybrid modelling approach that combines aspects of both process- and pattern-based modelling. PioLaG was originally designed as a complementary component of a spatially-explicit rangeland simulation model for southern African savannas (www.idessa.org) in order to provide vegetation input maps in a user-defined farm setting with realistic, yet modifiable piosphere formations and landscape representations. To our knowledge there is currently no such landscape generator available. Here we introduce and describe the model underlying PioLaG and illustrate model performance with respect to the formation of vegetation patterns around watering points. Results are discussed with respect to the general functionality of PioLaG, its behavior under varying pre-settings defining the environment and management regime, and how well real-world patterns are resembled. Complementary to the model description, we also outline steps in model validation by means of real piosphere formations identified on aerial images covering portions of our reference savanna system (Online Resource 1).

Methods

The ODD (Overview, Design concepts, Details) protocol standardizes the description of individual-based simulation models (Grimm et al. 2006, 2010). In absence of such a protocol for landscape generators, we follow the ODD protocol as closely as possible.

Overview

Purpose

The overall aim was to generate piosphere patterns under a range of environmental and management conditions. In order to serve vegetation modelling and decision support in rangeland management, generated landscapes should feature the number and spatial arrangement of management units (hereafter referred to as ‘camps’), the number and position of livestock watering points therein, as well as the spatial distribution and abundance of different vegetation components defining the overall rangeland condition. PioLaG was conceptualized in reference to Kalahari thornbush-type savannas typically found in northern semi-arid South Africa and bordering Botswana (Mucina and Rutherford 2006). Here, increased livestock activity and grazing pressure around watering points and kraals, and the concomitant reduction in perennial grass cover and suppression of fire is seen as a primary cause of present bush thickening (encroachment of indigenous species) over extensive areas (Perkins and Thomas 1993; Reed et al. 2015; Harmse et al. 2016). Though the application ability of PioLaG should not be restricted to this region, its parameterization was based on the ecology of this savanna type and typical forms of local land use and range management. The underlying model was kept as simple as possible in operating with dominant plant functional types, rule-based biological, pedological, and meteorological information, as well as local and expert knowledge.

Entities, state variables, and scales

PioLaG is spatially explicit and has no temporal progression. The modelled landscapes are grid-based, with cells 30 m × 30 m in size being the basic entity of the model. This cell size was chosen as a middle ground between a very detailed resolution (e.g. 0.5 m × 0.5 m) which requires extremely high computational power, and on the other hand a spatial scale that is too coarse to depict patterns and dynamics of different vegetation components or plant functional types, that are relevant for management purposes. Moreover, piosphere zones would become more difficult to differentiate (e.g. at 250 m × 250 m). Each cell is assigned several state variables, the most important ones being dominant vegetation type, piosphere zone, and herbivore-use intensity. For simplicity and broad applicability, the model only considers basic plant functional types and ‘bare ground’ (altogether referred to as ‘dominant vegetation type’, or short ‘vegetation type’) (Fig. 1). There are six dominant vegetation types: bare soil, annual grasses, unpalatable and palatable perennial grasses, as well as young and mature woody vegetation. Dominant woody vegetation corresponds to bush thickening and subsumes both trees and shrubs. A piosphere zone is the distinct, roughly circular, area of a given vegetation type. Jointly, the zones form the piosphere. The herbivore-use intensity quantifies the impact of livestock on the rangeland depending on the distance to watering point and the stocking ratio (see below). While cell properties define the conditions at the local scale, several cells can define spatially confined conditions of larger spatial extent, i.e. at the multi-cell scale (e.g. clumps of woody vegetation in a grass-dominated matrix). The extent of a landscape is not pre-determined and should be chosen to encompass at least one piosphere. Here we use 400 × 400 cells. Parameters related to climate, soil, and the overall grazing regime are to be set for the landscape scale.

Fig. 1
figure 1

Causal diagram of the PioLaG model showing global parameters (rounded boxes), local cell state variables (square boxes), and their positive (+) or negative (−) causal relationships. The target output variables are shown in the lower left. Primarily external drivers (right part) influence rangeland decision making (dashed arrow), i.e. they determine the internal land-use drivers (upper part). External and internal drivers jointly determine the system responses, i.e. the vegetation response to herbivory at a cell level (bottom left)

Full size image

Process overview and scheduling

The PioLaG landscape generator considers the complex structure of effects from different factors, state variables, and parameters that interact with each other through different processes (Fig. 1). The temporal sequence of the related computational steps (Fig. 2) is rather secondary, except for initialization and farm setup. In the following, the single steps are briefly described. Information about their implementation and related calculations are detailed in ‘Submodels’ section. The model was implemented in NetLogo 5.3.1., an agent-based programmable modelling environment (Wilensky 1999).

Fig. 2
figure 2

Flowchart depicting the process overview and scheduling within the PioLaG model

Full size image

A simulation begins with the setting of global parameters, especially regarding precipitation and soil texture, and the initialization of cell variables. In addition, farm size, exact positions of watering points, and camp layout (arrangement of camps and fences) are set. Next, stocking ratio and herbivore-use intensity (HUI) are calculated for each camp. Each watering point is assigned the highest HUI of all camps associated with a specific watering point. This is because the highest HUI will be used as the upper boundary for calculating the relative probability of the different vegetation types and with that of the piosphere zone widths (Table 2). Upper bounds for the frequency of ‘woody vegetation’ cells are calculated depending on mean annual precipitation and soil texture. Next, the user-defined initial abundance (cell frequency) of (young/mature) woody vegetation is implemented and considered in the calculation of the final overall abundance of woody vegetation. After these preparatory steps, the actual landscape generation begins. The extent of each piosphere zone is classified based on the herbivore-use intensity values of the cells (Table 2). Subsequently, the different vegetation types are distributed on the landscape randomly, based on piosphere zone-specific probabilities. This constitutes the final map of vegetation types. If rotational grazing is used, however, its positive effect on vegetation regeneration is additionally considered via an increase of cells dominated by perennial grasses. Finally, vegetation map, used parameters, and variables are saved for further analysis and use.

Design concepts

Basic principles

PioLaG is a spatially-explicit, hybrid landscape generator, combining properties of process-based modelling (e.g. the effect of precipitation and soil texture on woody vegetation abundance and increased herbivory effects close to watering points) and pattern-based modelling (e.g. zonation of vegetation and clustering of similar vegetation types) (cf. Pe’er et al. 2013). The basic principles underlying this model are the vegetation response to internal and external drivers. The main internal driver is the distance-dependent herbivore-use intensity around livestock watering points. External drivers include precipitation and soil texture, which determine the general environment and are unaffected by local vegetation properties (Fig. 1). At this, the piosphere formation or zoning is based on the assumption of different plant-specific abilities to cope with a certain level of livestock-induced disturbance, as well as differences in attractiveness to livestock (palatability). Vegetation patterning is further influenced by processes like favorable conditions for grass establishment depending on precipitation and scarcity of woody vegetation.

Emergence and interaction

By being spatially-explicit, PioLaG takes into account the fact that spatial patterns (at landscape scale) emerge from ecological processes that usually take place at other spatial scales (e.g. local interactions among vegetation, livestock, and watering points and interactions with global/external drivers) (Levin 1992). The main emergent properties and key outputs of the model are the herbivore-use intensity for each cell and the spatial distribution of dominant vegetation types. Herbivore-use intensity and spatial vegetation distribution shape the circular vegetation zoning around each watering point. At the landscape scale—intended to cover a whole farm or communal rangeland—specific piospheres emerge that can display patterns ranging from quite homogeneous to highly heterogeneous savanna vegetation. By assigning a dominant vegetation type to each cell, the outcome of competition between different vegetation types like tree-grass interactions is modelled indirectly. PioLaG also considers how the user-defined livestock-, farm infrastructure-, and management settings interact with environmental conditions, available resources, and natural processes.

Stochasticity

Vegetation types are randomly assigned to cells, based on piosphere zone-specific probabilities (Table 2). These probabilities depend on herbivore-use intensity, mean annual precipitation, and soil clay content. Thereafter, the resulting spatial distribution of vegetation types is slightly rearranged to achieve clustering of woody vegetation. In case of a rotational grazing system, the ratio of perennial grasses to other vegetation types is increased depending on the duration (and resting) of grazing and on precipitation.

Observation

A main outcome of PioLaG are landscape maps with information about distribution and abundance of vegetation types (in terms of cell frequencies). Further, grazing pressure is determined both for each single cell and per camp unit. The context-specific vegetation zoning and pattern formation can provide an informative basis about the causal relationships of management decisions and their environmental consequences.

Input data

Various input data are needed to run the model (see Table 1 and Online Resource 2). For the current parameterization, some data were directly extracted or calculated from available maps and databases based on GPS coordinates (i.e. precipitation and soil texture) or set manually by the user. Information and data used for parametrization available in peer-reviewed literature were often directly transferable to the model or could be used to develop and improve formulas, such as extents of piosphere zones, vegetation distributions or the effect of clay content and rainfall on shrub encroachment or bush thickening. Other data (e.g. realistic and common values for the number of camps, camp layout, number and position of watering points) were derived from discussions and expert interviews. We held semi-structured interviews with a total of eleven Kalahari farmers in March 2014, with in-depth discussions of specific management practices depending on area, environmental conditions, type of farming and animals kept, etc. We also discussed model features and management practices in the frame of a series of workshops, one with rangeland experts (July 2016, Kimberley, South Africa; 25 participants) and two with commercial farmers (February 2017, Ganyesa and Askham, South Africa; 29 and 28 participants, respectively).

Table 1 Overview of the input data needed to run a PioLaG simulation with description and source
Full size table

Submodels

Initialization: global parameters

For global parameters, the mean annual precipitation (MAP) and relative clay content (as an indicator of soil texture) are read from ASCII maps. The values for MAP and clay content were taken from the CHIRPS dataset (Funk et al. 2015) and the SOTER-based soil dataset for Southern Africa (Batjes 2004), respectively. The value for the grazing capacity is available in official grazing-capacity maps for the respective area or has to be estimated based on the condition of the grass sward at the location. In our scenarios we used typical grazing-capacity values for the study region (Republic of South Africa 1993).

Initialization: farm setup

The farm setup includes parameters for farm size and camp arrangement. The water requirements and water utilization of livestock is highly variable among types of animals and dependent on many factors such as foraging behavior and physiology (Derry 2004). If not being forced to walk longer distances, cattle can be assumed to forage in an area of 3 km around watering points (Tolsma et al. 1987). In South Africa, farm camps in commercial systems are often smaller in size than, for example, in Australia or compared to rangelands in communal areas of South Africa. Accordingly, grazing gradients are usually much shorter and between 1 and 2 km (Smet and Ward 2005), which corresponds to a typical Kalahari 4-camp grazing system with 800-ha camps (Keyser, pers. comm. 2017). However, camps may be smaller and possess two watering points less than 1 km apart, resulting in overlapping piospheres. In contrast, in unfenced Kalahari savannas, the often observed zone of bush thickening may extent over several kilometers (Reed et al. 2015).

PioLaG offers two basic options for camp arrangement according to what is commonly found in the semi-arid Kalahari region: (a) camps in grids or rows with watering points, often installed at central junction points; (b) pie-shaped camps in a wagon-wheel fashion with a central watering point (Online Resource 3, Figures OR3.1-2). However, the number of watering points is flexible and they can be placed anywhere within a camp.

Calculation: stocking ratio

The stocking rate describes the average area of rangeland made available to each livestock unit and is a key parameter of rangeland management (Trollope et al. 1990; Allen et al. 2011). The stocking rate affects vegetation condition and composition according to animal density and herbivore-use intensity. Hence, it depends on factors like season, precipitation, and herbage production. In this regard, within PioLaG, the stocking ratio sets the stocking rate into relation to the grazing capacity of the rangeland, where grazing capacity is the actual area needed to keep an animal without detrimental effects on the rangeland (Trollope et al. 1990). The stocking ratio is calculated as:

$$stocking\, ratio = \frac{{grazing \,capacity\, \left[ {{\text{ha LSU}}^{ - 1} } \right]}}{{stocking\, rate \,\left[ {{\text{ha LSU}}^{ - 1} } \right]}}$$
(1)

An LSU refers to a large stock unit defined as “an animal with a mass of 450 kg and which gains 0.5 kg per day on forage with a digestible energy percentage of 55%” [Meissner (1982) cited in Trollope et al. (1990)]. To keep the model simple, PioLaG differentiates classes of stocking ratios that are based on discussions with Kalahari farmers and rangeland scientists in 2017: ‘understocked’ (if stocking ratio < 0.9), optimally stocked (≥ 0.9 and < 1; i.e. stocking rate slightly below grazing capacity), ‘overstocked’ (≥ 1 and ≤ 1.5) and ‘severely overstocked’ (> 1.5) (see also Online Resource 2, Table OR2.1).

Calculation: herbivore-use intensity

The herbivore-use intensity (HUI) describes the livestock effect on a rangeland. It considers grazing and browsing effects, as well as trampling and excretions by the animals [Georgiadis (1987) cited in Perkins and Thomas (1993)]. In PioLaG, HUI is used as an auxiliary variable and is calculated using a decreasing logistic function of the cell’s distance to the nearest watering point (dToWP [m]), a stocking ratio classes-dependent stocking factor sF (adjusting the steepness; compare Online Resource 2, Table OR2.1), and the calibration constant c1 = 0.005 m−1:

$$HUI = \frac{1}{{(1 + \left( {0.1 \times e^{{( - 2 + \left( {c_{1} \times sF \times dToWP)} \right)}} } \right)}}$$
(2)

The lower the HUI, the lower the livestock impact on the cell. The adjustment of the parameters of the logistic curve alters the gradient of the slope according to how heavily stocked the farm camps are (Fig. 3). With increasing stocking rate (i.e. smaller stocking factors when grazing capacity stays the same), HUI decreases less steep with increasing distance to watering point, because livestock pressure would be higher and spread out over a larger area. Under an overstocking scenario (sF = 0.8), HUI approaches zero not before 2.5 km from the watering point. In contrast, understocking (sF = 1.5) results in a relatively early drop of HUI towards zero already at a distance of about 1.3 km (Fig. 3). When there is no livestock (stocking rate = 0), HUI is uniformly set to one across the modelled landscape to account for grazing wildlife.

Fig. 3
figure 3

Relationship between normalized herbivore-use intensity (HUI) and distance to watering point for stocking factors sF = 1.5, 1.0, 0.9, and 0.8 assigned to the stocking ratio classes understocked, optimally stocked, overstocked, and severely overstocked, respectively

Full size image

Initialization: mean annual precipitation and clay content

An important driver of woody vegetation:grass ratios in savanna systems is the amount of precipitation, which is crucial for seed production and seedling establishment, especially of woody species (Barnes 2001; Joubert et al. 2013). In semi-arid African savannas, maximum cover of woody vegetation increases linearly with mean annual precipitation (MAP), but seldom reaches its maximum potential due to frequent disturbances and interactions with other factors affecting the soil water regime (Sankaran et al. 2005). We use clay content as an indicator of soil texture as it was shown to be the most important soil property influencing bush thickening dynamics. Fine-textured clay soils have smaller pores and lower rain water infiltration compared to sandy soils, thus diminishing growth of trees and counteracting an increase in density or cover of woody species even under otherwise favorable conditions and especially during high rainfall years (Kgosikoma et al. 2012; Grellier et al. 2014). Accordingly, MAP and soil clay content are parameters used in PioLaG to set upper bounds of possible woody vegetation cell frequency. The effect of MAP on the upper bound of bush frequency (%) (maxBush) is calculated based on the linear regression model from Sankaran et al. (2005) with the constant c2 = 0.14 mm−1a:

$$maxBush = \frac{{\left( {c_{2} \times MAP\left[ {mm\, a^{ - 1} } \right]} \right) - 0.42}}{100}$$
(3)

The effect of soil clay content on vegetation is accounted for by differentiating the classes low (< 0.3%), medium (≥ 0.3% and < 0.5%), and high (≥ 0.5%) clay content. These classes are assigned different factors (compare Online Resource 2, Table OR2.2) to influence the vegetation type probability of a cell (see also below).

Initialization: increased woody vegetation frequency

Tree:grass ratios can vary considerably in time, so that a dense woody layer may characterize initial vegetation conditions of one or several whole camps irrespective of the actual piospheres. Accordingly, PioLaG offers the possibility to specify whether and how many camps are already affected by bush thickening at the expense of a perennial grass layer, and at which severity level (‘low’, ‘medium’, and ‘high’). The factors (Online Resource 2, Table OR2.3) are used in the function that assigns the vegetation types to modify the frequency of ‘woody vegetation’ cells (Table 2).

Table 2 Probability (P) of occurrence of the five vegetation types in PioLaG in each of four piosphere zones depending on herbivore-use intensity (HUI). ‘max_HUI’ is water-point specific, i.e. it is the highest cell-wise HUI in the model across all cells assigned to the same watering point. ‘bushFactor’ is the probability of (young and mature) woody vegetation dominating a cell
Full size table

Assignment of vegetation types and spatial distribution

In PioLaG, each cell is assigned to one of four piosphere zones based on HUI (Table 2). The zones are based on a basic piosphere pattern commonly described for Kalahari savannas (e.g. Perkins and Thomas 1993; Jeltsch et al. 1997; Moleele et al. 2002; Smet and Ward 2005), i.e. a low vegetated ‘sacrifice zone’ up to 100–400 m from the watering point, a ‘bush thickening zone’ with increasing woody vegetation density or cover up to 2000 m (depending on camp size), and the outer zone. In PioLaG, the vegetation gradient in the ‘bush thickening zone’ is accounted for by separating this zone into a grass- and subsequent woody vegetation dominated zone. The ‘grass dominated zone’ close to the watering point is characterized by annual and/or unpalatable species, whereas more disturbance-sensitive and perennial forage grasses dominate where HUI and competition by woody vegetation (i.e. the frequency of ‘woody vegetation’ cells) is low. The outer zone can represent conditions from a least impacted grazing reserve with a balanced tree:grass ratio to overall dominance of woody vegetation depending on the pre-set level of bush thickening in the selected camp.

Each cell in the simulated landscape is assigned a ‘bushFactor’ (Eq. 4), calculated as a multiplicative function of the mean annual precipitation (via ‘maxBush’, compare Eq. 3), clay content (‘soilBushFactor’) and HUI (after natural log transformation). The bushFactor describes the probability of woody vegetation dominating a cell:

$$bushFactor = maxBush \times soilBushFactor \times HUI$$
(4)

The vegetation type is then assigned according to the probabilities shown in Table 2. The probabilistic variation accounts for environmental noise. The cells containing a watering point are always set to ‘bare soil’. In a subsequent step, cells set to woody vegetation dominance are further categorized as dominated by either young or mature woody vegetation. ‘Young woody vegetation’ describes vegetation dominated by seedlings and saplings (< 2 m in height) and ‘mature woody vegetation’ describes vegetation dominated by larger trees (≥ 2 m in height). The ratio of cells dominated by young woody vegetation to cells dominated by mature woody vegetation is pre-set as 4:1 based on field evidence (Angassa and Oba 2010; Dreber et al. 2014; Harmse et al. 2016). Young woody vegetation is modelled in a clustered spatial pattern, because clustering of young woody vegetation is commonly observed, likely due to limited seed dispersal (Hesselbarth et al. 2018). Clustering is achieved through swapping of ‘young woody vegetation’ cells of the outer piosphere zone with randomly chosen neighboring cells until most ‘young woody vegetation’ cells have at least two other cells with the same vegetation type in their Moore neighborhood, i.e. the eight cells adjacent to the central cell (Balzter et al. 1998). If the option of bush-thickened camps is selected, the relative frequency of the vegetation types in the outer piosphere zone within the selected camps is altered by using a correction factor (except for bare ground and palatable perennials). The ‘bushFactor’ (see Table 2) is multiplied with this correction factor, which is set to 1, 0.5 or 0.33 for a high, medium, or low level of bush thickening, respectively. Consequently, with a correction factor of 1, bush thickening in the outer piosphere zone is as high as in the woody vegetation dominated zone.

Effect of grazing system

There are various land tenure types in Kalahari savannas including commonage and private as the most common [for an overview refer to Kong et al. (2014)]. In commonages, the rangelands are communally-owned and managed or used with often little to no camp infrastructure. In contrast, independently managed private farms are usually based on a multi-camp system with the associated grazing management varying according to farming objectives, animal type, climate-related constraints, and personal preferences. Camps are either grazed in a continuous or rotational manner, whereas the utilization pressure differs with the stocking ratio and timing of resting periods for the vegetation to recover (Tainton et al. 1999). In this respect, PioLaG enables the recreation of different landscapes such as large, single-camp rangelands with an open and continuous grazing (as often found in communal areas) or multi-camp systems that allow for a more sophisticated type of grazing management. However, the model uses simplified estimations for the calculation of the effect of the management system because in reality there are many individual and complex schemes that depend on the philosophy, flexibility, and immediate response of land users to current environmental, vegetation, and livestock conditions. PioLaG allows the user to disable the effect of rotational grazing. This functionality enables the exploration of effects of rotational grazing on the spatial distribution of vegetation types, by being able to compare rangelands resulting from conditions that differ only in presence and absence of rotational grazing. Farms with multiple camps are usually under rotational grazing by default.

In PioLaG, rotational grazing is modelled via its benefit for perennial grasses (Tainton 1999). Both precipitation and the resting period are key factors for regeneration of perennial grasses. Thus, a function increases the number of cells dominated by perennial grasses, based on the current number of cells being dominated by perennial grasses (‘numPGC’, [cells]), mean annual precipitation (‘MAP’, [mm a−1]), length of the vegetation resting period (‘RP’, [weeks]), and the constant c3 = 0.2 [a]:

$$Increase \,in \,perennial \,grass \,\,cells = floor\left( {numPGC \times \frac{{\left( {\frac{RP}{{4 \,{\text{weeks}}}}} \right) \times \left( {\frac{MAP}{{100 \,{\text{mm}}}}} \right) \times c_{3} }}{100}} \right)$$
(5)

The increase in the number of perennial grass cells is implemented by converting a corresponding number of randomly selected cells from their current vegetation type to ‘palatable perennial grass’ and ‘unpalatable perennial grass’ cells in a ratio of 60:40. In case of continuous grazing (the whole simulated area being a single camp), the model does not consider any further changes in the number of perennial grass cells.

Output

PioLaG output includes distribution maps of dominant vegetation types, locations of watering points, and camp and farm boundaries. The output is saved in separate ASCII grid files with a spatial resolution of 30 m × 30 m. This universal format is ideal both for further analyses and as an input for rangeland simulation models. Additionally, PioLaG is able to directly output the piosphere pattern for a single watering point, as well as HUI as a function of the distance to the watering point, and saves them as a text or image file.

Model performance

In order to test the performance of PioLaG, we generated landscapes 12 km × 12 km in extent with a single central watering point. This allowed us to analyze the dominance distribution of vegetation types over long distances from the watering point. Results for different environmental and management conditions are presented in the following. The complete parameterization for each scenario can be found in Online Resource 5 and the PioLaG model in NetLogo is available from the authors on request.

Sensitivity analysis

For the sensitivity analysis, we used a revised version of the Morris’s elementary effects screening (Campolongo et al. 2007) to identify the most influential parameters and their interaction effects (Morris 1991; Thiele et al. 2014). As central output variables, we analyzed the relative amounts of bare soil, annual grasses, unpalatable and palatable perennial grasses, as well as young and mature woody vegetation. As input parameters of interest, we varied number of bush-thickened camps, soil type (relative clay content), mean annual precipitation and stocking ratio. In Morris screening, μ is the overall effect (first-order effect) of a parameter (called ‘factor’ in Campolongo et al. 2007) on the respective output. Positive and negative values of μ indicate the direction of the effect. μ* denotes the absolute value of μ. The higher the value of μ*, the stronger the effect. Morris screening also reports the standard deviation of the elementary effects (σ) for each parameter. Elementary effects varying a lot (high σ) indicate that other parameters also have an effect on the output variable, i.e. there is an interaction between parameters but this may also be due to non-linearity (Menberg et al. 2016).

Results and discussion

General piosphere patterns

PioLaG recreated radial grazing gradients that revealed distinct patterns in the dominance distribution of vegetation types including bare ground. The corresponding piospheres were clearly visible in form of zones of differing width and vegetation composition around watering points (for an example see Fig. 4). The zone immediately surrounding the watering point, i.e. the ‘sacrifice zone’, was dominated by cells denoted as ‘bare ground’ with few interspersed ‘annual grass’ and ‘unpalatable perennial grass’ cells. The latter two dominated the next zone together with some ‘young woody vegetation’ cells. This zone was followed by a belt of increased frequency of cells denoted as ‘young woody vegetation’ and ‘mature woody vegetation’, indicating a transitional state towards bush thickening. At a recommended stocking rate (e.g. 13 ha LSU−1 at a grazing capacity of 12 ha LSU−1 in Fig. 4), this zone extended to 2 km from the watering point with woody vegetation reaching highest values within a distance from around 0.8 to 1.5 km. ‘Palatable perennial grass’ cells dominated the outer piosphere, as being characteristic for a grazing reserve (Fig. 4). Overall, these basic spatial characteristics and the vegetation zoning resembled patterns described in studies of Kalahari savannas (for references see above).

Fig. 4
figure 4

a Generated landscape of an example simulation run with an optimally stocked rangeland, divided into four camps (white lines) with one central watering point. b Piosphere pattern graph of the above landscape (averaged over the 4 camps), illustrating the change in relative frequency of vegetation types with distance to the watering point. Parameter settings: MAP = 405 mm/a, stocking rate = 13 ha/LSU, grazing capacity = 12 ha/LSU, stocking ratio = 0.92 (optimally stocked), clay content = 0.03 (further details in Online Resource 5, Table OR5.1)

Full size image

Different herbivore-use intensities

The effect of an increasing stocking rate on the simulated rangelands was observed as the spatial expansion of piosphere zones (Fig. 5; grazing capacity: 12 ha LSU−1). In general, the total area dominated by woody vegetation (except for unpalatable perennial grasses) increased with increased stocking. At under-stocked conditions (20 ha LSU−1), the zone with cells dominated by woody vegetation extended to 1.5 km from the watering point (Fig. 5a). At a recommended stocking rate (13 ha LSU−1) and under severely over-stocked conditions (6 ha LSU−1), this zone extended to a distance of 2 km and 2.5 km, respectively (Fig. 5b, c). Despite these differences, the maximum relative frequency of ‘woody vegetation’ cells (young and mature) was similar across all stocking rates (Fig. 5). The area covered by bare soil generally increased with increasing stocking rate. While bare soil was concentrated less than 100 m from the watering point at under-stocked conditions (Fig. 5a), it extended to a distance of about 200 m under severely over-stocked conditions (Fig. 5c). In the under-stocked scenario, ‘annual grass’ cells were most frequent in the area of 50 m to 500 m from the watering point. With increasing stocking rate, this zone extended to 1000 m with peaks at around 40% of relative frequency. Generally, ‘annual grass’ and ‘woody vegetation’ frequency showed a directly opposed pattern, even though recovery of annual grasses was limited following the decrease of woody vegetation further away (Fig. 5). At a low stocking rate, palatable perennial grasses showed an almost linear increase from the watering point up to 1300 m distance before leveling off at 55% relative frequency (Fig. 5a). With increasing stocking rates, this off-leveling shifted until a distance of 2500 m under over-stocked conditions (Fig. 5b, c). The relative frequency of cells dominated by unpalatable perennial grasses was not much influenced by the stocking rate and fluctuated after 50 m to 100 m from the watering point at a frequency between 30 and 50% (Fig. 5).

Fig. 5
figure 5

Comparison of model outputs with different pre-set stocking rates (SR): a below recommended SR, b recommended SR and c above recommended SR. For all parameter settings refer to Online Resource 5, Tables OR5.2-4

Full size image

This response of the output landscapes to varying herbivore-use intensities concurs with studies from arid and semi-arid savanna systems showing that (1) higher grazing pressure promotes the growth of annual generalist species and reduces the competitiveness of especially palatable perennial grasses (van Rooyen et al. 1994; Fynn and O’Connor 2000; Dreber et al. 2011), and (2) changes in the competitive environment can facilitate the increase of woody species (Skarpe 1990; Weber and Jeltsch 2000; Harmse et al. 2016). Further, the different output landscapes and vegetation responses show that PioLaG is able to account for the fact that shape and spatial extent of piospheres are not static (compare Introduction).

Mean annual precipitation and clay content

With increasing mean annual precipitation (MAP), the piosphere pattern and zoning became more distinct (Fig. 6). Mature and especially young woody vegetation were more frequent in a scenario with high MAP (490 mm a−1) and low relative soil clay content (5%) compared to low MAP (124 mm a−1) and high clay content (27%). In the first scenario, the frequency of ‘woody vegetation’ cells (young and mature) peaked at 1200 m distance from the watering point (59%; Fig. 6a). The frequency of ‘unpalatable perennial grass’ cells was highest between 300 and 600 m and slightly leveled off at 30% to 50% thereafter. Palatable perennial grasses were mostly suppressed (up to 20% frequency) but increased towards dominance beyond 1700 m from the watering point (Fig. 6a).

Fig. 6
figure 6

Piosphere graphs of a scenario with a high mean annual precipitation (490 mm/a) and low clay content (5%) in the soil and b with low mean annual precipitation (124 mm/a) and high clay content (27%). For all parameter settings refer to Online Resource 5, Tables OR5.5-6

Full size image

In the low MAP–high clay scenario, ‘woody vegetation’ frequency reached only a maximum of 36% at a distance of 990 m from the watering point (Fig. 6b). ‘Unpalatable perennial grass’ cells dominated at a distance between 100 m and 1700 m. Beyond 600 m, their frequency varied quite constantly between 30% and 45%. The palatable perennial grasses increased almost linearly along the gradient until a constant frequency around 55% was reached at a distance of 2000 m.

As expected, the opposing clay and precipitation scenarios led to basically unchanged piosphere patterns compared to the default simulation (cf. Figure 4 and Fig. OR4.1). The observed zoning and vegetation responses in relation to precipitation is in line with model results of Jeltsch et al. (1997) from another Kalahari savanna. They also found higher precipitation to result in more distinct piosphere zones, as well as an extended zone of bush thickening (Jeltsch et al. 1997). This can be expected as the establishment of woody vegetation in arid to semi-arid savannas is primarily limited by moisture availability (Sankaran et al. 2005). Our results are also in accordance with studies showing that higher soil clay contents can counteract the positive effect of high precipitation on recruitment and growth of woody vegetation (Kgosikoma et al. 2012; Grellier et al. 2014).

Effect of grazing system

In comparison to the use of a rotational grazing system (Fig. 4), the setting of an overstocked continuous (open) grazing system (Fig. 7) led to an increased number of ‘palatable perennial grass’ and ‘unpalatable perennial grass’ cells. Under continuous grazing, the relative frequency of ‘woody vegetation’ cells was only slightly higher but extended over a larger area than under rotational grazing with recommended stocking rate (Fig. 4b). However, these effects were not as strong as expected. Continuous grazing with high animal numbers usually result in more widespread bare soil and a higher cover and/or density of woody vegetation (Teague et al. 2013). However, since only the dominant vegetation type is represented per cell, low to moderate changes in absolute vegetation may be masked, and thus are not quantified and visible.

Fig. 7
figure 7

Piosphere pattern of a landscape with an open grazing system (no camps). For all parameter settings refer to Online Resource 5, Table OR5.7

Full size image

Bush thickening in outer piosphere

Factors such as mismanagement over longer time periods and under extreme weather conditions (e.g. droughts) can result in permanent states of dominant woody vegetation irrespective of the actual grazing gradient (Pickup et al. 1994; Harmse et al. 2016). Accordingly, pre-defining certain camps as having an overall higher woody vegetation density resulted in the outer piosphere being dominated by ‘woody vegetation’ cells at the expense of ‘perennial grass’ cells (Fig. 8). Also, ‘annual grass’ cells could be observed in greater numbers. This is reasonable because these grasses are often better adapted to intense disturbances and may cope better with increased bush cover than perennial grasses (Dreber et al. 2011; Harmse et al. 2016).

Fig. 8
figure 8

Simulated landscape with four watering points and 16 camps, of which six were pre-set as severely bush thickened (lower part of figure). For all parameter settings refer to Online Resource 5, Table OR5.8

Full size image

Sensitivity analysis

In the sensitivity analysis, most parameters behaved as expected, even though not necessarily exactly as in reality because of the abstract, non-temporal approach of PioLaG. The parameters mean annual precipitation (MAP) and number of bush-thickened camps had the largest effect on the frequency of ‘palatable perennial grass’ and ‘woody vegetation’ cells. Overall, the input parameters showed a lot of interaction which is expected for such a complex system.

First-order effects

The relative amount of bare ground in the landscape was only affected by stocking ratio (Fig. 9), and only weakly so, likely because bare ground appeared as the dominant vegetation type only within the sacrifice zone. Regarding annual grasses, the number of bush-thickened camps had the largest elementary effect (Fig. 9). This reflects the association of ‘annual grass’ and ‘woody vegetation’ cells in the outer piosphere zone at pre-set higher levels of bush thickening, as described above. In contrast, MAP had a slightly negative effect on annual grasses, which may result from its simultaneous positive effect on woody vegetation and the resulting competition (Fig. 9). Likewise, palatable perennial grasses were negatively affected by MAP but also by the number of bush-thickened camps, as both increase the woody vegetation at the expense of this vegetation type. Clay content counteracted increases in woody vegetation [as described by Kgosikoma and Mogotsi (2013) and Grellier et al. (2014)] and consequently had a positive effect on (palatable and unpalatable) perennial grasses. In addition to competition also disturbance played a role: the higher the overall grazing pressure (stocking rate), the fewer ‘palatable perennial grass’ cells occurred (Fig. 9). For unpalatable perennial grasses, similar effects of the input parameters were observed with the difference that the initially set number of bush-thickened camps had a weak effect, and direction of this effect was equivocal (Fig. 9). This is because the number of unpalatable perennial grass cells is not directly set by the function that assigns the vegetation types (see Table 2). Its value rather depends on the frequency of the other vegetation types present in the outer piosphere zone. The parameter effects on young and mature woody vegetation reflected the positive relationship with the number of bush-thickened camps and beneficial conditions for woody recruitment and establishment created by an increase of MAP (e.g. Sankaran et al. 2005) and increased grazing pressure (stocking ratio) on competitive strong grasses (e.g. Skarpe 1990), especially in non-clayey soils (see above; Fig. 9).

Fig. 9
figure 9

Results of the sensitivity analysis illustrating the absolute mean values (μ*) on the x-axis and real mean values (μ) on the y-axis of the elementary effects which show the magnitude and direction of the first-order effect of the four input parameters

Full size image

Interaction effects

There were no parameters important only because of their first-order effect (σ/μ* ≤ 0.1), i.e. no (almost) linear responses with respect to the parameters. Hence, the effects of the analyzed parameters on the output variables were interlinked. For the relative amount of bare ground, the stocking ratio showed an almost monotonic behavior and a small interaction effect with the other parameters (Fig. 10). In case of annual grasses, both MAP and number of bush-thickened camps displayed a monotonic behavior (σ/μ* ratio between 0.1 and 0.5). While MAP exhibited only very little interaction, the number of bush- thickened camps interacted more strongly with the other parameters (Fig. 10). For unpalatable perennial grasses, most parameters showed a σ/μ* ratio > 1, indicating large interaction effects and/or a non-linear behavior. MAP behaved differently in having an almost monotonic effect (Fig. 10). For palatable perennial grasses, the number of bush-thickened camps had a stronger overall effect than MAP. However, both parameters had about the same σ, i.e. they experienced the same interaction with other parameters. When set in relation to the total effect on the output variable, MAP had a larger relative interaction effect compared to the number of bush-thickened camps (larger direct effect with the same sigma) (Fig. 10). All parameters for young and mature woody vegetation showed an almost monotonous behavior (σ/μ* ratio between 0.5 and 1) and an increasing absolute interaction effect in the order of stocking ratio, relative amount of clay, number of bush-thickened camps and MAP. However, when set in relation to the total effect on the output variable, the relative interaction effect was similar among the parameters (Fig. 10).

Fig. 10
figure 10

Absolute mean values (μ*) on the x-axis and the standard deviation of the elementary effects (σ) for each parameter on the y-axis show the interaction of the input factors i.e. the second-order effects. The auxiliary lines show the ratios of σ/μ* = 1, 0.5 and 0.1

Full size image

The parameters number of bush-thickened camps, MAP, and clay content showed strong interaction in affecting the output variables palatable perennial grasses and young and mature woody vegetation. This is because the availability of soil moisture plays a major role for bush thickening and is highly dependent on the clay content in the soil and on MAP (see also first-order effects). The input parameter number of bush-thickened camps strongly influences the vegetation composition through its effect on the ‘bushFactor’ (Eq. 4), whose product is used to calculate the frequency of the different vegetation types.

Similarity to observed piospheres

To verify how realistic the output of the landscape generator PioLaG is, we compared piosphere patterns created with PioLaG with patterns identified on aerial images from our reference savanna system in the southern Kalahari. This revealed a good match of the basic patterns (for detailed results and discussion see Online Resource 1). The general PioLaG outputs were also comparable to piosphere patterns described in scientific literature for similar savanna systems as outlined below. However, it should be noted that the appearance of real piospheres can vary greatly in a local area, which, apart from patchy rainfall, may primarily relate to differences in land management and grazing history. Accordingly, the basic vegetation zoning described for piospheres in the savanna literature varies to some extent.

Peaks of increased density or cover of trees and shrubs along grazing gradients are reported to occur at quite different distances from watering points in Kalahari savanna rangelands: 20–200 m (Tolsma et al. 1987), 50–800 m (Perkins and Thomas 1993), 150–250 m (Moleele et al. 2002) and 200–400 m (Smet and Ward 2005). At the recommended stocking rate for the Kalahari reference system, PioLaG generated a frequency-peak of woody vegetation at a distance around 700–1200 m from the watering point (in a scenario with 405 mm MAP), which corresponds to findings of Tobler et al. (2003) from a more humid savanna system. However, other scenario settings would result in different patterns with dominant woody vegetation also closer to the watering point. This shows that woody vegetation-affecting submodels in PioLaG deliver patterns that are absolutely within the spatial range of actually observed zones with increased bush cover or density. The hump-shaped response of the relative woody frequencies in PioLaG was also comparable to simulation outputs of a vegetation model for another Kalahari savanna by Jeltsch et al. (1997). Accordingly, the simulations of the woody vegetation can be considered solid.

Except for some disturbance-tolerant species, herbaceous vegetation often shows a hump-shaped distribution or decreasing trend along grazing gradients towards watering points, especially many perennial grasses decrease under increasing herbivore use intensity (Van Rooyen et al. 1991; Todd 2006; Wesuls et al. 2013). These response patterns were similarly evident in the PioLaG landscapes. The additional distinction between palatable and unpalatable perennial grasses also revealed that the latter were generally more abundant close to the watering point. This resembles preferential grazing patterns and concurs with descriptions by Perkins and Thomas (1993), Smet and Ward (2005), and others (Thrash and Derry 1999). Clearly it was not possible to make such a distinction of the herbaceous layer depicted on the aerial images, but overall herbaceous vegetation increased with distance from the watering point before leveling off (Online Resource 1 and Fig. OR1.1). Therefore, as for the woody vegetation, the output for simulated herbaceous vegetation is reliable in mimicking observed patterns.

It is not trivial to compare the zoning and composition of vegetation at specific distances from watering points between PioLaG outputs and observed patterns. The causes leading to the specific patterns described in the cited studies and depicted on aerial images are complex and highly site-specific and context-dependent. Information about the factors determining shape and extent of the piosphere are often not available in detail or completely unknown, as in the case of the random sample of piospheres from aerial images. Further, PioLaG simulates a single dominant vegetation type for a 30 m × 30 m cell, which certainly leads to scaling issues since vegetation types (including bare ground) with low to moderate frequency are not output by PioLaG. It is therefore not possible to recreate piospheres exactly matching a case from the real-world. However, thanks to its flexibility, PioLaG allows to simulate diverse scenarios representing a solid approximation of what can be found in savanna rangeland systems.

Limitations of the model

Real landscapes are often structurally diverse and show much environmental heterogeneity caused by a multitude of partly interacting system-internal and external drivers across spatial scales. Moreover, extreme climatic events and management decisions of the past may still be reflected in current vegetation patterns, both at the level of species and communities. This complexity cannot be accounted for by a landscape generator like PioLaG and needs to be considered when interpreting generated simplified landscapes in comparison to real landscapes. Our chosen modelling approach, which integrates scientific expert knowledge and perspectives from practitioners, does not fail to produce useful landscape representations; on the contrary. However, there is room for improvement and PioLaG offers the flexibility to add parameters and to adapt processes and routines as needed. A model is only as good as the input data that is used to parameterize and run a simulation. The science of savanna vegetation is a recursive process and with the help of models like PioLaG data gaps become evident and can guide researchers in data collection as required. New data can then be used to refine models and advance outputs toward even more realistic savanna vegetation patterns at local to landscape scale.

Conclusions

PioLaG proved to be a multifunctional landscape generator for southern African savanna rangelands. The created vegetation patterns and the influence of cause-and-effect relationships were in line with the findings of other modelling studies and field observations. PioLaG is optimized for climates with a mean annual precipitation between 100 and 650 mm. However, it can be altered to various environmental conditions and landscape properties and considers the consequences of different management settings. Thus, it can be used for or adapted to a wide range of grazing systems in different biotic and abiotic contexts, such as Australian or South American savannas but also other arid to mesic grassland ecosystems with an increasing or invading tree or shrub component. The simple output feature allows to create various ASCII grids of vegetation type, farm and camp borders, as well as watering point positions showing the effects of environmental change and different management decisions on the spatial vegetation distribution. Such output can be used for environmental education and training, understanding of complex interactions and insights into cause-and-effect relationships. Another use of the ASCII maps is as input data for rangeland models that usually need qualitatively precise, yet user-defined and context-specific input data about initial vegetation conditions. This is especially beneficial for scientists since field data are often limited and repeated simulations and analyses require model runs with many equivalent but not identical realistic landscapes. Overall, these features make PioLaG a highly flexible landscape generator for different ecosystems and regions and a powerful tool for many research questions.