Monitoring Land Degradation through Vegetation Dynamics Mathematical Modeling: Case of Jornada Basin (in the U.S.)

Abstract

Arid ecosystems are known to be sensitive to climate change. The Jornada Basin in the USA, as one representative of arid land, has suffered from land degradation in recent decades. In order to disentangle the climate–vegetation feedback, we analyzed the vegetation dynamics under the effects of climate change via a mathematical model based on the reaction–diffusion mechanism. Using this model, we conducted a sensitive analysis of climate factors and concluded that the ecosystem might experience a catastrophic shift with the climatic deterioration. We considered the non-local interaction term to explain the competition among plants. Additionally, the PLR (power law range) metric was used to quantify the extent of the degradation and to compare the results of the vegetation patterns from the remote sensing data and the simulations. From the results, this model could simulate the trends of land degradation in this area. We found that the land degradation could be mainly attributed to climate changes in recent years. This approach suggests that vegetation patterns can provide hints as to whether the ecosystem is approaching desertification. These results can help with mapping vulnerable arid areas around the world through model simulation and satellite images.

1. Introduction

The increase in atmospheric greenhouse gases in recent decades has significantly impacted global and regional temperatures with concomitant modifications of precipitation patterns [1]. Arid and semiarid ecosystems cover about 40% of terrestrial land and are very sensitive to climate changes. There are growing concerns that these ecosystems have a high risk of land desertification, which may affect approximately 25% of the population around the world [2].

The Jornada Basin (103,373 ha), having a typical closed-basin topography found in many arid areas of the Southwestern United States [3], has suffered from land degradation due to climate change. Research has shown that recent global changes, due to both environmental and human influence, combined with climatic stress, especially a prolonged drought, have impacted the structure and function of the ecosystem in southern New Mexico [4]. Grassland areas with weak resistance to external factors (such as grazing and climate change) tend to shift to alternate states in which desert and shrubs dominate [5]. Given its special research history, the Jornada Basin is one of the natural research sites for the Long Term Ecological Research (LTER) program organized by the National Science Foundation (its location is shown in Figure 1). It is highly possible that many factors, including the rising concentration of atmospheric carbon dioxide [$C{O}_2$], and changes in the seasonal distribution of precipitation, have led to large changes in the structure and function of the Jornada Basin’s ecosystem [6]. Vegetation, as the ecosystem’s engineer, can capture, in its patterns, changes above ground and act as an indicator of the state of the ecosystem [7].

Current vegetation patterns may reflect historical legacies, dynamic patterns of climate variables, resource redistribution across landscapes, and different cross-scale nonlinear interactions [8]. Many researchers agree that dry land ecosystems, such as the Jornada Basin, are not in an equilibrium state. These ecosystems easily experience regime shift driven by the stochasticity and variability in external fluctuations, such as rainfall, temperature, fire, etc. [9,10,11,12]. Inter- and intra-annual rainfall patterns change in arid and semiarid areas due to climate change, resulting in an increased number of extreme events (IPCC AR6) [13,14,15]. Additionally, it is predicted that the mean temperature and atmospheric carbon dioxide will increase, in addition to the changes in the mean annual rainfall, in the Jornada Basin. However, it is difficult to disentangle the state of the current ecosystem. Most studies quantifying the ecosystem state have only used simulation data. Yet, with the increasing availability of remotely sensed data, satellite-driven approaches have been established. Currently, the development of remote sensing and the freely available data have made it convenient and feasible to research the vegetation state of the Jornada Basin [8,16,17]. It has been suggested that vegetation patterns can be applied as an indicator to infer the underlying mechanisms and environmental conditions of an ecosystem [18,19,20,21]. Coupling remote sensing data and methods of pattern dynamics can play a critical role in linking ecosystem structure and function and can advance our understanding of how ecosystems respond to environmental changes.

In this study, a mathematical model was applied to analyze the dynamics of vegetation with the effects of climate change. Then, remote sensing data were used to evaluate the changes in the plants in recent years and to compare the results of the vegetation patterns from the remote sensing data and the simulated model. In Section 2, we described the data and methods used in the analysis. A mathematical model with a non-local term was put forward, which considered the mechanism of vegetation growth in the Jornada Basin to simulate the vegetation patterns. In Section 3, some results of the study were provided. We analyzed the dynamics of the model and obtain the vegetation pattern from the model. Additionally, this study applied the PLR (power law range) to detect the land degradation and compare the results from the observation and simulation. It was found that climate change played a critical role in the land degradation of the Jornada Basin. In Section 4, some of the advantages and limitations of our approach were given. Finally, conclusions were drawn in Section 5.

2. Materials and Methods

2.1. Study Area

The Jornada Basin is located along the northern edge of the Chihuahuan Desert in the southwestern USA. The area includes both the Jornada Experimental Range (JER, 78,413 ha) and the New Mexico State University-owned Chihuahuan Desert Rangeland Research Center (CDRRC, 24,960 ha; Figure 1). The mean maximum monthly surface temperature ranges from 13.7 °C in January to 34.9 °C in July. The mean annual precipitation is about 203 mm [22]. In total, 64% of the rainfall is accumulated in summer. The elevation range of this area is about 1214–1768 m above sea level (excluding the mountain areas). The Jornada Basin’s climate is characterized by cool dry winters and warm summers with monsoon rain. Recent changes in the structure and function of the ecosystem in the Jornada Basin may represent a degradation process that is driven by both environmental and human influence, in combination with climate deterioration, especially prolonged drought [22]. In recent decades, the research on the Jornada Basin has mainly included the following themes: community ecology, land management, animal husbandry, ecosystem sciences, and rangeland improvement [23]. Notably, a significant motivation for most of the research has been the remarkable vegetation change in this area [24].

2.2. Data

The daily temperature and precipitation data were sourced from the NOAA weather station at the Jornada Experimental Range headquarters, southern New Mexico, USA, 1914–2019 (https://www.ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USC00294426/detail (accessed on 6 February 2023)). Daily data have been collected from the Jornada Experimental Range since 1914 for the minimum and maximum air temperature and the daily accumulated precipitation, using standard American climatological service instrumentation and procedures. The included data were transcribed directly from the original data sheets and underwent quality control and assurance procedures. The remote sensing images of the Jornada Basin were from Google Earth. Google Earth images come from various sources, such as the QuickBird (0.6 m/pixel), GeoEye (0.4 m/pixel), and WorldView (0.4–0.5 m/pixel) satellite sensors. In this study, we gathered the images that provided a sufficient resolution to visually identify the vegetation patterns.

2.3. Image Analysis

We used the images from the remote sensing of the Jornada Basin sourced from Google Earth with enough resolution to visually identify the vegetation patches. At each site, we applied a 50 m × 50 m plot representative of the vegetation pattern present in the area [25,26]. We applied the K-mean classification approach implemented in R. This classification method could classify the pixels of images into clusters based on their luminance intensity. We used 30 clusters, ranging from 1 (dark) to 30 (whole image) [27,28]. We classified the images according to the luminance threshold, which could detect all the plant pixels in the image using the graythresh and im2bw functions from Matlab [29]. The graythresh function applies Otsu’s method to identify the threshold, which is used as one approximation for autoclassification of the image [30]. We then verified whether the identified threshold was visually correct.

2.4. Vegetation Patch Size Metric

In order to show the state of the ecosystem, we applied the method of calculating the patch size of the vegetation pattern. Research has shown that vegetation patch size distributions (PSD) fit a power law function (Equation (1)) when the ecosystem is robust [31]. However, they deviate from this distribution when the ecosystem is degraded [31,32].

$$P\left(s\right)=\frac{d-1}{s_{min}}{\left(\frac{s}{s_{min}}\right)}^{-d},$$

(1)

where s is the vegetation patch size. $P\left(s\right)$ represents the frequency of certain sized vegetation patches. The minimum size of a patch is $s_{min}$ (patches that are smaller than this are discarded). Parameter d describes the rate of the frequency decay of the patch size. We used the maximum likelihood method to estimate the parameter d. $P\left(s\right)$ and s are log-transformed. The inverse cumulative distribution (the frequency of patches larger than a certain size as a function of size) of a pure power law would appear as a straight line with a negative slope. The approach allows fitting a power law function to those distributions with a heavy tail, including lognormal and truncated power laws, although only a fraction of the data distribution fit the power law. The range of patch sizes that fit the power law distribution is defined as the power law relative range (PLR). The PLR equation is expressed as:

$$\mathrm{PLR}=1-\frac{lo{g}_{10}\left(s_{min}\right)-lo{g}_{10}\left(s_{small}\right)}{lo{g}_{10}\left(s_{max}\right)-lo{g}_{10}\left(s_{small}\right)}.$$

(2)

In Equation (2), $s_{small}$ represents the smallest sized patch, and $s_{max}$ is the largest patch in the image. The range varies from 0 (no data fit the power law distribution) to 1 (all data fit the power law distribution). We could regard the decrease in the PLR as the degradation of the ecosystem. More details can be found in reference [28].

In Figure 2, the flowchart of this study is depicted. The left part of the flowchart describes the model. The right part describes the observation. The historical precipitation, temperature, and $C{O}_2$ data were applied to the model. The PLR was used in both the vegetation patterns from the model and the observations to disentangle the ecosystem state of the Jornada Basin. Finally, we compared the results from the model and observations. We found that the approach could not only simulated the vegetation structure of the Jornada Basin but also simulated the trend of the ecosystem state.

2.5. Mathematical Model

We applied the climate–vegetation model, mainly based on Rietkerk et al. [33] and Kefi et al. [34]. The vegetation pattern is generally linked to the mechanism by which rainfall infiltrates into soil, in combination with the low annual rainfall climate conditions. It can be modeled as:

$$\left\{\begin{array}{c}\frac{\partial P}{\partial t}=g\frac{W}{W+k_1}P-R_{e}P+D_p\mathrm{\Delta }P,\hfill \\ \frac{\partial W}{\partial t}=R\frac{P+k_2{W}_o}{P+k_2}-c\frac{W}{W+k_1}hP-r_{w}W+D_w\mathrm{\Delta }W.\hfill \end{array}\right.$$

(3)

The model mainly includes two parts: water resources (W) and plant density (P). It is based on the assumption that the infiltration of water is related to the plant density of one area. Furthermore, the plant growth includes the effects of $C{O}_2$ fertilization and photosynthesis, which are explained in Appendix A, as well as in ref. [34]. The dynamics of the water density and vegetation biomass are modeled in Equation (3), where R is the rainfall, the first term in $\frac{\partial W}{\partial t}$ represents the distribution of the water, the second term represents the transpiration, which explains the difference between the saturated and actual specific humidity, and the third term represents the water loss due to evaporation and drainage. The first term in $\frac{\partial P}{\partial t}$ represents the part of water absorbed by the vegetation, and the second term represents the effect of photosynthesis and respiration. $\mathrm{\Delta }$ is the Laplacian operator, which equals $\frac{{\partial }^2}{\partial x^2}$ in one dimension and $\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$ in two dimensions. The operator describes the random diffusion of plant and water. The system considers climate factors, such as precipitation, temperature, and $C{O}_2$ concentration [$C{O}_2$]. Further descriptions of the model formation and parameters are in Appendix A, as well as in reference [34].

It is known that the Jornada Basin is mainly covered by shrubs [3]. Shrubs have the effect of being “islands of fertility”, which can gather nutrients under the canopy [35]. When rainfall occurs, it is first absorbed through the canopy, and the excess gathers at the roots. The interaction of competition between vegetation is enhanced by this effect. Based on this effect, we applied one non-local interaction and integral term to describe this phenomenon, as put forward by Zaytseva et al. [36]:

$$\left\{\begin{array}{c}\frac{\partial P}{\partial t}=g\frac{W}{W+k_1}P-R_{e}P+\gamma P\left(x\right){\int }_{-\infty }^{+\infty }V(x-h)P\left(h\right)dh+D_p\mathrm{\Delta }P,\hfill \\ \frac{\partial W}{\partial t}=R\frac{P+k_2{W}_o}{P+k_2}-c\frac{W}{W+k_1}hP-r_{w}W+D_w\mathrm{\Delta }W,\hfill \end{array}\right.$$

(4)

where ${\int }_{+\infty }^{-\infty }V(x-h)P\left(h\right)dh$ is the non-local interaction and integral term representing the competition among the plants. $V\left(x\right)$ is a kernel function [36,37]:

$$V\left(x\right)=\frac{1}{\sqrt{2\pi }}[\frac{1}{{\delta }_1}{e}^{-\frac{x^2}{2{\delta }_1^2}}-\frac{1}{{\delta }_2}{e}^{-\frac{x^2}{2{\delta }_2^2}}],({\delta }_1<{\delta }_2),$$

(5)

where ${\delta }_1^2$ is the activated scale, and ${\delta }_2^2$ is the inhibited scale. In addition, when position h is close to position x, the interaction is positive, indicating that a short distance can benefit plant growth. However, the interaction is negative when position h is far from position x, meaning the long distance inhibits plant growth. Figure 3a shows the characteristics of the kernel function. From the mathematical analysis, we assume that $V\left(x\right)$ has following feature:

$${\int }_{-\infty }^{+\infty }V\left(x\right)dx=0.$$

(6)

From the biological view, vegetation interaction occurs on a relatively small scale. For simple calculation, we substituted Taylor’s expansion for the integral term at $x_0=x$:

$${\int }_{-\infty }^{+\infty }V(x-h)P\left(h\right)dh={\int }_{-\infty }^{+\infty }V\left(z\right)P(x-z)dz={\int }_{-\infty }^{+\infty }V\left(z\right)[P\left(x\right)-z\frac{\partial P\left(x\right)}{\partial x}+\frac{z^2}{2!}\frac{{\partial }^{2}P\left(x\right)}{\partial x^2}-\frac{z^3}{3!}\frac{{\partial }^{3}P\left(x\right)}{\partial x^3}...]dz.$$

(7)

We defined $V_n=\frac{1}{n!}{\int }_{-\infty }^{+\infty }{z}^{n}V\left(z\right)dz$. We could obtain $V_n$ as an odd function of z. $V_n=0$, when n is an odd number. In addition, $V_0={\int }_{-\infty }^{+\infty }V\left(z\right)dz=0$.

Based on the above analysis, we obtained: $V_0=V_1=V_3=0$, $V_2=\frac{{\delta }_1^2-{\delta }_2^2}{2}$, and $V_4=\frac{{\delta }_1^2-{\delta }_2^2}{8}$.

Consequently, Equation (4) is formulated as:

$$\left\{\begin{array}{c}\frac{\partial P}{\partial t}=g\frac{W}{W+k_1}P-R_{e}P+\gamma P\left(x\right)[V_2\frac{{\partial }^{2}P\left(x\right)}{\partial x^2}+V_4\frac{{\partial }^{4}P\left(x\right)}{\partial x^4}]+D_p\mathrm{\Delta }P,\hfill \\ \frac{\partial W}{\partial t}=R\frac{P+k_2{W}_o}{P+k_2}-c\frac{W}{W+k_1}hP-r_{w}W+D_w\mathrm{\Delta }W.\hfill \end{array}\right.$$

(8)

3. Results

3.1. System Dynamics Implementation

For the sake of analyzing the model dynamics, we obtained the equilibrium points of Equation (4) without a spatial term. We denoted $M=g\frac{W}{W+k_1}-R_{e}P,$ $N=R\frac{P+k_2{w}_o}{P+k_2}-c\frac{W}{W+k_1}P-r_{w}W$, and we let M = N = 0. Then, we obtained three equilibria: $E_0=(0,W_0^{*})=(0,\frac{R{w}_o}{r_w}),E_1=(P_1^{*},W_1^{*})$, and $E_2=(P_2^{*},W_2^{*})$, where $W^{*}=W_1^{*}=W_2^{*}=\frac{R_e{k}_1}{g-R_e}$, $P_1^{*}=\frac{R-r_w{W}^{*}-\frac{c{k}_2{W}^{*}}{W^{*}+k_1}+\sqrt{{(R-r_w{W}^{*}-\frac{c{k}_2{W}^{*}}{W^{*}+k_1})}^2-\frac{4c{W}^{*}}{W^{*}+k_1}(r_w{k}_2{W}^{*}-R{k}_2{W}_o)}}{\frac{2c{W}^{*}}{W^{*}+k_1}}$, and $P_2^{*}=\frac{R-r_w{W}^{*}-\frac{c{k}_2{W}^{*}}{W^{*}+k_1}-\sqrt{{(R-r_w{W}^{*}-\frac{c{k}_2{W}^{*}}{W^{*}+k_1})}^2-\frac{4c{W}^{*}}{W^{*}+k_1}(r_w{k}_2{W}^{*}-R{k}_2{W}_o)}}{\frac{2c{W}^{*}}{W^{*}+k_1}}$.

From the ecological view, we selected the equilibrium that was strictly positive under the condition of ${(R-r_w-\frac{c{W}^{*}{k}_2}{k_1+W^{*}})}^2\ge \frac{4c{W}^{*}}{k_1+W^{*}}(k_2{r}_w{W}^{*}-R{k}_2{W}_o)$. So, $E_1$ was used to study the dynamics of the system, hereafter named $E^{*}=(P_1^{*},W_1^{*})$.

3.2. Sensitivity Analysis of the Climate Factors

We analyzed the equation without the spatial term to conduct the bifurcation analysis, which was as follows:

$$\left\{\begin{array}{c}\frac{dP}{dt}=g\frac{W}{W+k_1}P-R_{e}P,\hfill \\ \frac{dW}{dt}=R\frac{P+k_2{W}_o}{P+k_2}-c\frac{W}{W+k_1}hP-r_{w}W.\hfill \end{array}\right.$$

(9)

In the system, when R continues to decrease, the ecosystem is suffering from climate change and would approach the critical threshold shown in Figure 4. These results are in line with Rietkerk et al. [33]. The solid line represents the stable state of the system. When the rainfall keeps decreasing, the system will transform from a stable state (solid line) to an unstable state (dotted line). From an ecological view, the ecosystem might go through a regime shift with the climate’s deterioration.

From the bifurcation diagrams of system (9), we show the sensitive analysis of temperature and [$C{O}_2$] with the change in rainfall. Compared with the [$C{O}_2$], the system more easily reaches the tipping point of a state change from stable to unstable with the increase in the temperature. This means that, under the background of global warming, a regime shift would take place more easily. However, with the increase in [$C{O}_2$], the system is more robust. This is because of the effects of the $C{O}_2$ fertilization that promotes photosynthesis, with evidence from free-air $C{O}_2$ FACE experiments and satellite observations [38,39,40].

3.3. Linear Stability Analysis with a Spatial Term

In order to disentangle the system dynamics from the spatial term, stability analysis was used to discover whether Turing instability would occur in the system (8) [41]. Specifically, the principle of the analysis was as follows. We began from a spatially homogeneous steady state of plant density and water in the conditions, where this equilibrium was stable with spatially homogeneous perturbation. Then, we added or removed small, but spatially heterogeneous, plant biomass and water. When the perturbations grew spatially, the system could develop into a spatial pattern named the ‘Turing pattern’. The range of parameters in a system that has this pattern formation is known as the Turing instability range [41]. The linear formation of the system (8) near equilibrium is:

$$\left\{\begin{array}{c}\frac{\partial P}{\partial t}=a_{11}P+a_{12}W+\gamma P\left(x\right)[V_2\frac{{\partial }^{2}P\left(x\right)}{\partial x^2}+V_4\frac{{\partial }^{4}P\left(x\right)}{\partial x^4}]+D_p\mathrm{\Delta }P,\hfill \\ \frac{\partial W}{\partial t}=a_{21}P+a_{22}W+D_w\mathrm{\Delta }W,\hfill \end{array}\right.$$

(10)

where

$a_{11}=\frac{W^{*}g}{W^{*}+k_1}-R_e$, $a_{12}=\frac{-P^{*}{W}^{*}g}{{(W^{*}+k_1)}^2}+\frac{P^{*}g}{W^{*}+k_1}$,

$a_{21}=\frac{R}{P^{*}+k_2}-R\frac{P^{*}+w_0{k}_2}{{(P^{*}+k_2)}^2}-\frac{W^{*}c}{W^{*}+k_1}$, $a_{22}=\frac{P^{*}{W}^{*}c}{{(W^{*}+k_1)}^2}-\frac{P^{*}c}{P^{*}+k_1}-r_w$.

Let $\left(\genfrac{}{}{0pt}{}{P}{W}\right)=\left(\genfrac{}{}{0pt}{}{P^{*}}{W^{*}}\right)+\left(\genfrac{}{}{0pt}{}{c_1}{c_2}\right)exp(\lambda t+ik\overrightarrow{r})+c.c+O\left({\epsilon }^2\right)$, where k is the wavenumber, $\lambda$ is temporal growth rate of perturbation, and $i^2=-1$. Then, we obtain the characteristic formula:

$$|J-\lambda E-k^{2}D+k^{4}H|=\left|\begin{array}{ccc}{a}_{11}-(a+D_p)k^2+b{k}^4-\lambda & a_{12}& \\ a_{21}& a_{22}-D_w{k}^2-\lambda \end{array}\right|=0,$$

(11)

where $a=\gamma P^{*}{V}_2$, and $b=\gamma P^{*}{V}_4$.

We can simplify Equation (11) as:

$${\lambda }^2-t{r}_k\lambda +{\mathrm{\Delta }}_k=0,$$

(12)

where

$t{r}_k=a_{11}+a_{12}-(a+D_p+D_w)k^2+b{k}^4$,

${\mathrm{\Delta }}_k=-D_{w}b{k}^6+(D_w+D_{w}b+a{a}_{22})k^4-(a_{11}{D}_w+a_{22}a+a_{22})k^4+a_{11}{a}_{12}-a_{12}{a}_{21}$.

Then, we obtain the root of Equation (12): ${\lambda }_{1,2}=\frac{t{r}_k\pm \sqrt{t{r}_k^2-4{\mathrm{\Delta }}_k}}{2}$.

According to the Turing instability criteria [41], the characteristic equation must at least have one positive eigenvalue or complex eigenvalue with a positive real part, as shown in Figure 5; that is, $Re\left({\lambda }_1\right)>0$ or $Re\left({\lambda }_2\right)>0$.

3.4. Vegetation Pattern Formation

Based on the above analysis, we obtain the vegetation patterns in two dimensions. From the simulated vegetation patterns, we could intuitively see the vegetation change in the Jornada Basin with the climate changes (Figure 6).

We used the climate parameters in the Jornada Basin in 2005 and 2017 to see ecosystem change due to climate change. The simulated model results showed that climate change played a critical role in the degradation. The vegetation pattern showed normal patch sizes in 2005. However, in 2017, the vegetation density and vegetation patch size had decreased. The vegetation area (green area) became less, and the bare soil area (yellow area) occupied the largest area of the simulation. In order to compare the simulated results of our model with actual vegetation states from the remote sensing data, we applied the PLR (power law range) metric to illustrate the degradation in Jornada Basin, as described in the next section.

3.5. Degradation Detection

3.5.1. Detection from the Remote Sensing Data

The vegetation’s patch size distribution (PSD) can be used as an indicator of an ecosystem [18]. If the vegetation patch size of an area fits a power law function, we can say the ecosystem is robust and less likely to change when facing disturbance. In this study, the PLR (power law range) was used to analyze the remotely sensed images.

First, two sites at two time periods (2005.7 and 2019.7) in the Jornada Basin were selected for this research. One site was in the JER, and the other was in the CDRRC. We took a snapshot of the sites from Google Earth. Second, we selected the vegetation from the raw images using the K-mean cluster method. Next, the images were smoothed by a Gaussian filter, so the images could be more easily analyzed. Then, the images were binary processed, where the vegetation areas were marked as 1, and the bare soil areas were marked as 0. Finally, we calculated the relationship between the patch frequency and the patch size and used different function relationships to fit it. The best fitting was used to calculate the PLR. In order to select the best fitting function for the data, we applied the exponential, lognormal, power law, and truncated power law function to fit the data from the images and compared them all. We found that the lognormal tended to fit the data best among these functions (Figure 6, Figure 7 and Figure 8).

From the analysis of the two sites, we found that the ecosystem was undergoing a shift to an unstable state. The PLR of site 1 changed from 47% to 37% over the 15 years. The PLR of site 2 changed from 51% to 36%. Both of the two sites had a decreasing PLR, indicating that the state of the ecosystem was worsening.

3.5.2. Detection from the Model Simulation

This study used the vegetation model based on climate variables to simulate the vegetation patterns. Through this model, we attributed the degradation of the Jornada Basin mainly to climate change, which was consistent with the results from Kidron et al. [42]. This approach used the climate parameters in 2005 and 2019. The results are shown in Figure 9. The simulations had the same results as in the above section. The PLR metric was used on the simulated results to calculate the relationship between the patch frequency and the patch size. The decreasing PLR trend implied our model can simulate the degradation. The degradation was observed in Figure 7 and Figure 8. Our model only simulated the vegetation pattern types. When using the climate condition in 2005, the model showed better conditions. The area was mainly occupied by vegetation (green area), and the PLR is 54%. The vegetation type was dense dot. However, in 2019, the vegetation density decreased, and the PLR decreased to 16%. The vegetation type changed to sparse dot. The PLR calculated from our model had the same trend as the observational images, although the image did not coincided with the satellite images. This was mainly because the degradation was not only caused by climate change but also some other factors such as land use and grazing, etc., which our model did not include. The effect of these factors on ecosystems will be researched in our future work. Our model will be improved to research the degradation in the future.

4. Discussion

Globally, dry land ecosystems are extremely sensitive to climate change. The degradation of dry land ecosystems has been an important research topic. Previous research has provided a simplified framework to connect key climate factors with ecosystem patterns by using basic physiological principles [34,43]. However, most research has only studied vegetation patterns in theory, because of the shortage of observational technology [4,33]. With the development of remote sensing technologies, we now have the ability to monitor the vegetation patterns at an adequate spatial resolution [34]. The Jornada Basin, as a typical research site, has received growing concern, due to its changes in recent decades [24]. This study coupled a mathematical model, based on the reaction–diffusion mechanism, and climate variation to simulate the vegetation growth dynamics. It is a novel method to obtain the state of the ecosystem in the Jornada Basin in recent decades. This model simulated the vegetation patterns, and the results were in line with those we obtained from satellite images. Based on this model, we understand more about arid and semiarid areas under future climate change scenarios.

However, other factors, such as land use and radiation effects, etc., may influence vegetation growth. These factors were not included in our model. We need to improve our model and consider the effects of more factors that contribute to vegetation growth in future work. Additionally, more data are needed for our conclusions to be more robust.

5. Conclusions

This study applied a climate–vegetation model, coupled with climate elements and vegetation growth characteristics, to illustrate the vegetation dynamics of the Jornada Basin in relation to climate change. According to the vegetation characteristics, mainly shrubs affected by the ’island of fertility’ effect, we applied the non-local competition term to illustrate the competition based on the distance. We conducted the sensitivity analysis of the climate factors in a non-spatial model and analyzed the bifurcation diagrams with the change in the rainfall. It was concluded that the system more easily reached the tipping point where the state changed from stable to unstable with increase in temperature. This meant that, with global warming, a regime shift would take place more easily. However, with the increase in the [$C{O}_2$], the system was more robust, due to the effects of the $C{O}_2$ fertilization. $C{O}_2$ fertilization was beneficial for the effects of photosynthesis, as evidenced with evidence from free-air $C{O}_2$ FACE experiments and satellite observations [38,39,40].

In order to obtain the parameter range, stability analysis of the model with the non-local term and the spatial term was conducted. Our model had Turing instability with the spatial term. From the stability analysis, we obtained the Turing instability range for the vegetation pattern. The model used the climate data from 2005 and 2019 and simulated the vegetation pattern. The results showed that the area covered by vegetation was decreasing. This evidenced the negative effects of climate change on this ecosystem. It was concluded that climate change has played a critical role in the land degradation in the Jornada Basin in recent years. In order to quantify the land degradation, we applied the PLR to the satellite images, which indicated the extent of land degradation, We selected two sites in the Jornada Basin. Comparing the results between 2005 and 2019, the PLR decreased. This result coincided with the simulation. Our model results provide an interesting hypothesis for field testing. These results could help in future work to map vulnerable arid areas around the world using aerial photographs or satellite images.