This document was made in the RStudio environment, and expects that you have installed the deSolve, rootSolve, coda and FME R-packages beforehand (which can be done by Install Packages in the Tools menu of RStudio). To start a session in Grind first source the grind.R script:

source("grind.R")
## grind.R (03-05-2024) was sourced

Grind needs to sourced to define its functions. This document uses run(), plane(), newton() and continue().

Saddle-node bifurcations and hysteresis in an arid vegetation

A simple model for the dynamics of water uptake in arid zones (Rietkerk & Van de Koppel, 1997; HilleRisLambers et al., 2001) revealed that the sparse vegetation may collapse when people keep too many cattle in these areas, and that it may very difficult to recover from such a collapse. The main idea of this model is that in arid areas with a low vegetation coverage most of the water that falls on the soil by the sometimes heavy rainfall fails to penetrate into the soil. The surface water is rapidly washed off into rivers and disappears. A vegetation increases the penetration of water into the soil but also increases the evaporation of water from the leaves. Since in arid areas the limiting factor for vegetation growth is the availability of water in the soil, they wrote \[W'=R\Bigl(w_0+\frac{V}{k_2+V}\Bigr) - r_WW -\frac{gVW}{k_1 + W} \quad\hbox{and}\quad V'=i+\frac{cgVW}{k_1+W} - dV -GV\ ,\] where \(R\) is the rainfall (in mm per day), \(V\) is the vegetation biomass (in g per m\(^2\)), and \(W\) is the amount of water in the soil (mm). The model has two saturation constants, \(k_1\) (mm) defines the amount of water at which the vegetation grows at half its maximal rate, and \(k_2\) (g per m\(^2\)) is the vegetation cover at which the penetration of water into the soil is \(R(w_0+1/2)=1.4\) mm per day. The parameter \(G\) defines the grazing rate by the cattle.

Here is the model in Grind with the parameters they estimated

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    Gt <- G + eps*t
    dW <- R*(w0 + V/(k2 + V)) - rW*W -  g*V*W/(k1 + W)
    dV <- i + c*g*V*W/(k1 + W) - d*V - Gt*V
    return(list(c(dW, dV)))  
  }) 
}  
p <- c(eps=0,R=2,c=10,d=0.05,G=0.2,g=0.05,i=0.01,k1=5,k2=5,rW=0.2,w0=0.2)
s <- c(W=5,V=100)

Note that the maximum growth rate of the vegetation is \(cg=0.5\) per day, and that the parameters define a total vegetation loss rate of \(d+G=0.25\) per day. In the Grind model we add the parameter eps that allows us to increases the grazing rate over time. For eps=0 the grazing rate remains fixed at G=0.2 per day.

Draw nullclines and determine the stability of the three steady states

plane(xmin=0,xmax=8,ymin=0.01,ymax=1000,log="y")
low <- newton(c(W=2,V=0.1),plot=TRUE)
##         W         V 
## 2.1930781 0.1025049 
## Stable node, eigenvalues:
## [1] -0.21592831 -0.08212333
mid <- newton(c(W=4,V=2),plot=TRUE)
##        W        V 
## 4.825137 2.247506 
## Saddle point (unstable), eigenvalues:
## [1] -0.24581026  0.03554034
hig <- newton(c(W=5,V=50),plot=TRUE)

##         W         V 
##  4.991785 48.649989 
## Stable spiral point, eigenvalues:
## [1] -0.1610153+0.0184104i -0.1610153-0.0184104i

which reveals that the point in the middle is a saddle point, and that the two other steady states are stable. The lower one corresponds to a poor vegetation coverage, the upper one to a lush vegetation. These two states co-exist for this grazing rate of \(G=0.2\) per day, which would explain patches with a low and a high vegetation cover in the same region with the same amount of cattle.

Two saddle-node bifurcations

Make a bifurcation diagram by following the state in the middle while varying the grazing rate

continue(mid,x="G",y="V",xmin=0.1,xmax=0.3,ymin=0.01,ymax=200,log="y")
## Starting at G = 0.2 with:
##        W        V 
## 4.825137 2.247506

## Turning point point at G = 0.2474141 
## Turning point point at G = 0.1505234
## NULL

Slowly increase the grazing rate

Let the grazing rate increase from Gstart to Gend and start the model at the steady state corresponding to p["G"] <- Gstart. Compute the time required to stop at Gend.

Gstart <- 0.1; Gend <- 0.3; p["G"] <- Gstart
s <- run(state=mid,log="y",timeplot=FALSE) 
p["eps"] <- 0.0001
timespan <- as.numeric(abs((Gend-Gstart)/p["eps"]))

Allow for noise on the rainfall, run the model while slowly increasing the grazing rate, and store the output in the dataframe data:

after <- "parms[\"R\"]<-abs(rnorm(1,2,0.2))"
data <- run(timespan,log="y",table=TRUE,after=after)

Add the numerical solution to the bifurcation diagram

continue(s,x="G",y="V",xmin=0.05,xmax=0.35,ymin=0.01,ymax=200,log="y")
## Starting at G = 0.1 with:
##          W          V 
##   2.141244 126.444920
## Turning point point at G = 0.2474102 
## Turning point point at G = 0.1505313
## NULL
lines(p["G"]+p["eps"]*seq(0,timespan),data$V,lwd=2)

which reveals that the vegetation suddenly collapes close to the upper saddle-node bifurcation (without any apparent warning).

Try to recover by decreasing the grazing rate

Start at the steady state corresponding to p["G"] <- Gend.

Gstart <- 0.3; Gend <- 0.1; p["G"] <- Gstart
s <- run(state=mid,log="y",timeplot=FALSE)  
p["eps"] <- -0.0001
data <- run(timespan,log="y",table=TRUE,after=after)

Add the trajectory to the bifurcation diagram

continue(s,x="G",y="V",xmin=0.05,xmax=0.35,ymin=0.01,ymax=200,log="y")
## Starting at G = 0.3 with:
##          W          V 
## 2.09420089 0.04940699
## Turning point point at G = 0.1505117 
## Turning point point at G = 0.247418
## NULL
lines(p["G"]+p["eps"]*seq(0,timespan),data$V,lwd=2)

which reveals that the vegetation only recovers at the lower saddle-node vegetation. This late switching is called hysteresis.

References

Rietkerk, M., Dekker, S. C., De Ruiter, P. C. & Van de Koppel, J. (2004). Self-Organized Patchiness and Catastrophic Shifts in Ecosystems. Science 305, 1926-1929.

HilleRisLambers, R., Rietkerk, M., Van den Bosch, F., Prins, H. H. T. & De Kroon, H. (2001). Vegetation pattern formation in semi-arid grazing systems. Ecology 82, 50-61.

More information

Check Grind’s homepage for more tutorials