Nonsmooth Dynamic Behaviors Inherited from an Ecohydrological Model: Mutation, Bifurcation, and Chaos

Abstract

The existence of nontrivial dynamic behaviors in a hydrological system is intensively discussed in the literature. However, most of the work has been done from the nonlinear data analysis perspective, with only a few exceptions, due to the mathematical difficulties for theoretical analysis. In this study, a simple but comprehensive enough ecohydrological model with the pulsed atmospheric forcing was developed from the process analysis perspective. The model was then utilized to analyze the non-trivial dynamic behaviors in a coupled ecohydrological system qualitatively and numerically. Our results confirm the existence of multiple stationary states discussed by many researchers. Furthermore, parameter bifurcation was studied and the phenomenon of mutation is found to be rather common. Also, the chaotic characteristic of the system state is obtained under some specific parameters. Parts of these behaviors were seldom reported through the deterministic dynamic analysis done previously.

1. Introduction

As stated by Peterson et al. [1], the dynamic behavior of a hydrological system is always assumed to be primary in traditional hydrology; that is, at the cessation of a transient hydrological disturbance of any magnitude, the system will return to the same stable state and thus show an infinite resilience. This hypothesis is challenged by the dynamic theory, especially the nonlinear theory raised in meteorology [2]. In fact, quite a few studies have been devoted to the “nontrivial” dynamic characters in hydrological systems in the past several decades [3–5]. For example, chaos has been intensively studied in many hydrology-related processes including rainfall, river flow, rainfall-runoff, lake volume, temperature, pressure, and wind velocity (see the summary by Sivakumar [6]). These studies are, however, principally based on the nonlinear data analysis method alone, which can obtain the statistical features from the existing data but cannot indicate the causal relationship between physical processes and dynamic phenomena.

Another point of view, as this study mainly concerns, uses process-based model incorporating intertwined interactions and feedbacks to explore the nontrivial dynamic behaviors. For example, considering the coupled system including the upper soil layer and the atmosphere planetary boundary layer, D'Andrea et al. [7] found the bimodal distribution of soil moisture resulting from the existence of multiple equilibria in the continental water balance. Ridolfi et al. [8] developed a conceptual framework to study the interaction and positive feedback in wetland forests and riparian ecosystems, which shows that the feedbacks may lead to the emergence of multiple stable states. Based on the field observations from savanna ecosystems in Kalahari, D'Odorico et al. [9] found that the difference between the soil moisture in the canopy and intercanopy spaces is related to the aridity level, and the system may exhibit two stable states corresponding to the conditions with and without tree canopy cover. With positive feedbacks between vegetation and soil moisture in arid ecosystems, Borgogno et al. [10] found that random rainfall fluctuation may turn the bistable deterministic system into a stochastic system with only one statistically stable state in an analytical model. In order to challenge the potentially significant assumption of the single stable state in the ecohydrological models, Peterson et al. [1] developed a numerical distributed ecohydrological model and found that the multiple stable states may arise from the reduction of leaf area index as the saline water table approaches the surface. Also, they employed limit cycle continuation to locate the unstable stationary state in the phase space. Runyan and D'Odorico [11] developed a modeling framework relating vegetation-soil salinity feedbacks to the emergence of multiple stable states in the underlying dynamics and the results show that the presence of a strong feedback can result in bistable dynamics for a wide range of environmental conditions.

The above process-based work demonstrates the existence of multiple stationary states in some specific ecohydrological systems, which depends on no particular data or conditions. But more sophisticated dynamic behaviors, such as mutation and chaos, have not been explored yet. Nevertheless, Rodriguez-Iturbe and the collaborators [12, 13] conducted a study with the similar purpose but for the hydrometeorological system operating at the climate scales (i.e., continental spatial scale and annual temporal scale). They analyzed the soil moisture balance equation and demonstrated that soil moisture states could present abundant nontrivial dynamic behaviors including bimodal states, fixed points, limit cycles, and chaos due to the feedback between soil moisture and precipitation (i.e., local recycling of moisture). For the ecohydrological system running at much smaller spatial and temporal scales interested in this study, the dominating feedback is not the local recycle of moisture, but the vegetation growth fed by soil moisture. The evolution of ecohydrological system has recently drawn a lot of attentions and its dynamic behavior deserves further study, which is the purpose of this study.

The remainder of this paper is organized as follows. In Section 2, a simple but comprehensive enough ecohydrological model is established under pulsed atmospheric forcing. Dynamic analysis follows in Section 3, and the so-called nontrivial dynamic behaviors are explored by the way of numerical simulation. The paper is concluded in Section 4.

2. Framework of the Ecohydrological Model

2.1. Functional Expression of Nonsmooth Rainfall Process

As mentioned above, Rodriguez-Iturbe et al. conducted the notable work in the dynamic analysis of hydrometeorologic process. However, the annual rainfall amount in their analysis is regarded as a simple function (in the mathematical sense) of advective precipitation, local evaporation, and local soil moisture, and hence the smaller scale rainfall variability and nonsmooth feature (see below) are ignored. Such description of rainfall process is inapplicable to ecohydrological system running at a smaller spatiotemporal scale. For the detail study of hydrological processes at the finer scale, the Poisson process with white noise is often utilized to describe the stochastic properties of rainfall [14, 15]. However, such stochastic manner will go against the qualitative analysis of the dynamic process due to the mathematical complexity.

A deterministic functional expression of rainfall with moderate complexity is necessary for dynamic analysis of hydrological process. From the mathematical perspective, the rainfall intensity is a nonsmooth function of time [16]; that is, a pulsing function which means the time derivative is not continuous everywhere. In general, nonsmooth modeling has been widely used in many different areas [17–19] for its intuitional advantage; although nonsmooth terms can introduce additional mathematical difficulties for dynamic analysis [20, 21]. In hydrology, a simplified nonsmooth function, that is, a rectangular pulsing function, is frequently used for theoretical analysis of ecohydrological or hydrometeorologic processes. For example, a simple step function was used to describe the pulsing rainfall forcing by Collins and Bras [22] to examine how the gradient of mean annual precipitation is weaved in the topography of water-limited ecosystems; an ideal quasiperiodic pulsing precipitation was utilized by Lu [23] to study the physical basis of daily flood and drought index; a periodic step function was adopted by Robinson and Sivapalan [15] to identify different hydrological regimes according to the interactions between rainfall variability and catchment response time within a theoretical framework. Following these studies, a quasiperiodical rectangular pulsing function was chosen to represent real stochastic rainfall process by keeping between-storm temporal structure while ignoring more variability for a preliminary dynamic analysis of ecohydrological system.

Remarkably, the chosen pulsing function is just like the water application manner in an automatic irrigation system under the equilibrium state, see also Section 3 and Figure 8(c). The so-called automatic irrigation [24] means the duration of irrigation is controlled by the real-time observation on soil moisture content, with the switching bounds predetermined. As the dynamic analysis in this study focuses on the equilibrium state of ecohydrological system, we will use automatic irrigation scheme to represent the quasiperiodical pulsing function. Mathematically, such substitution of automatic irrigation scheme for stochastic rainfall process can ensure that the model is deterministic and autonomous in the resultant ordinary differential equations (ODEs), and therefore the dynamic phenomena discovered in the model depends not on any particular selection of external rainfall time series but on the internal system configuration. Practically, the automation has been widely used in water-saving irrigation systems [25, 26] under dry as well as wet conditions [27], and therefore we can assume a no-rainfall area and focus on the automatic irrigation system only, for the irrigation amount is much larger than the rainfall in the arid especially hyperarid regions [28].

2.2. Ecohydrological Model

An ecohydrological model is then introduced as follows. The schematic diagram is shown as Figure 1. As mentioned above, it can be assumed to be in the hyperarid area where water supply comes from irrigation alone. Our work imitates the tank model [29] to describe the storage and release of aquifer.

Figure 1 

The schematic diagram of the model.

As shown in Figure 1, the height ℎ(𝑡)(cm) stands for the storage per unit area at time 𝑡(d). Roughly speaking, it is equivalent to water content in the soil. The tank has only one side outlet with height 𝐻(cm). The side runoff from tank is assumed to be linear with respect to tank water storage excess relative to initial height 𝐻, defined by a positive coefficient 𝑐. The surface runoff is ignored. The term 𝑣(𝑡) (quantity of grasses per square meter) stands for biomass, that is, the vegetation density per unit surface area at time 𝑡. Generally, the evaporation and transpiration are related to biomass by 𝛼(ℎ,𝑣)𝑣. Assume that 𝛼 is constant for a preliminary study here. Then, the dynamical equation of water storage can be written as follows: ̇𝑐()(𝑡)̇𝑐()(𝑡)ℎ=𝑝−𝛼𝑣−ℎ−𝐻,ifℎ≥𝐻,0,ifℎ<𝐻,if𝑖=1,ℎ=−𝛼𝑣−ℎ−𝐻,ifℎ≥𝐻,0,ifℎ<𝐻,if𝑖=−1. (1) Here we divide the whole process into two phases and introduce a variable𝑖(𝑡) with value domain {−1,+1} to depict them; that is, 𝑖(𝑡)=1 means irrigation at time 𝑡, with constant intensity 𝑝(cm/d), while 𝑖(𝑡)=−1 means no irrigation at time 𝑡, that is, irrigation intensity is 0.

Another important dynamical equation is the evolution equation of the vegetation, which can be described as 𝑣̇=𝛽(ℎ,𝑣)𝑣, (2) where 𝛽(ℎ,𝑣) indicates the varying rate per unit vegetation. It is a complicated problem to describe 𝛽 function in ecology. However, a qualitative function can be described approximatively as follows: ()=with𝛽ℎ,𝑣−𝛽₀ℎ,if𝑣∈0,𝜇₁−𝛿₁,𝛽₀+𝛽₁𝛿₁+𝛿₂ℎ𝑣−𝜇₁ℎ,if𝑣∈𝜇₁−𝛿₁,𝜇₁+𝛿₂,𝛽₁ℎ,if𝑣∈𝜇₁+𝛿₂,𝜇₂−𝛿₃,−𝛽₀+𝛽₁𝛿₃+𝛿₄ℎ𝑣−𝜇₂ℎ,if𝑣∈𝜇₂−𝛿₃,𝜇₂+𝛿₄,−𝛽₀ℎ,if𝑣∈𝜇₂+𝛿₄,+∞. (3) The framework is similar to (the linearization of) the work in [30]. It means that the plants would die down when ℎ/𝑣 is too large or too small, that is, 𝑣̇<0 [8]. A sketch map of relation between “breed” rate and available water amount (for unit plant) is shown as Figure 2. Here 𝛽₁ (d⁻¹) is the maximum growth rate per unit grass, while −𝛽₀ (d⁻¹) is the rate of natural death. Also, there are two equilibrium states ℎ/𝑣=𝜇₁(cm/n), and ℎ/𝑣=𝜇₂(cm/n) during which the grass density keeps constant. Positive coefficients 𝛿_𝑘(cm/n), 𝑘=1,2,3,4, are used for separating the suitable growth parameter region.

Figure 2 

The schematic diagram of an ideal 𝛽 function.

2.3. Automatic Irrigation Scheme

To summarize, the model has three critical variables: ℎ(𝑡), 𝑣(𝑡), and 𝑖(𝑡). The value sets of ℎ(𝑡) and 𝑣(𝑡) are continuous with their dynamic functions (1)–(3), while the value set of 𝑖(𝑡) has only two discrete elements 1 and −1. We only have to define a dynamic function of 𝑖(𝑡) to complete our model. As mentioned in Section 2.1, we use a scheme of automatic irrigation switching rule based on the soil moisture content indicated by ℎ: 𝑖𝑡⁺=−𝑖(𝑡),if𝑖(𝑡)=−1,ℎ(𝑡)≤𝑤₀,or𝑖(𝑡)=1,ℎ(𝑡)≥𝑤₁,𝑖(𝑡),else, (4) where 𝑖(𝑡⁺)=lim_Δ𝑡→0⁺𝑖(𝑡+Δ𝑡). Bounds of water content 𝑤₀(cm) and 𝑤₁(cm) are introduced, 𝑤₀<𝑤₁. In agricultural water management practice, 𝑤₀ could be chosen near the wilting point, and 𝑤₁ could be chosen near the field capacity [31]. Automatic control can be realized by combining the equipments of soil moisture monitoring, pipe irrigation, and switching [25]. Finally, an executable conceptual ecohydrological model is obtained by combining (1)–(4) together.

2.4. Qualitative Analysis on the Model

Plotting a trajectory {(ℎ(𝑡),𝑣(𝑡),𝑡)∣𝑡⩾0} starting from (ℎ(0),𝑣(0)) in ℎ-𝑣 space, we get a phase diagram; for example, see Figure 3. Thereinto, the state point is the pair (ℎ(𝑡),𝑣(𝑡)), and the limit stationary state (the state when 𝑡→+∞) appears as a period “circle” motion between ℎ∈[𝑤₀,𝑤₁]. Generally, there are two stable stationary states. 𝑣(𝑡)≡0 is always a stationary state, which means that given plants cannot be generated without progenitor. Note the motion {(𝑤₀,0)→(𝑤₁,0)→(𝑤₀,0)} as the zero circle. However in most cases there should be another periodic orbit with 𝑣≠0 for some parameters and initial values. It means that plants can be grown steadily under certain conditions. Note this stable stationary state as the nonzero circle. Because the vegetation density changes periodically, we concern the average vegetation density in a period of nonzero circle, noted as 𝑣. Furthermore, note the average water supply and the average outflow in a nonzero circle as 𝑝 and 𝑄. They are two important efficiency indices.

Figure 3 

Typical multiple stationary states.

Efforts for theoretical analysis can be done after sufficient simplification of (1)–(3). For example, the Fourier series of the irrigation intensity of the limit cycle can be computed as follows (assume that 𝑖(0⁻)=−1 and 𝑖(0)=1): 𝑝(𝑡)=𝑝𝜂+2Σ^∞_𝑛=1sin𝑛𝜋𝜂cos𝑛𝜋𝜂−2𝑛𝜋𝜂𝑡𝑇, (5) where 𝑇 is the time length of a period, and 𝜂 is the proportion of irrigating time length in a period. Also restricting ℎ<𝐻,𝑣<𝜇₁, for all ℎ,𝑣, the first-order derivative of the so-called Poincaré map [2] at (𝑤₀,𝑣) is 𝑓^′(𝑣)|_(𝑤0,𝑣)=𝑤₁𝛼𝑤₁−𝑤₀𝛽₀+𝛽₁𝛿₁+𝛿₂+𝑤₁𝛼𝑝−𝛼𝑣²𝛽₀+𝛽₁𝛿₁+𝛿₂²𝛼𝑤₀𝑣−𝜇₁𝑝. (6) This map indicates the iterative mapping of vegetation density after a loop of irrigation-nonirrigation. It transforms the different equations into the iterative equations and may be convenient to solve some problems; for example, the stationary state of the differential equations corresponds to the equilibrium of 𝑓(𝑣). But generally, solving 𝑓(𝑣) without any constraint is quite difficult and beyond the scope of this study. Numerical experiment is then carried out to explore the potential nontrivial dynamic behaviors.

3. Analysis of Nontrivial Dynamic Behaviors through Numerical Experiments

Generally speaking, for a given research area and particular vegetation, 𝑐,𝐻,𝛼,𝛽,𝛿_𝑖,𝑖=1,2,3,4,𝜇₁,𝜇₂ are constants, and 𝑝,𝑤₀,𝑤₁ are artificial parameters. Furthermore, the adjustment range of 𝑤₀ and 𝑤₁ is quite small for they should be near to the wilting point and field capacity, respectively, in the irrigation application. Of course, they can also be extended when we move beyond the irrigation problem, see also the discussion in the conclusion section. During the following discussion, we choose 𝛼=1(cm⋅n⁻¹⋅d⁻¹),𝑐=0.5(d⁻¹),  𝐻=5(cm),𝛽₀=1(d⁻¹),𝛽₁=0.5(d⁻¹),𝛿₁=𝛿₃=0.5(cm⋅n⁻¹),𝛿₂=𝛿₄=1(cm⋅n⁻¹),𝜇₁=2(cm⋅n⁻¹), and 𝜇₂=5(cm⋅n⁻¹) as a numerical testing environment, tune 𝑤₀(cm) and 𝑤₁(cm) coarsely, and consider the dynamic behaviors as depending on fine tuning of 𝑝(cm/d). (Here 𝑛 stands for the quantity of grasses.)

There are three stationary states in Figure 3. Here (𝑤₀,𝑤₁,𝑝)=(1,7,4). The two thick “circles” are the two stable stationary states, and the thin line divide the first quadrant into two domains of attraction. Points starting from the upper region converge to the upper thick circle, while points starting from the lower region converge to the zero circle. Furthermore, the thin line between [𝑤₀,𝑤₁] is also a half part of a stationary state (from ℎ=𝑤₀ to ℎ=𝑤₁), and it is unstable, that is, any tiny disturbance would lead the trajectory to leave this stationary state.

Figure 4 shows some other types of the nonzero stable states for different parameters. The hysteresis is similar with the famous LotkaVolterra predator-prey models raised in the First World War [32]; that is, during different regions, vegetation intensity and soil moisture content can either be of positive correlation or negative correlation, similar to the relationship between cartilaginous fishes and food fishes in the Mediterranean.

(a) (𝑤0,𝑤1,𝑝)=(1,4,4)

(b) (𝑤0,𝑤1,𝑝)=(1,7,3)

Figure 4 

Different types of multiple stationary states.

Also for some arbitrary parameters 𝑤₀,𝑤₁,𝑝, the nonzero stable circle vanishes, which means all plants will die after sufficiently long time. The alteration of the number of stationary states is identified as typical bifurcation. As other parameters fixed, the bifurcation map with 𝑝 and 𝑣 can be drawn as in Figure 5. Bifurcations with other variables are omitted for conciseness. From the bifurcation map, we can see, as 𝑝 increased, that the steady vegetation coverage vanishes at 𝑝≈4.4 (Figure 5(a)) and 𝑝≈8.3 (Figure 5(b)), respectively. That means overlarge irrigation intensity will lead plants to die out.

(a)

(b)

(c)

Figure 5 

Bifurcation map with 𝑣 and 𝑝. Here (𝑤0,𝑤1)=(1,7) in (a) and (𝑤0,𝑤1)=(1,4) in ((b), (c)). (c) is a locally magnified (b).

Another important dynamic phenomenon presented in Figure 5 is the mutation, or so-called jump. The long-term behavior of our ecohydrological system is continuous with the parameter 𝑝 for most of the time but jumps at some individual points. For example, in Figure 5(a), 𝑣 varies continually with 𝑝 for most time and jumps near 𝑝=3.7,3.9,4.4. So along with the continuous increase of irrigation intensity, the limit vegetation density may suddenly go up or drop down discretely. A rough interpretation of this phenomenon is the mutation parameters collide the border of continuity intervals of the vegetation growth Equations (2) and (3) [21].

Further experiments are worthy to be done to verify the reliability of numerical analysis. In nature, the independent evolution events with the same parameters will never be repeated, and therefore these dynamical phenomena could not be directly verified by field observations. Nevertheless, the independent repeated artificial experiments could be conducted for the demonstration of these phenomena including the automatic irrigation experiments.

Analysis on the nontrivial dynamic behaviors can be used to guide human practice. For example, under the text condition, we plot the bifurcation maps 𝑣,𝑝,𝑄 (stand for the average vegetation coverage, irrigation quantity, and lateral flux, resp.) with 𝑝 in a same figure, see Figure 6. From the computation we can see that the maximal stable vegetation coverage is at 𝑝≈3.8, and the unit water consumption is also the least economical at this time.

Figure 6 

Bifurcation map of 𝑣,𝑝,and𝑄 with 𝑝.

Another typical example in which dynamic research plays a key role in practical problems is to qualitative study the influence of catastrophic mutation. For example, a fire can push the system from the attracting basin of the nonzero circle to the attracting basin of the zero circle, such as Figure 7.

Figure 7 

Influence of catastrophic mutation—a fire case.

(a)

(b)

(c)

Figure 8 

(a) The chaotic phase diagram; (b) The Lyapunov exponent. Here (𝑤0,𝑤1,𝑝)=(1,4,8.2); (c) the irrigation time series.

Moreover, deterministic chaos also occurs in this model, for example, when 𝑝∈[7.8,8.3], as shown in Figures 5(b) and 5(c). The chaotic limit set in ℎ-𝑣 space when 𝑝=8.2 is shown in Figure 8(a). In the scenario, with the constant irrigation strategy, the vegetation would not converge to a unique circle and the vegetation is “fluctuating” forever. Meanwhile, some chaotic characteristic can be computed. Consider the Poincaré map at the point (𝑤₀,𝑣), the Lyapunov exponent [2] of the Poincaré map is larger than zero, as shown in Figure 8(b), which confirms the existence of chaos. The time series of irrigation can be drawn as in Figure 8(c). To be noted, the graph of irrigation time series indicates a quasiperiodic rectangular pulsing function (see Section 2.1). As mentioned above, planar chaos is a specific behavior of nonsmooth dynamic system, which maybe more natural to depict hydrology process. The possible cause of this behavior is that the observation of soil water content of automatic irrigation system is discrete. The scenario should be designed (or avoided) to depend on the purpose of the irrigation system in the practice.

4. Conclusion

An ideal ecohydrological model convenient for dynamic analysis was established in this paper. Compared with traditional models, the nonsmooth property is emphasized. Pulsing function was used to describe the rainfall intensity with a practical interpretation of automatic irrigation. Under this approach, the deterministic process-based ecohydrological model was analyzed theoretically. Based on the model, some deterministic nontrivial dynamic phenomena are revealed. For example, the existence of multiple stationary states, mutation, bifurcation, and chaos was discussed under numerical simulation, with an overview on their possible application in optimization, avoiding the influence of catastrophic mutation, and so on.

This work can be seen as a kind of demonstration to show the existence of some complicated dynamic phenomena in ecohydrological process. Further observational study should be done to determine the scope of parameters in the real world which has nontrivial dynamic phenomena, probably starting from hyperarid area just as Liu et al. [33] have done to validate ecohydrological feedbacks in hyperarid Tarim River. Noticing that for the convenience of theoretical study, we used a step function to indicate the atmospheric forcing term, instead of considering the stochastic features of rainfall. Actually, one can change the switching policy (4) to a more complicated one, for example, through fixing the time length of 𝑖(𝑡)=1 corresponding to a timing irrigation, or randomizing the time length of 𝑖(𝑡)=1 corresponding to a stochastic irrigation. A bold hypothesis is that if we understand the rainfall mechanism thoroughly, a proper switching policy can be invented to simulate the fully stochastic features of real rainfall process by introducing additional oceanic and climatic variables, as implicated by the simple formulation of local recycling rainfall in [12, 13]. Also, the model could be further improved by refining the existing feedbacks between hydrological and ecological processes and by incorporating more feedbacks among hydrological, atmospheric, and other related processes.

Nomenclature

𝑐:	Side outlet coefficient (d−1)
ℎ:	Storage of aquifer (cm)
𝐻:	Height of the side outlet of the tank (cm)
𝑖:	Irrigation state 𝑖∈{±1}
𝑝:	Irrigation intensity (cm/d)
𝑄:	Side outlet intensity (cm/d)
𝑡:	Time (d)
𝑇:	Time length of a period (d)
𝑣:	Biomass (n) ((n) is short for quantity of grasses)
𝑤0,𝑤1:	Bounds of storage of aquifer (cm)
𝛼:	Transpiration coefficient (cm/(d⋅n))
𝛽0,𝛽1:	Bounds of evolution rate (d−1)
𝛿1∼4:	Distance between kickpoints and equilibria in (3) (n)
𝜇1,𝜇2:	Equilibria of evolution (n)
𝜂:	The proportion of irrigating time length in a period (100%)
𝑥̇:	d𝑥/d𝑡,𝑥 stands for any variable
𝑥:	Average 𝑥 in a nonzero circle, 𝑥 stands for any variable.

Acknowledgments

This research was funded by the National Science Foundation of China (NSFC51179084, 51190092, 51222901, and 11072274), the foundation of State Key Laboratory of Hydroscience and Engineering of Tsinghua University (2012-KY-03), and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.