Pattern Dynamics of Vegetation Growth With Saturated Water Absorption

Regular pattern is a typical feature of vegetation distribution and thus it is important to study the law of vegetation evolution in the fields of desertification and environment conservation. The saturated water absorption effect between the soil water and vegetation plays an crucial role in the vegetation patterns in semi-arid regions, yet its influence on vegetation dynamics is largely ignored. In this paper, we pose a vegetation-water model with saturated water absorption effect of vegetation. Our results show that the parameter 1/P, which is conversion coefficient of water absorption, has a great impact on pattern formation of vegetation: with the increase of P, the density of vegetation decrease, and meanwhile it can induce the transition of different patterns structures. In addition, we find that the increase of appropriate precipitation can postpone the time on the phase transition of the vegetation pattern. The obtained results systematically reveal the effect of saturated water absorption on vegetation systems which well enrich the findings in vegetation dynamics and thus may provide some new insights for vegetation protection.


INTRODUCTION
In nature, vegetation is very widely distributed in different places all over the world. At the same time, vegetation, as a producer in nature, converts carbon dioxide into carbohydrates through photosynthesis, which ensures the food source for humans and animals and keeps the content of carbon dioxide and oxygen in the environment relatively stable [1,2]. Moreover, the soil and water conservation function of vegetation is also very significant. For example, vegetation can reduce the loss of rainwater on the surface and the erosion of the surface soil, and protect the sloping land. Vegetation stems and leaves release water vapor into the atmosphere by transpiration, so that water vapor emitted into the atmosphere and condensed water alleviates drought. Based on the above functions of vegetation, it is particularly necessary to study vegetation dynamics [3][4][5][6].
In recent years, duo to the impact of the greenhouse effect on human life and climate, vegetation plays an indispensable role in climate regulation [7,8]. As for vegetation, people are always concerned about its growth and distribution. There are many factors affecting vegetation distribution, among which climate, geographical conditions and human factors are the most important. Moreover, different conditions will form different vegetation structure, and inhomogeneous distribution of vegetation is called vegetation pattern [9]. Pattern is a kind of non-uniform macroscopic structure with some regularity in space or time, which is ubiquitous in nature, such as stripes of clouds in the sky, waves on the water, figures on the animals and regular spatial pattern which observed in spatiotemporal systems far from equilibrium states [10,11].
Patterns have been extensively studied and a wide range of patterns are found including vegetation patterns [12], infectious disease patterns [13], and patterns on predator-prey systems [14][15][16]. They are induced by different mechanisms and it is vital to understand these mechanisms. Mathematical modeling has become one of the most useful tools in exploring the mechanisms on vegetation dynamics including pattern formation and ecological functions [17]. There are many studies on vegetation pattern. In 1997, Lefever and Lejeune established a single-variable model, which revealed a resource competition mechanism among vegetation communities, namely promotion at short distance and inhibition at long distance [18]. In 1999, Klausmeier firstly proposed the classical vegetated-water model, explaining the regular stripes on the slopes and irregular mosaics on the ground, and pointed out that nonlinear mechanisms play a major role in determining the spatial structure of plant communities [19]. In 2013, Sun et al. revealed the relationship between precipitation and pattern formation: when rainfall is small, the vegetation will form spot pattern; when precipitation increases, the density of the spot pattern will increase, and vegetation appears as spot-stripes mixed pattern with low density [20]. In 2018, Liu et al. proposed a cross-diffusion vegetation system, in which the phenomenon of spot pattern transition was found [17]. In addition, cross-diffusion increased the vegetation density. In 2017, Zhang et al. proposed a vegetation-soil model and explained that wind can induce the generation of vegetation spot pattern. These models do not take into account that vegetation water absorption is not immoderate [21]. When vegetation water absorption reaches a certain degree, vegetation water absorption rate will decrease, which is called the saturation effect of vegetation water absorption. Yuval revealed that high water absorption and rapid diffusion of water in perennial herbs [22]. Of particular interest, this work showed that the pattern transition between multi-steady states is not necessarily catastrophic, yet it can be gradually phasechanged. Based on the observation data of mathematical model, the cause of fairy circles vegetation patch is explained as intra specific competition and the scale dependent effect of vegetation between animals that capture from plants [23]. There are also some work on the early warning signal of desertification [24][25][26].
Water absorption by vegetation is an important process of vegetation growth. The existed work assumed that water absorption is a linear function of vegetation biomass [6,9,19]. However, many types of vegetation have a saturation effect when absorbing water [27][28][29][30], which is generally not well studied by scientists. In fact, this saturated water absorption may have great influences of the vegetation pattern. In this sense, we will show the effect of saturation on the dynamical behavior of vegetation system.
The paper is organized as follows. In Section 2, we pose a vegetation-water model with saturated water absorption of vegetation and mathematical analysis on the emergence of Turing patterns is presented. In Section 3, we reveal the influences of saturated water absorption of vegetation on the patterns and persistence of vegetation system. In the last section, we give some discussion and conclusion.

MATHEMATICAL ANALYSIS
In this section, we will introduce two-dimensional model to descrbe the interactions of vegetation and water, which is posed by Klausmeier [19]: In the above model, there are seven parameters, which are used to depict vegetation physiological phenomena and the change of water. They are all positive depend on what the parameters mean. The first equation of the model (1) represents the change of water. A represents precipitation and water is reduced by evaporation at rate LW. Vegetation absorbs water at rate RG(W)F(N)N, where G(W) W is saturated water absorption on vegetation, and take F(N) N. The second equation of the model (1) is used to simulate the growth process of vegetation, where J is the conversion rate of vegetation into biomass through water absorption and M is lost through mortality. Water flow downhill at speed V and vegetation dispersal is modeled by a diffusion term with diffusion coefficient D.
In this work, we introduce a model contains two variables with saturated vegetation water absorption. This model is more reasonable compare with that model, in which the saturation of vegetation in absorbing water is not taken into account. It is because that the physiological process by which vegetation absorbs water from the soil and forms vegetation biomass is not inordinate, instead, as water increases, it is absorbed by FIGURE 1 | Saturated water absorption on vegetation. When the water concentration is small, the water absorption of vegetation will increase with the increase of water concentration. However, as the water concentration continues to increase, the water absorption of vegetation tends to be constant.
vegetation to a state of saturation. Therefore, we take F(N) aN 1+bN to model the saturated water absorption of vegetation ( Figure 1).
At the same time, due to the diffusion of water, we will add D 1 Δ W to our system, namely: where the biological meanings and units of the parameters in system (2) can be found in Table 1. Let After the original system (2) is dimensionless, the following system is obtained: The initial conditions and boundary conditions are as follows: where L x and L y give region size in the directions of x and y respectively, n ⃗ is the outward unit normal vector of the boundary zΩ, here we consider the boundary zΩ with no flux, namely, Neumann boundary [31][32][33].
In the absence of diffusion, we consider the following system: It is easy to gain that system (3) has a boundary equilibrium E 0 (S, 0) and two positive equilibriums provided that Under conditions (1), (2), (3), we focus on the stability of three equilibriums E 0 , E 1 , and E 2 . The Jacobian matrix corresponding to equilibrium (W p , N p ) as follows: Then we can gain the linearized system: And characteristic equation is: where i) When we consider E 0 (S, 0), one can obtain and thus it is clear that E 0 (S, 0) is stable.
ii) When we consider E 1 (W 1 , N 1 ), then one can obtain the Jacobian matrix of system (5) at equilibrium E 1 :  (2).

Parameter Units Description
Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 721115 Therefore the necessary and sufficient conditions for the equilibrium E 1 being stable is b 11 > 0 and b 12 > 0. Now combining biological significance of each parameter, S > B holds. Then In the following, we consider the sign of b 11 , iii) Next, we investigate the stability of E 2 (W 2 , N 2 ). Similarly, we note the above equation as: Substituting N 2 for b 21 and b 22 , then we can obtain: Then analyzing the sign of b 21 and b 22 . Because Therefore the system has only one stable positive equilibrium E 1 . From biological perspective, we are interested in studying the stability behavior of E 1 . The Jacobian matrix corresponding to E 1 is as follows: In the absence of diffusion, E 1 is stable, whereas become unstable when diffusion is added, which is called Turing instability.
Nonuniform perturbation near the equilibrium point E 1 : where λ is the growth rate of perturbations in time t, and i is the imaginary unit, k is the wave number, r ⃗ (x, y) is the spatial vector in two dimensional space and c. c. stands for the complex conjugate. Substituting (Eq. 9) into (Eq. 8), we obtain characteristic equation: It is easy to get tr k < 0 for any k due to that tr 0 < 0, while the sign of △k is indeterminate. Hopf bifurcation occurs when Im (λ 0 ) ≠ 0, Re (λ 0 ) 0, that is a 11 + a 22 0, a 11 a 22 − a 12 a 21 > 0, then we obtain critical Hopf bifurcation curve a 11 + a 22 0. Then choosing S as Hopf bifurcation parameter, then Im (λ k ) Re (λ k ) 0 at k k T ≠ 0, that is Δ k T 0. And critical wave number satisfies .
We take S as Turing bifurcation parameter, and its critical value S T satisfies the following equation: Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 721115 In this paper, taking D 30, B 5, we gain the Turing region of system (3). In this region, stationary patterns can be observed ( Figure 2). In addition, we obtain the dispersion relation, and find that the real part of the eigenvalue R(λ) increases as the value of P increases. Moreover, Turing pattern will appear within the appropriate parameter range. Dispersion relation shows that when there is no space, the equilibrium point E 1 is stable. When combined with space, loss of stability occurs in relation to the wave numbers. These curves reveal that saturated water absorption induces the instability of system (3).

MULTIPLE SCALE ANALYSIS FOR TURING PATTERNS
The standard multiple-scale analysis yields the well-known amplitude equations. Close to the onset S S T , the eigenvalues associated to the critical modes are close to zero, and they are slowly varying modes, whereas the off-critical mode relax quickly [15,34]. Consequently, the whole dynamical behaviors can be mainly determined by the dynamics of the active slow modes. The stability and the selection of the different patterns close to onset can be derived from the amplitude equations that govern the dynamics of these active modes. Turing patterns (e.g., hexxagon and stripe patterns) are thus well described by a system of three active resonant pair of modes (k j ,−k j ) (j 1, 2, 3) making angles of 2π 3 and |k j k T |. We obtain the linearized form of model (3) at the equilibrium point E 1 as follows: We note and we will gain: zx zt a 11 x + a 12 y + a 13 xy + a 14 y 2 + a 15 xy 2 + a 16 y 3 + DΔx, zy zt a 21 x + a 22 y + a 23 xy + a 24 y 2 + a 25 xy 2 + a 26 y 3 + Δy.
Close to onset S S T , the solutions of model (5a) (5b) can be expanded as At the same time, the solution of model (12) can be expanded as where U S represents the uniform steady state. A j and the conjugate A j are the amplitudes associated with the modes k j and −k j , respectively. The amplitude equations are described through the equations: Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 721115 where μ (S T − S)/S T is a normalized distance to onset, τ 0 is a typical relaxation time. In the following, we will give the exact expressions of the coefficient τ 0 , h, g 1 and g 2 . Setting X (x, y) T , N (N 1 , N 2 ), model (12) can be converted to the following system: where L a 11 + DΔ a 12 a 21 a 22 + Δ , During the calculation, we just analysis the behavior of the parameter close to onset S S T . With this method, we can expanded S in the following term: where ε is a small parameter. Expanding the variable X and the nonlinear term N according to this small parameter, we have the following results: where h 2 and h 3 are corresponding to the second and the third order of ε in the expansion of the nonlinear term N. At the same time, the linear operator L can be expanded as follows: where Here one can have the expression of a * ij by substituting A T for A in a ij and b ij is easy to be obtained. As for multiple-scale analysis, what really pivotal is that we can separate the dynamic behavior according to different time or spatial scale. We only need to separate the time scale for model (16) (i.e., T 0 t, T 1 εt, T 2 ε 2 t). Each time scale T i can be considered as independent variable. The derivative with respect to time becomes the following form: Since that amplitude A is a variable that changes slowly, the derivative with respect to time z zT0 , which changes fast does not effect on the amplitude A. As a result, we have the following result: By using the Eq. 19, Eq. 20, Eq. 21, Eq. 22, and expanding Eq. 15. according to different orders of ε, we can obtain three equations as follows: The first order of ε: The second order of ε: The third order of ε: As for the first order of ε: as L T is the linear operator of the system close to the onset, (x 1 , y 1 ) T is the linear combination of the eigenvectors that corresponding to the eigenvalue 0. Solving the first order of ε, we can obtain: where |k j | k * T , l a * 11 −a * 22 D 2a * 12 D . W j is the amplitude of the mode e ik j r⃗ when the system is under the first order perturbation.
For the second order of ε, we can obtain: According to the Fredholm solubility condition, the vector function of the right hand of Eq. 25. must be orthogonal with the zero eigenvectors of operator L + c , where L + c is the adjoint operator of L + c . In this system, the zero eigenvectors of operator The orthogonality condition is where F i x and F i y , separately, represent the coefficients corresponding to e ikjr ⃗ in F x and F y . Taking e ik1r ⃗ for instance, we will gain The coefficient in Eq. 30. are obtained by solving the sets of the linear equations about exp (0), exp(ik j r ⃗ ), exp(i2k j r ⃗ ), exp(i(k j − k k )r ⃗ ).
With this method, we have  Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 721115 For the third order, we can gain Using the Fredholm solubility condition again, we can obtain By transformation of W, the other two equations can be obtained and the amplitude A i can be expanded as For the order ε 2 and ε 3 , we can obtain the amplitude equation corresponding to A 1 as follows: where The other two equations we can gain by transforming the subscript of A. Based on Ref. [10], one can calculate the values of μ i (i 1, 2, 3, 4). When the controlled parameter μ increase to the critical point μ 2 0, the stationary state of the system begins to lose stability. If μ 1 < μ < μ 2 , then the system exists a bistable region in the range of the controlled parameter. The emergence of Stripe patterns derives from supercritical bifurcation which are unstable  Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 721115 8 for μ < μ 3 and stable for μ > μ 3 . When the controlled parameter μ exceed μ 4 , there is coexistence of hexagon and stripe pattern.

MAIN RESULTS
In order to verify the above theoretical results, we carry out numerical simulation by taking B 5, D 30, and S 10. This paper focuses on the spatial distribution of vegetation with the change of parameter P.
Currently, the desertification phenomenon is particularly austere, so it is the fundamental way for people to understand the cause of desertification correctly to master the law of vegetation evolution [35,36]. Figure 3 shows the evolution of the spatial pattern of vegetation at 0, 1,000, 2000, 4,000, 10,000 and 20,000 iterations with p 0.15. In (A), the vegetation distribution presents an irregular uniform mixing state, and the vegetation is uniformly distributed in the two-dimensional space. Then in (B), the vegetation can be observed to begin to gather into spot and stripe, but the density is not high. With time going by, in (C), the vegetation density increased. In (D), the spatial distribution of vegetation changes, forming a low-density stripe state. However, in (E), the structure of vegetation changes from low density stripe to higher stripe state gradually. At last, in (F), the density of the mixed pattern increases to form a clearer mixed pattern, and it doesn't change for a long time. The above figures show the change and evolution of vegetation structure over time. It can be concluded from this figure that random distribution can result in mixed patterns.
In Figure 4, we can find that the pattern is self-organizing. When the parameter(P) describing the saturated water absorption of vegetation changes, self-organizing patterns of different states can be obtained. Figure 4 Figure 5, we can conclude that vegetation density get smaller and smaller with P increasing, which is consistent with the fact that 1/P is proportional to the rate at which vegetation absorbs water to generate vegetation. Meanwhile, we can find that the pattern changes from honeycomb pattern to mixed pattern, then changes from mixed pattern to labyrinth pattern, at last pattern changes from labyrinth pattern to spot pattern. As a result, we can conclude that the change of P induces pattern phase transition. Many researchers have proposed that spot pattern is the early warning of desertification [31,37,38], therefore we can get that P is of great significance in indicating desertification. Moreover, from the tendency of pattern phase transition, it can be found that P also has a vital influence on the ecosystem robustness. The larger P is, the more unstable the system tends to be. In general, the bigger vegetation density corresponds to a more robust ecosystem.
It is well known that precipitation plays a very significant role in vegetation growth. However, the relationship between precipitation and saturated water absorption of vegetation is still unclear. To reveal the relationship between them, we take S 10.5, B 5 and D 30 and perform simulations. Then we gain the effect of P on the pattern phase transition and find several types of typical patterns in Figure 6.
On the circumstance of the increase of S, the pattern phase transition is insensitive to P. That is to say, when P changes from 0.1428 to 0.15 with S 10, the pattern structure changes from honeycomb pattern to mixed pattern, however, when S 10.5, the same transition tendency of which pattern structure from honeycomb pattern to mixed pattern will need a change of P from 0.16 to 0.23. Of course, the increase of the value of precipitation also increases the density of vegetation (Figure 7).

CONCLUSION
In this work, we show the effect of saturated water absorption on the vegetation dynamics based on a mathematical model in the form of reaction-diffusion equations. We gain rich pattern structures including spotted, mixed, stripe, honeycomb, and labyrinth patterns. It is revealed that there is a negative correlation between P and vegetation density. That is to say, the vegetation biomass decreases as the increase of P, and saturated water absorption can induces the pattern transition of vegetation structures in two-dimensional space. In addition, we also conclude that appropriate precipitation increase can postpone the pattern phase transition.
In this work, we focused our attention on the influences of parameter change on the dynamical behaviors. The results showed that small change may induce the behavior shift between different dynamical regions [39]. The findings can also be applied in other related fields, such as ecosystems, disease transmission, evolutions and so on.
It needs to point out that climatic factors are important impact factors for vegetation dynamics [40,41]. In this sense, we need to combine these factors including temperature, illumination and wind to the mathematical models. Furthermore, big data analysis is useful to explore the inherent law of vegetation evolution in both space and time [42,43]. These topics will be well addressed in the further study.