Efficient spectral collocation method for nonlinear systems of fractional pantograph delay differential equations

1.   Introduction

Modern climate change is a world problem that is paid general attention to by all of mankind. Global climatic variability not only affects the living environment of human beings, but also affects the world economic development and social progress. Over the past 100 years, the global climate has been facing a dramatic change distinguished by temperature increasing ^[1,2,3,4], which was related to both natural factors and human activities. The global surface temperature between 2011 and 2020 was 1.1($^\circ$C) greater than that between 1850 and 1900. From 1975 to 2014, $CO_2$ concentrations [$CO_2$] increased from 280 ppm to 387ppm ^[5]. Observational data from China indicates that the increase rate of the annual average temperature in China was much faster than the world during the last 50 years, especially the Tibetan Plateau ^[6,7].

The Qinghai-Tibetan Plateau (QTP) has become one of the most typical areas affected by climate change ^[8]. That is because the QTP is the highest altitude region in the middle latitudes of the world, and the high latitudes and high elevations are more susceptible to global warming. In recent decades, the QTP has undergone rapid warming, and the warming rate is almost twice that of the world ^[9]. The precipitation has increased over the decades, which contributes to the rising of water levels and lakes expanding ^[10]. The climate change of the QTP has a great impact on the adjacent areas and serves as an indicator of global climate change ^[11]. The researches on climate change in the (QTP) have achieved abundant results ^{[12,13,14,15]}.

Qinghai Lake, which is the biggest salt lake in China, is located in the northeast edge of the QTP at an altitude of about 3, 200m, with longitude from $99^\circ36'$E to $100^\circ47'$E and latitude from $36^\circ32'$N to $37^\circ15'$N ^[16], as shown in Figure 1. Qinghai Lake is a natural barrier for controlling the eastward spread of desertification in the western region, while it lies in a monsoon transitional zone, and it is a famous tourism resort in China ^[17]. The average depth of the lake is 21 meters and the maximum depth exceeds 29 meters ^[18]. In the past 50 years, the average temperature in Qinghai Lake area increased by 0.319($^\circ$C) every 10 years ^[19]. It is a typical semi-arid region, Its heterogeneous environments are vulnerable to global climate variability and the ecosystem is fragile ^[20].

In recent decades, Qinghai Lake has attracted increased attention ^{[21,22,23,24,25,26]}. Additionally, there have been some studies on the vegetation of Qinghai Lake. Zhang et al. studied vegetation change in the Qinghai Lake watershed by conducting pollen-based vegetation reconstruction at an archaeological site ^[27]. Wang et al. studied the relationship between grassland vegetation and climatic parameters in Qinghai Lake. The result showed that the main reason for the improvement of vegetation cover in the Qinghai Lake basin was the increase of precipitation ^[28]. The plant community characteristics of different sand-forming communities in the largest desert area on the East Coast of Qinghai Lake were studied and the results showed that species diversity of plant community and herb coverage were positively correlated with dune stability ^[29]. Cai et al. studied the effects of human activities and climate change on vegetation in Qinghai Lake basin ^[30]. We can see that most of the current researches on vegetation of the Qinghai Lake are based on observation data and statistical methods, paying little attention to spatial distribution and the growth of vegetation based on pattern dynamics. The evolution law of vegetation pattern can be qualitatively analyzed based on dynamic equation ^{[31,32,33,34,35]}. Therefore, in this paper, we present a reaction-diffusion equation and apply the pattern dynamics theory to reveal characteristics of temporal and spatial distribution of vegetation.

The study aims to address the question as follows: (1) How to establish a suitable vegetation-cilimitic dynamics model? (2) How do different climitic factors affect the growth of vegetation? (3) How will the vegetation pattern transform under different climate scenarios? (4) What are the measures to prevent desertification? In this study, the conditions of steady-state bifurcation are obtained via theoretical analysis. Moreover optimal control theory provides a framework for avoiding desertification of the ecosystem. Finally, we employ numerical simulations to verify the response of the vegetation system to climate change that aim to try to avoid desertification and enhance the robustness of the ecosystem.

2.   The model

2.1.   Model derivation

Water is an essential condition for maintaining the normal physiological function of vegetation. Water resource for vegetation growth mainly comes from precipitation. When rain falls to the ground, some water seeps into the soil and is absorbed by vegetation, then some forms surface runoff. Taken together, Klausmeier established a vegetation-water model in 1999 ^[36]:

where $N$ and $W$ represent the biomass of vegetation and water, respectively. $A$ is the precipitation, the evaporation rate of water is $L$, vegetation takes up water at rate $RWN^2$, $J$ is the rate of conversion of biomass per unit of water consumption, the natural mortality rate of vegetation is $M$, $D_N$ is the diffusion coefficient of vegetation and water flows downhill at speed $V$.

It is worth mentioning that the shading effect of vegetation can reduce the evaporation rate of water. Here, we mainly consider the growth of vegetation on flat ground. At the same time, the carbon gain generated by photosynthesis promotes plant growth and the carbon loss generated by respiration consumes vegetation biomass. Most of the water absorbed by vegetation is lost to the atmosphere in the form of water vapor through transpiration. The major factors affecting these three physiological processes are [$CO_2$] and temperature. Therefore, based on model (2.1) and the above facts, the dynamic model of vegetation and water is established as follows:

where $\rho$ is the reduced evaporation rate of vegetation due to shading. $D_W$ is the diffusion coefficient of water. The amount of vegetation growth $C_g$ due to photosynthesis can be given by the following expression ^[37]:

where $C_a$ is environmental $CO_2$ concentration, $C_1$ is the photosynthetic conversion coefficient of plant biomass and $C_i$ is the available $CO_2$ concentration between canopy cells. The rate of vegetation loss due to respiration $R_{esp}$ can be approximated by a Michaelis $M_{10}$ function ^[37]:

where $B_R$ describes the basic respiration per unit biomass.

$E_r$ stands for transpiration of vegetation, which can describe the difference between saturation specific humidity, and actual specific humidity and the expression for $E_r$ is as follows ^[37]:

where $g_{canopy}$ is for canopy water transfer, which is related to water absorbed by vegetation, and $q$ is the dimensionless specific humidity. Based on the above analysis, let

where $g_{H_2O}$ is the maximum conductivity to $H_2O$ and $CO_2$, respectively and $\gamma$ is the conversion coefficient of diffusivity difference between $CO_2$ and $H_2O$.

In (2.3), specific humidity is defined as follows: $q = \cfrac{\rho_v}{\rho_d}$, where $\rho_v(kgm^{-3})$ and $\rho_d(kgm^{-3})$ represent the density of water vapor and dry air, respectively. According to Dalton's law, there are

where $P$ is atmospheric pressure, $s$ is the pressure of steam, $R_d$ is a constant and $T_a$ is the absolute temperature. We assume that $p$ is large enough and has $q^* = \frac{0.622s^*}{P}$. According to the above analysis, $E_r$ can be obtained:

According to the Clausius-Clapeyron function, the saturated vapor pressure is determined:

Let relative humidity be $R_h = \cfrac{s}{s^*}$ and one has

Based on the above analysis, we obtained a bivariate dynamics model to study the vegetation growth in Qinghai Lake:

with $J = C_a(1-\frac{C_i}{C_a})C_1$ and $q = \frac{0.622}{P}e^*(1-R_h))$. See appendix A for explanations of parameters in the model (2.4).

2.2.   Stability analysis

In this subsection, we shall demonstrate the occurrence of the Turing pattern by stability analysis for system (2.4). The steady states of (2.4) are

where $\Phi = (ARJg_{co_2}+\rho R_{esp})^2-4R\gamma g_{co_2}qR_{esp}^2$, $E_1$ and $E_2$ only exist if $\Phi > 0$. $E_0$ is the bare ground equilibrium.

In what follows, the stability of the steady states will be discussed. We assume the condition $\Phi > 0$ holds so that $E_1$ and $E_2$ are biologically meaningful.

Let

The linearization of (2.4) at $E^*$ is

with

where

Consider the spatially heterogeneous perturbations ^[38,39]:

where $\kappa$ is a wave-number and $\lambda$ reprsents a growth rate of perturbation in $t$. Substituting the above formula into (2.5), the characteristic equation is given:

It follows from (2.6) that the characteristic equation of (2.5) is:

where

In accordance with the above derivation, our result reads as follows.

Theorem 2.1. Suppose that $\Phi > 0$, then the bare-soli steady state $E_0$ is always stable and the positive steady state $E_2$ is unstable.

Proof. The characteristic equation corresponding to the bare-soli steady state $E_0$ is

where

It is easy to see that $\beta_1(\kappa) > 0$ and $\beta_2(\kappa) > 0$ $(\kappa = 0, 1, 2...)$. Therefore, $E_0$ is always stable. Analogously, the characteristic equation is as follows for $E_2$:

where

It is easily seen that $\beta_2(\kappa) < 0$ when $\kappa = 0$, then the positive steady state $E_2$ is unstable. □

In what follows we shall analyze the dynamic behavior of $E_1$. First, the characteristic equation is given:

where

It is easy to check that $\beta_2(0) > 0$ when $\kappa = 0$. Based on the above discussion, the result reads as follows.

Theorem 2.2. Suppose that $\Phi > 0$ holds. If $D_N = D_W = 0$, then $E_2$ is stable for $\beta_1(0) > 0$ and unstable for $\beta_1(0) < 0$.

On the basis of the Turing theory, we can conclude that system (2.4) induces Turing pattern under the two conditions: First, $E_2$ is stable without diffusion; Second, $E_2$ is unstable in the presence of diffusion. As a result, Turing instability occurs only provided that $\beta_1(0) > 0$ and $\beta_1(\kappa) < 0$ for some $\kappa\in\mathbb{N}^+$.

3.   The optimal control problem

Based on the condition for Turing instability deduced in part two, the vegetation patterns with different structures can be presented by numerical experiments. With the increase of precipitation $A$, vegetation patterns change from spot structure to stripe structure (shown in Figure 2), which implies that the robustness of the ecological system is enhanced. Therefore, we can prevent the degradation of the vegetation ecosystem by controlling pattern formations. Here, we aim to get the stripe structure under the case of low precipitation. The optimal control problem provides a powerful tool to realize the aim. We regard the artificial planting rate $r(x, t)$ as a control parameter and rewrite system (2.4) as follows:

The set of admissible controls for $r(x, t)$ is ^[40]:

The objective functional expresses a trade-off between the desired precision and a cost of achieving such precision. Specifically, optimal control aims to lower costs (artificial planting amount) while making the uncontrolled pattern $(N(x, T), W(x, T))$ approach the target pattern $(N_T(x), W_T(x))$. Consider the following optimal control problem:

subject to

where

$J$ is the objective functional, $N_T(x)$ and $W_T(x)$ are the objective patterns and $N(x, t)$ and $W(x, t)$ are state variables. $r(x.t)$ is the control variable and $b_1, b_2, c$ are the constant.

Next, we discuss the expression of an optimal solution.

Constructing Lagrange functional ^{[41,42,43,44,45]}:

Here, the local optimal solution of the optimal control problem is $(N^*, W^*, r^*)$, for any small enough and smooth function $N(x, t)$ with $N(x, t) = 0$. By calculation, one has the directional derivative of the Lagrange functional at $(N^*, W^*, r^*, v_1, v_2)$, which satisfies the following equation:

It follows from the arbitrariness of $N(x, t)$ that $v_1$ satisfies

Analogously, one has

and

Substituting $f_{1, N}, f_{1, W}, f_{2, N}$ and $f_{2, W}$ into Eqs (3.4) and (3.5), the adjoint equation of $(v_1, v_2)$ is:

Note that the allowed control set is a closed convex set. It is clear that the directional derivative of the Lagrange functional at $(N^*, W^*, r^*, v_1, v_2)$ along $r-r^*$ satisfies:

Since $r$ is arbitrary, we substitute $f_{1, r}$ and $f_{2, r}$ into the above inequality, then the following variational inequality can be obtained:

By the variational inequality (3.7), it follows that

where we define the projection $P$ as

4.   The simulation results

4.1.   Response of vegetation pattern to different climatic factors

In this section, we apply the biologically realistic parameters to perform the numerical simulations and research the response of vegetation pattern to climate change, which is based on the climatic data from 1969–2019 of Qinghai Lake. The average values of the climatic factors (precipitation, temperature and [$CO_2$]) are obtained and they are 1.05(mm/d), 0.9879($^\circ$C) and 396(ppm), respectively. The other parameters are fixed: $B_R = 1, Rh = 0.4, g_{co_2} = 10*10^{-3}, M_{10} = 1.6, R = 2.6*10^{-2}, \gamma = 2.5*10^3, C_1 = 12, \frac{C_i}{C_a} = 0.6, \rho = 0.24, D_N = 0.1, D_W = 100$.

Figures 2 and 3, respectively, show the effects of precipitation and [$CO_2$] on vegetation pattern. We can observe that the vegetation patterns change from spot structure to stripe structure as increase of precipitation or [$CO_2$]. The highest density decreases gradually, the lowest density increases gradually and the distribution of vegetation is more uniform. Figure 4 illustrates the variation of vegetation patterns in regard to temperature. With increase of temperature, the transition of vegetation pattern experiences stripe and spot. In contrast with the first two meteorological factors, the highest density increases and the lowest density decreases and the distribution is more uneven, which is not conducive to the robustness of ecosystem. The three climatic factors have different effects on the mean density of vegetation. More precisely, the mean vegetation density is positively associated with rainfall and [$CO_2$], which is in contrast to temperature. This is visualized in Figure 5.

4.2.   Prediction of vegetation pattern under different climate scenarios

In this section, we devote to forecast the future vegetation growth in Qinghai Lake area under three different climate scenarios. The three climate scenarios are simulated data selected from the Coupled Model Intercomparison Project Phase 5 (CMIP5) which has three representative concentration paths (RCP2.6, RCP4.5 and RCP8.5)^[46,47,48]; see Table 1 for an introduction of CMIP5.

We adopt linear regression analysis to statistical temperature, [$CO_2$] and rainfall data under three different scenarios to obtain climate change trends. The results are visualized in Table 2. Temperature and [$CO_2$] increase under the three climate scenarios, and rainfall increases in RCP4.5 and RCP8.5. To predict the future evolution of vegetation pattern in regard to different climate scenarios, Figure 6 shows the variation of vegetation pattern. We can observe that vegetation pattern transitions with the increase of time in RCP2.6: Stripes $\rightarrow$ gap $\rightarrow$ uniform, which indicates that the robustness of the vegetation system is increasing. This is reasonable to infer that the increased robustness of the vegetation system is due to the increase of precipitation and [$CO_2$]. Compared with RCP2.6, the spatial distribution structure changes from stripe to spot in RCP4.5 and RCP8.5, which implies that the ecosystem may undergo degradation, which can finally lead to desertification.

4.3.   The optimal control

In order to increase the resilience of degraded ecosystem and avoid the occurrence of desertification, we offer the optimal conservation strategy$-$human activity (e.g. artificial planting). As illustrated in Figure 2, the transformation from stripe pattern to spot pattern occurs with decreasing precipitation, and spot structure can serve as an early warning signal of the catastrophic shift ^[49]. Moreover, the precipitation has a significant impact on the spatial distribution structure of vegetation. Choose a pattern structure corresponding to $A = 1.3$ as the target pattern, which represents a robust ecological structure. In Figure 7, the snapshots of the optimal control $r(x, t)$ and the associated state variable $N(x, t)$ at rainfall $A = 0.8$ and $1.15$ are presented, respectively. Consequently, artificial planting is an effective way to transform the vegetation pattern into an ideal state.

5.   Conclusions and discussion

For the last few decades, global climate has been facing a great change. The vegetation system exhibits a sensitive response to climate change, especially in arid and semi-arid regions. In this paper, we chose Qinghai Lake as the study area, which is a typical semi-arid area, and studied the response of vegetaion pattern to climate change. We developed a vegetation-water system (2.4) with climatic factors by the zero-flux Neumann conditions and investigated the dynamical behavior. First, we showd the stability of the constant equilibria for (2.4) without diffusion. Moreover, we analyzed the spatiotemporal dynamics at the constant equilibria with diffusion. The conditions for Turing instability of the positive constant equilibrium were obtained in the framework of Turing principle.

Our findings in numerical results revealed that the variation of data for climatic factors had a major impact on vegetation pattern. As precipitation or [$CO_2$] increased, the robustness of the vegetation system was enhanced and the mean density increased. Unlike the first two climatic factors, rising temperatures not only lead to the emergence of spot pattern, but also reduced the mean density of vegetation and accelerated the degradation of the vegetation system. Furthermore, in order to forecast the future vegetation growth in Qinghai Lake, we presented evolution of the vegetation system under different climate scenarios. The results showed that the vegetation system collapses to desertification state in the future in RCP4.5 and RCP8.5 scenarios. This results from the synergy of precipitation, temperature and [$CO_2$]. Compared to RCP4.5 and RCP8.5, due to the dual effects of the increase of precipitation and [$CO_2$], RCP2.6 is a desired climate scenario for Qinghai Lake.

A direction for further study is how to timely avoid desertification. Owing to the method of optimal control, we can induce phase transitions between different vegetation pattern through human activities, such as artificial planting. More precisely, any presented pattern structure (e.g. spot pattern, strip pattern, and gap pattern) can be transformed into a desired pattern (see the Figure 7). As a result, artificial planting contributes to guarding against desertification of vegetation systems, even in low-rainfall regions. The result certifies the effectiveness of the optimal-control method in respect of prevention and control desertification.

This study highlighted the response of vegetation to climate change (precipitation, temperature and [$CO_2$] from modeling point of view. It is necessary to take other climatic factors into account, for instance, light, evaporation and humidity. Based on the seasonality of climate change, another aspect for future work would be to consider the nonautonomous systems; that is, all climatic factors considered in the model are coupled as a function of time. Furthermore, previous work has revealed that spot structure can provide early warning signal for catastrophic shift ^[49], which is from a qualitative point of view. From the quantitative perspectives, we expect to propose a quantifiable indicator for desertification evaluation, and this is also one of the key points for future research.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Fundamental Research Program of Shanxi Province (Grant Nos. 202203021212327 and 202203021211213), Program for the (Reserved) Discipline Leaders of Taiyuan Institute of Technology, Taiyuan Institute of Technology Scientific Research Initial Funding.

Conflict of Interest

The authors declare no conflict of interest.

Appendix A.   Model notation and interpretation