Research article Special Issues

An integrated approach to modeling changes in land use, land cover, and disturbance and their impact on ecosystem carbon dynamics: a case study in the Sierra Nevada Mountains of California

  • Received: 23 February 2015 Accepted: 14 June 2015 Published: 16 June 2015
  • Increased land-use intensity (e.g. clearing of forests for cultivation, urbanization), often results in the loss of ecosystem carbon storage, while changes in productivity resulting from climate change may either help offset or exacerbate losses. However, there are large uncertainties in how land and climate systems will evolve and interact to shape future ecosystem carbon dynamics. To address this we developed the Land Use and Carbon Scenario Simulator (LUCAS) to track changes in land use, land cover, land management, and disturbance, and their impact on ecosystem carbon storage and flux within a scenario-based framework. We have combined a state-and-transition simulation model (STSM) of land change with a stock and flow model of carbon dynamics. Land-change projections downscaled from the Intergovernmental Panel on Climate Change’s (IPCC) Special Report on Emission Scenarios (SRES) were used to drive changes within the STSM, while the Integrated Biosphere Simulator (IBIS) ecosystem model was used to derive input parameters for the carbon stock and flow model. The model was applied to the Sierra Nevada Mountains ecoregion in California, USA, a region prone to large wildfires and a forestry sector projected to intensify over the next century. Three scenario simulations were conducted, including a calibration scenario, a climate-change scenario, and an integrated climate- and land-change scenario. Based on results from the calibration scenario, the LUCAS age-structured carbon accounting model was able to accurately reproduce results obtained from the process-based biogeochemical model. Under the climate-only scenario, the ecoregion was projected to be a reliable net sink of carbon, however, when land use and disturbance were introduced, the ecoregion switched to become a net source. This research demonstrates how an integrated approach to carbon accounting can be used to evaluate various drivers of ecosystem carbon change in a robust, yet transparent modeling environment.

    Citation: Benjamin M. Sleeter, Jinxun Liu, Colin Daniel, Leonardo Frid, Zhiliang Zhu. An integrated approach to modeling changes in land use, land cover, and disturbance and their impact on ecosystem carbon dynamics: a case study in the Sierra Nevada Mountains of California[J]. AIMS Environmental Science, 2015, 2(3): 577-606. doi: 10.3934/environsci.2015.3.577

    Related Papers:

    [1] Shengyu Huang, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Dynamics of a harvested cyanobacteria-fish model with modified Holling type Ⅳ functional response. Mathematical Biosciences and Engineering, 2023, 20(7): 12599-12624. doi: 10.3934/mbe.2023561
    [2] Moitri Sen, Malay Banerjee, Yasuhiro Takeuchi . Influence of Allee effect in prey populations on the dynamics of two-prey-one-predator model. Mathematical Biosciences and Engineering, 2018, 15(4): 883-904. doi: 10.3934/mbe.2018040
    [3] Shengyu Huang, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Dynamic analysis of a modified algae and fish model with aggregation and Allee effect. Mathematical Biosciences and Engineering, 2022, 19(4): 3673-3700. doi: 10.3934/mbe.2022169
    [4] Md Nazmul Hassan, Kelsey Thompson, Gregory Mayer, Angela Peace . Effect of Excess Food Nutrient on Producer-Grazer Model under Stoichiometric and Toxicological Constraints. Mathematical Biosciences and Engineering, 2019, 16(1): 150-167. doi: 10.3934/mbe.2019008
    [5] Zhenliang Zhu, Yuming Chen, Zhong Li, Fengde Chen . Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey. Mathematical Biosciences and Engineering, 2023, 20(6): 10977-10999. doi: 10.3934/mbe.2023486
    [6] Santanu Bhattacharya, Nandadulal Bairagi . Dynamic optimization of fishing tax and tourism fees for sustainable bioeconomic resource management. Mathematical Biosciences and Engineering, 2025, 22(7): 1751-1789. doi: 10.3934/mbe.2025064
    [7] Wenjie Yang, Qianqian Zheng, Jianwei Shen, Linan Guan . Bifurcation and pattern dynamics in the nutrient-plankton network. Mathematical Biosciences and Engineering, 2023, 20(12): 21337-21358. doi: 10.3934/mbe.2023944
    [8] Huanyi Liu, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Dynamical analysis of an aquatic amensalism model with non-selective harvesting and Allee effect. Mathematical Biosciences and Engineering, 2021, 18(6): 8857-8882. doi: 10.3934/mbe.2021437
    [9] Eric M. Takyi, Charles Ohanian, Margaret Cathcart, Nihal Kumar . Dynamical analysis of a predator-prey system with prey vigilance and hunting cooperation in predators. Mathematical Biosciences and Engineering, 2024, 21(2): 2768-2786. doi: 10.3934/mbe.2024123
    [10] Yong Yao . Dynamics of a delay turbidostat system with contois growth rate. Mathematical Biosciences and Engineering, 2019, 16(1): 56-77. doi: 10.3934/mbe.2019003
  • Increased land-use intensity (e.g. clearing of forests for cultivation, urbanization), often results in the loss of ecosystem carbon storage, while changes in productivity resulting from climate change may either help offset or exacerbate losses. However, there are large uncertainties in how land and climate systems will evolve and interact to shape future ecosystem carbon dynamics. To address this we developed the Land Use and Carbon Scenario Simulator (LUCAS) to track changes in land use, land cover, land management, and disturbance, and their impact on ecosystem carbon storage and flux within a scenario-based framework. We have combined a state-and-transition simulation model (STSM) of land change with a stock and flow model of carbon dynamics. Land-change projections downscaled from the Intergovernmental Panel on Climate Change’s (IPCC) Special Report on Emission Scenarios (SRES) were used to drive changes within the STSM, while the Integrated Biosphere Simulator (IBIS) ecosystem model was used to derive input parameters for the carbon stock and flow model. The model was applied to the Sierra Nevada Mountains ecoregion in California, USA, a region prone to large wildfires and a forestry sector projected to intensify over the next century. Three scenario simulations were conducted, including a calibration scenario, a climate-change scenario, and an integrated climate- and land-change scenario. Based on results from the calibration scenario, the LUCAS age-structured carbon accounting model was able to accurately reproduce results obtained from the process-based biogeochemical model. Under the climate-only scenario, the ecoregion was projected to be a reliable net sink of carbon, however, when land use and disturbance were introduced, the ecoregion switched to become a net source. This research demonstrates how an integrated approach to carbon accounting can be used to evaluate various drivers of ecosystem carbon change in a robust, yet transparent modeling environment.


    In eutrophic lakes, cyanobacterial Microcystis colonies usually aggregate and float upwards to form nuisance mucilaginous Microcystis blooms, the development of mucilaginous cyanobacterial Microcystis blooms is a serious environmental and ecological problem, and information on the bloom-formation mechanism has been lacking until now [1]. Based on their long-term research work on cyanobacteria bloom in Taihu Lake, Qin et al. [2] believed that the outbreak of cyanobacteria bloom actually experienced four stages: 1) cyanobacteria cell proliferation stage; 2) cyanobacteria cell groups formation stage; 3) cyanobacteria cell groups floating stage; 4) cyanobacteria bloom outbreak stage. At the same time, they also deemed that in the process of cyanobacteria cell aggregation, the collision probability of cell groups could increase, and smaller cell masses were easier to bond, which could form larger groups and accelerate the emergence of cyanobacteria bloom. Furthermore, Chen et al. [3] pointed out that the aggregation and migration of Microcystis had many ecological advantages, which could promote the formation of Microcystis bloom. Kong et al. [4] proposed that the most direct reason for the formation of cyanobacteria bloom was the massive aggregation, growth, migration and floating to the water surface stage of Microcystis. Obviously, the formation of cyanobacteria cell aggregation is one of the keys to the outbreak of cyanobacteria bloom. Therefore, how to prevent and disperse cyanobacteria cell groups is one of the key control measures to prevent cyanobacteria bloom.

    Cyanobacterial blooms are commonly referred to as harmful algal blooms (HABs) due to their vast ecological impacts and ability to generate metabolites and taste, a variety of in-lake/reservoir control measures are implemented to reduce the abundance of nuisance cyanobacteria biomass[5]. Kibuye [5,6] gave a critical review on operation and performance of source water control strategies for cyanobacterial bloom: chemical control methods, mechanical/physical control methods and biological control methods. And they pointed out that each control method had site-specific strengths, limitations, and ecological impacts, which should hint that none of the reviewed control strategies can provide a comprehensive solution to cyanobacterial blooms. In view of the cyanobacteria blooms in Taihu Lake, the competent authorities mainly adopt the comprehensive method of combining mechanical control method and biological control method, this is because that Taihu Lake is a drinking water source. Although most mechanical and biological control strategies can offer long-term control, their application can be cost-prohibitive and treatment efficacy is influenced by source water geometry and continual nutrient inputs from external sources[6]. Therefore, how to coordinate mechanical control methods with biological control methods is very important, and it is also one of the problems worthy of our in-depth study.

    In order to curb the scale of cyanobacteria bloom in Taihu Lake, the management department of Taihu Lake should realize the transformation from focusing on source control and pollution interception to paying equal attention to source control, pollution interception and ecological regulation[7]. At the same time, there are many kinds of cyanobacteria in Taihu Lake, usually Microcystis[8], and allowing cyanobacteria cells to aggregate is one of the necessary conditions for cyanobacteria bloom [9]. Furthermore, once the cyanobacteria bloom breaks out, the regeneration from cell turnover and nutrient recycling by closely associated heterotrophic bacteria and microzooplankton grazers can help sustain bloom biomass, so the key biological factors to control cyanobacteria bloom must include grazing of zooplankton (possibly benthos and fish)[10]. However, some literatures have studied the dynamic relationship between cyanobacteria cell aggregation and zooplankton/fish grazing. Yang et al.[11] investigated the effect of grazing by different sorts of zooplankton on the induction of defensive morphology in the cyanobacterium Microcystis aeruginosa, and pointed out that cyanobacterium Microcystis aggregation could effectively block zooplankton predation and thus increase the survival of Microcystis aeruginosa. Burkert et al.[12] inquired into the interactions between the cyanobacteria and the potential grazer, pointed out that the predator came in direct contact with Microcystis, aggregates were formed. Chow-Fraser et al.[13] had evidence that grazers in oligotrophic lakes exert a greater impact on algae than those in eutrophic lakes. To summarise, cyanobacteria species can also benefit by their tendencies to congregate as large filamentous and colonial colonies, which reduces zooplankton predation.

    At present, the cyanobacteria bloom research team of Wenzhou University is trying to use a novel comprehensive method to deal with the problem of cyanobacteria bloom in Taihu Lake, which is a combination of mechanical control method and biological control method. Firstly, in the initial stage of cyanobacteria bloom, a part of cyanobacteria agglomerates were broken by physical jet treatment device (a cyanobacteria bloom treatment ship is developed by Wenzhou University), its purpose is to destroy the aggregation stage of cyanobacteria cells and improve the chance of individual cyanobacteria cells being grazed. Secondly, some streamline ecosystems and food web structures will be reconstructed to create unfavorable conditions for cyanobacteria, some filter-feeding organisms (such as zooplankton, mussels, bivalves and fish) can be added to feed on cyanobacteria [14,15,16]. Finally, with the help of the internal operation characteristics of aquatic ecosystem, we expect that aquatic ecosystem can quickly change from algae bloom ecosystem to algae-zooplankton-fish cycle ecosystem. However, whether or not cyanobacterial bloom can be successfully controlled, the novel comprehensive method can influence changes in competition, predation, and behavioral response to predation and competition in the restored habitat. Therefore, after the cyanobacteria bloom is treated by physical jet technology, the restoration and increase of zooplankton and filter-feeding fish how to affect the self purification ability and balance of aquatic ecosystem are worthy of our further study, whether algal population can coexist with zooplankton or filter-feeding fish in the range of relatively low algal population is also worth exploring.

    With the rapid development of numerical analysis and simulation technology, dynamic models have gradually become one of the powerful tools for studying natural science phenomena and problems, which can promote the relevant scientific problems to obtain some excellent research results, such as relevant literatures [17,18,19,20,21,22,23,24,25,26,27]. In references [17,18,19] applied reaction-diffusion equation to examine the calcium diffusion in the cells, investigated some nonlinear dynamics problems and got some interesting results. In reference [20] explored some dynamical behaviors of forced KdV equation to describe the free surface critical flow over a hole, and some meaningful results were given. In reference [21] inquired into chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells, some meaningful results showed that order of the fractional derivative had a significant effect on the dynamics process. The paper [22] studied the dynamics of infectious diseases within a human host and in the population, and gave some meaningful theoretical and numerical simulation results. In reference [23] gave some good results to verify the effectiveness of linear feedback control for chaos. In reference [24] investigated a newly developed model for Hepatitis-B infection in sense of the Atangana-Baleanu Caputo fractional-order derivative, and obtained some good research conclusions. In references [25,26,27] looked into multiple bifurcations of discrete-time dynamic model, and some important theoretical and numerical results are obtained. In a word, the research on dynamic models has been developed rapidly, and a large number of excellent research results have emerged.

    Now, the disaster of aquatic ecosystem caused by cyanobacteria bloom is one of the hot issues in the world. After the cyanobacteria bloom is treated by chemical control methods, mechanical control methods and biological control methods, the essential characteristics of the internal self recovery and strengthening mechanism of aquatic ecosystem is one of the problems that need to be studied urgently. Furthermore, the algae bloom research team of Wenzhou University has built a physical spraying treatment ship, which mainly deals with the algae blooms in subtropical lakes and reservoirs, and reconstructs the biological food chain for ecological algae control after treatment. Thus, in order to understand the internal self recovery and strengthening mechanism of aquatic ecosystem, we must deeply explore the interactions between cyanobacteria (Microcystis aeruginosa) and potential grazer after cyanobacteria bloom treatment, and make clear the coexistence mode and change trend of cyanobacteria (Microcystis aeruginosa) and potential grazers based on water ecological model and dynamic simulation. Therefore, we believe that the originality of this paper is to illustrate the dynamic evolution relationship between cyanobacteria and potential grazers by means of the bifurcation dynamic evolution process of the model $ (2.1) $, and give the evolutionary characteristics of the coexistence mode of cyanobacteria and potential grazers.

    Bertalanffy [28] proposed a method to study the biological system using mathematical model in 1932, this approach is a combination of coordination, order and purpose, which can form three basic ideas of studying the biological system, namely: trophic-level analysis, system perspective and dynamic view. Since then, the research on ecological model and its related dynamics has developed unprecedentedly, and a large number of literature and works have appeared, such as fundamentals of ecological modelling [29], a nitrogen-phosphorus-blue-green algae model [30], some nitrogen and phosphorus dose-response models [31], a mathematical model including three types of zooplankton and two types of fish as well as two types of algae and nutrients [32], an aquatic bacteria-algae amensalism model [33], an ecological model considering two types of phytoplankton, three types of zooplankton, planktivorous fish, detritus and dissolved matters [34], an aquatic algae-fish model [35], one dimensional hydrodynamic-ecological model [36], a two-component model with nutrients and cyanobacteria [37], a four-component model including nutrient, unicellular algae, colonial algae and herbivorous zooplankton [38], an algal bloom mathematical model involving zooplankton-phytoplankton population [39]. Obviously, the ecological model has been developed rapidly in water eutrophication and algal bloom, and some meaningful theoretical results have been obtained. Although there is a certain understanding of the mechanism and control effect of cyanobacteria bloom in recent years, there are still many unsolved mysteries in the outbreak and control of cyanobacteria bloom.

    In order to make clear the coexistence mode and change trend of cyanobacteria (Microcystis aeruginosa) and potential grazers (mussels, bivalves, and silver carp) after cyanobacteria bloom is treated by comprehensive control strategy, which can help us explore the interactions between cyanobacteria and potential grazers, and explicit the effectiveness and feasibility of the comprehensive control strategy, we will construct a new aquatic ecological model based on the comprehensive control strategy and internal operation law of aquatic ecosystem. Therefore, some model assumptions and cyanobacteria bloom background will be given as follows:

    1) The growth environment of cyanobacteria and potential grazers is assumed to be a semi closed aquatic ecosystem, species biomass is always evenly distributed in space and changed instantaneously with time $ t $. $ x(t) $ and $ y(t) $ represent respectively the biomass of total cyanobacteria (unicellular cyanobacteria and colonial cyanobacteria) and potential grazers. After cyanobacteria bloom is treated by mechanical control method (physical jet processing device), total cyanobacteria will be divided into colonial cyanobacteria $ \alpha_{1} x(t) $ and unicellular cyanobacteria $ \alpha_{2} x(t) $ with $ \alpha_{1}+ \alpha_{2} = 1 $. If $ \alpha_{1} = 0 $ and $ \alpha_{2} = 1 $, this shows that cyanobacteria population mainly exists in the form of unicellular cyanobacteria, and the cyanobacteria population has no aggregation behavior and can not form a group aggregation state. If $ \alpha_{1} = 1 $ and $ \alpha_{2} = 0 $, this indicates that cyanobacteria bloom has erupted in a large area, especially serious, and all cyanobacteria populations are clustered, which is an extreme state of existence and almost does not exist in reality.

    2) Because when unicellular cyanobacteria gather into groups, the light utilization rate of algae cells can be enhanced [40], the salt stress resistance and high light stress of algal cells can be improved [41,42], the ability to avoid the invasion of viruses, bacteria and algae phages can be strengthened [43], nutrient and inorganic carbon absorption capacity can be greatly promoted[44,45], the risk of predation can be greatly reduced [46]. In light of the above, cyanobacteria aggregation can find more suitable niche, obtain greater ecological advantages and enhance growth rate. Therefore, we assume that the growth kinetics function of cyanobacteria $ x(t) $ is $ f(\frac{\alpha_{1}}{\alpha_{2}})x(t)(1-\frac{x(t)}{k}) $ with intrinsic growth rate function $ f(\frac{\alpha_{1}}{\alpha_{2}}) $ and maximum environmental capacity $ k $,

    $ f(\frac{\alpha_{1}}{\alpha_{2}}) = {r1α1α2,0<α1<1,0<α2<1,rmin,α1=0,α2=1,rmax,α1=1,α2=0,
    $

    where $ r_{1} $ is a median growth rate with $ \alpha_{1} = \frac{1}{2} $ and $ \alpha_{2} = \frac{1}{2} $, $ r_{min} $ is a growth rate without aggregation, and, $ r_{max} $ is a growth rate with overall aggregation state.

    3) During cyanobacteria blooms, cyanobacteria can induce morphological changes to colonial forms in order to avoid zooplankton grazing, this mechanism can be an adaptive reaction for survival that has developed through evolutionary processes[47,48]. For example, Microcystis exists as colonial forms, the relevant research results are in the papers[35,38,43]. However, unicellular cyanobacteria do not have this protective mechanism. Hence, we use Holling type II functional responses $ a_{2}\frac{\alpha_{2}x(t)y(t)}{b+\alpha_{2}x(t)} $ to describe the predation mechanism of potential grazers on unicellular cyanobacteria, this is because the predation process of potential grazers not only take time to capture, but also depend on the abundance of unicellular cyanobacteria, where $ a_{2} $ is a grazing coefficient and $ b $ is a half-saturation constant. At the same time, we use mathematical functional $ a_{1}(\alpha_{1}x(t)-g(\frac{\alpha_{1}}{\alpha_{2}}))y(t) $ to describe the predation mechanism of potential grazers on colonial cyanobacteria, $ a_{1} $ is a grazing coefficient and predation function $ g(\frac{\alpha_{1}}{\alpha_{2}}) $ avoided by colonial cyanobacteria, and

    $ g(\frac{\alpha_{1}}{\alpha_{2}}) = {mα1α2,0<α1<1,0<α2<1,0,α1=0,α2=1,mmax,α1=1,α2=0,
    $

    where $ m $ is a median biomass of colonial cyanobacteria (which can never be hunted) with $ \alpha_{1} = \frac{1}{2} $ and $ \alpha_{2} = \frac{1}{2} $, $ 0 $ indicates that unicellular cyanobacteria do not have this ability, $ m_{max} $ is a maximum biomass of colonial cyanobacteria (which can never be hunted) with overall aggregation state.

    4) It is widely known that potential grazers can not eat only cyanobacteria, and must forage other species as alternate foods. Thus, the growth kinetics function of potential grazers is $ r_{2}y(t)+e_{1}a_{1}(\alpha_{1}x(t)-g(\frac{\alpha_{1}}{\alpha_{2}}))y(t)+e_{2}a_{2}\frac{\alpha_{2}x(t)y(t)}{b+\alpha_{2}x(t)} $ with alternate coefficient $ r_{2} $ and absorption conversion coefficients $ e_{1} $ and $ e_{2} $. Furthermore, cyanobacteria and potential grazers have their own other losses (such as natural death), we might as well set them to $ m_{1}x(t) $ and $ m_{2}y(t) $ with loss coefficient $ m_{1}, m_{2} $.

    Based on the above assumptions and analysis, we can construct a new aquatic ecological control model and two special aquatic ecological models with non-negative initial conditions $ x(0) \equiv x_0 \geq 0 $ and $ y(0) \equiv y_0 \geq 0 $, which can be described as following:

    $ {dx(t)dt=r1α1α2x(1xk)a1(α1xα1α2m)ya2α2xyb+α2xm1x,dy(t)dt=r2y+e1a1(α1xα1α2m)y+e2a2α2xyb+α2xm2y,
    $
    (2.1)
    $ {dx(t)dt=rminx(1xk)a2xyb+xm1x,dy(t)dt=r2y+e2a2xyb+xm2y,
    $
    (2.2)
    $ {dx(t)dt=rmaxx(1xk)a1(xmmax)ym1x,dy(t)dt=r2y+e1a1(xmmax)ym2y.
    $
    (2.3)

    The model (2.1) represents the mixed state of unicellular cyanobacteria and colonial cyanobacteria, the biomass of unicellular cyanobacteria and colonial cyanobacteria is controlled by critical parameter $ \alpha_{1} $ or $ \alpha_{2} $ because of $ \alpha_{1}+\alpha_{2} = 1 $. If the value of $ \alpha_{1} $ belongs to $ (0, 0.2] $, the model (2.1) represents cyanobacteria cell proliferation stage. If the value of $ \alpha_{1} $ belongs to $ (0.2, 0.4] $, the model (2.1) represents cyanobacteria cell groups formation stage. If the value of $ \alpha_{1} $ belongs to $ (0.4, 0.6] $, the model (2.1) represents cyanobacteria cell groups floating stage. If the value of $ \alpha_{1} $ belongs to $ (0.6, 1) $, the model (2.1) represents cyanobacteria bloom outbreak stage. Overall, the model $ (2.1) $ can describe that the outbreak of cyanobacteria bloom in Taihu Lake actually experienced four stages [2]. The models (2.2) and (2.3) represent the stage of unicellular cyanobacteria and colonial cyanobacteria respectively, neither of these two states will exist in reality, therefore, we only study them as two limit states.

    With the help of bifurcation dynamic analysis and dynamic simulation of the model (2.1), the main purpose of this paper is to explore how the aquatic ecosystem can recover and strengthen itself after the comprehensive treatment of cyanobacteria bloom. To solve this problem, we firstly make a qualitative analysis of the model (2.1), deduce some critical threshold conditions to ensure the existence and stability of some equilibrium points, and give some critical conditions under which the model can induce Transcritical and Hopf bifurcation. Secondly, we will comprehensively investigate how comprehensive management measures affect the dynamic characteristics of cyanobacteria bloom and evaluate the effectiveness of comprehensive management measures in controlling cyanobacteria bloom. Finally, we carry out relevant dynamic experiments to verify the feasibility of the theoretical derivation results and identify the possible coexistence mode of cyanobacteria and potential grazers, then we can analyze interaction mechanism between cyanobacteria and potential grazers in the process of cyanobacteria bloom outbreak and control.

    In order to investigate the superiority of the model $ (2.1) $, the dynamic relationship between the biomass of cyanobacteria and potential grazers was qualitatively analyzed and numerically simulated with $ r_1 = 0.56 $, $ r_2 = 0.2 $, $ k = 5 $, $ a_1 = 0.6 $, $ a_2 = 0.4 $, $ b = 0.15 $, $ m_1 = 0.18 $, $ m_2 = 0.28 $, $ e_1 = 0.6 $, $ e_2 = 0.28 $, $ \alpha_{1} = 0.35 $, $ \alpha_{2} = 0.65 $, $ r_{min} = 0.2 $, $ r_{max} = 0.8 $, $ m_{max} = 3 $, $ m = 0.1 $.

    From the model (2.2), we can find that the relation expression of cyanobacteria $ x $ and potential grazers $ y $ is given by

    $ y = \frac{(r_{min} \ \ (1-\frac{x}{k})-m_{1})(b+x)}{a_{2}}. $

    It is easy to see from Figure 1(a) that the dynamic relationship of cyanobacteria $ x $ and potential grazers $ y $ is a concave function as the value of cyanobacteria $ x $ increases, and there is a unique maximum value, and the maximum value of cyanobacteria population is still low. Unfortunately, the value range of cyanobacteria $ x $ is relatively small, and the biomass of potential grazers $ y $ becomes negative with the increase of cyanobacteria $ x $. Obviously, this dynamic result can not better reflect the running law of cyanobacteria bloom control test in the laboratory.

    Figure 1.  The dynamic evolution relationship between cyanobacteria and potential grazers.

    From the model (2.3), we can find that the relation expression of cyanobacteria $ x $ and potential grazers $ y $ is given by

    $ y = \frac{a_{1}x+m_{1}x-r_{max} \ \ \ (1-\frac{x}{k})}{a_{1}m_{max}}. $

    It can clearly find from Figure 1(b) that the dynamic relationship of cyanobacteria $ x $ and potential grazers $ y $ is a convex function as the value of cyanobacteria $ x $ increases, and there is a unique minimum value $ 0 $, which is not in line with the wild growth law of cyanobacteria population because cyanobacteria population will never die out. Furthermore, it is even more regrettable that the value of potential grazers $ y $ approaches positive infinity with the increase of cyanobacteria $ x $ value, which completely violates the dynamic relationship between cyanobacteria $ x $ and potential grazers $ y $ in the process of cyanobacteria bloom.

    From the model (2.1), we can find that the relation expression of cyanobacteria $ x $ and potential grazers $ y $ is given by

    $ y = \frac{\frac{r_{1}\alpha_{1}x}{\alpha_{2}}-m_{1}x}{a_{1}(\alpha_{1}x-\frac{\alpha_{1}m}{\alpha_{2}})+\frac{a_{2}\alpha_{2}x}{b+\alpha_{2}x}}. $

    It is worth pointing out that potential grazers $ y $ is monotonically increasing and has supremum $ \frac{r_{1}}{a_{1}\alpha_{2}} - \frac{m_{1}}{a_{1}\alpha_{1}} $ as the value of cyanobacteria $ x $ increases to positive infinity. Furthermore, the better result is that cyanobacteria $ x $ approaches a small fixed value when potential grazers $ y $ approaches positive infinity, this result is more in line with the natural law (see Figure 1(c)). Moreover, what is more gratifying is that there is a unique minimum value, through which we can regulate the recyclability of cyanobacteria bloom test, which is also our favorite routine manipulation.

    In a word, through the comparative analysis of numerical simulation results, it can be seen that the model $ (2.1) $ can more accurately describe the dynamic relationship between cyanobacteria and potential grazers under the mechanical and ecological integrated control, so we mainly study the relevant dynamic problems of the model $ (2.1) $.

    The existence and local stability of all possible equilibrium points will be investigated in detail. The equilibrium point is a kind of special solution of the model (2.1), which can determine some coexistence modes of cyanobacteria and potential grazers, and it is also the basis of subsequent bifurcation dynamics research. Thus, the existence and local stability of all possible equilibrium points will be futher investigated in this section.

    Now, we consider the algae isocline vertically and potential grazers isocline horizontally, respectively:

    $ { r1α1α2x(1xk)a1α1(xmα2)ya2α2xyb+α2xm1x=0, r2+e1a1α1(xmα2)+e2a2α2xb+α2xm2=0&y=0.
    $
    (2.4)

    By solving the above equations with zero dash, there are two classes of equilibrium points: boundary equilibrium points $ E_{0} $ and $ E_{1} $, the coordinates of internal equilibria $ E_1^* $ and $ E_2^* $.

    It is obvious that the equilibrium points are the intersections of these nullclines. Thus, we easily see that the model (2.1) possesses two boundary equilibrium points $ E_{0}(0, 0) $ and $ E_{1}(x_1, 0) $. For the possible positive equilibria, we only need consider the positive solutions of the following equations:

    $ { a1e1α1α2x2+(α2r2α2m2+a2e2α2+a1be1α1a1e1α1m)x+b(r2m2a1e1α1bmα2)=0, r1α1α2x(1xk)a1α1(xmα2)ya2α2xyb+α2xm1x=0.
    $
    (2.5)

    For positive equilibrium points, the algae population $ x $ must satisfy $ \frac{m}{\alpha_{2}} < x < k $, Let $ \Delta(x) $ denote the discriminant of the first equation of (2.5) and express $ \Delta(x) $ in terms of, i.e,

    $ Δ(x)=B24AC,
    $

    where, $ B = (r_{2}-m_{2}+a_{2}e_{2})\alpha_{2}+a_{1}e_{1}\alpha_{1}b-a_{1}e_{1}\alpha_{1}m $,

    $ A = a_{1}e_{1}\alpha_{1}\alpha_{2} $, $ C = b(r_{2}-m_{2}-\frac{a_{1}e_{1}\alpha_{1}bm}{\alpha_{2}}) $.

    Hence, we have

    $ Δ=((r2m2+a2e2)α2+a1e1α1ba1e1α2m)24a1e1α1α2b(r2m2+a1e1m).
    $

    With respect to the conditions of the equilibrium points, we obtain the following theorem:

    Theorem 1. The equilibrium points of the model (2.1) are as follows: the model(2.1) always has a boundary equilibrium point $ E_{0}(0, 0) $. If $ \frac{\alpha_{1}}{\alpha_{2}} > \frac{m_{1}}{r_{1}} $, then the model (2.1) has a boundary equilibrium point $ E_{1}(x_{1}, 0) $ with $ x_{1} = \frac{k(r_{1}\alpha_{1}-m_{1}\alpha_{2})}{r_{1}\alpha_{1}} $.

    Theorem 2. For the positive equilibrium point, we have:

    1) When $ m\le I_{1} $, the model (2.1) has no positive equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $,

    2) When $ r_{2}-m_{2}+a_{2}e_{2} > 0 $ and $ r_{2}-m_{2} < 0 $, if $ max \{ 0, I_{1}\} < m < min \{ I_{2}, I_{3}, I_{4}\} $ or $ max \{ 0, I_{1}, I_{2}\} < m < min \{ I_{3}, I_{4}\} $, then the model (2, 1) has a unique positive equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $,

    3) When $ r_{2}-m_{2}+a_{2}e_{2} < 0 $, if $ max \{ 0, I_{1}, I_{3}\} < m < min \{ I_{2}, I_{4}\} $ or $ max \{ 0, I_{1}, I_{2}, I_{3}\} < m < I_{4} $, then the model (2.1) has a unique positive equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $,

    where

    $ x1=((r2m2+a2e2)α2+a1e1α1ba1e1α2m)+Δ2a1e1α1α2,y1=(r1α1α2m1)x1r1α1kα2x21a1(α1x1α1mα2)+a2α2x1b+α2x1,
    $
    $ I1=(r2m2)α2a1e1α1,I2=α2(r2m2+a2e2)α1a1e1,I3=b(m2r2)r2m2+a2e2,I4=(kα22+bα2)(r2m2+a1e1α22ka1e1α1(α2k+b).
    $

    Proof. Let the solution of the nullclines $ g(x) = Ax^2+Bx+C = 0 $, we know that the model (2.1) has the positive roots only if $ r_{2}-m_{2}-\frac{a_{1}e_{1}\alpha_{1}m}{\alpha_{2}} < 0 $, which means $ m\ge I_{1} $. At the same time, we also show that if $ m\ge I_{1} $, then $ \Delta > 0 $ and $ C < 0 $, which means that $ g(x) $ has a unique positive root $ E_1^*(x_{1}^*, y_{1}^*) $. Hence, if $ g(\frac{m}{\alpha_{2}}) < 0 $ and $ g(k) > 0 $, there is a positive real root in $ \frac{m}{\alpha_{2}} < x < k $. Due to the fact that $ \Delta > 0 $, hence, if $ g(\frac{m}{\alpha_{2}}) < 0 $ and $ g(k) > 0 $, we have the conditions with $ I_{3} $ and $ I_{4} $. If $ g(\frac{m}{\alpha_{2}}) < 0 $, it follows that $ r_{2}-m_{2}+a_{2}e_{2} > 0 $ for $ m < I_{3} $ or $ r_{2}-m_{2}+a_{2}e_{2} < 0 $ for $ m > I_{3} $. If $ m < I_{4} $, we have $ f(k) > 0 $. Next, we just discuss the symbols of B.

    If $ m = I_{2} $, we have $ B = 0 $, then $ g(x) = 0 $ has a unique positive root $ x_{1} = \sqrt{\frac{b(\alpha_{2}m_{2}-\alpha_{2}r_{2}+a_{1}e_{1}\alpha_{1}m)}{a_{1}e_{1}\alpha_{1}\alpha_{2}^2}} $. Similarly, if $ m < I_{2} $ or $ m > I_{2} $, then $ B > 0 $ or $ B < 0 $, we must know that $ g(x) = 0 $ has a unique positive root $ x^* $. From the above analysis, the model (2.1) has a unique positive equilibrium point if and only if one of the following situation:

    1) When $ m = I_{2} $, the model (2.1) has a unique equilibrium point $ E_{1}^*(x_{1}^*, y_{1}^*) $.

    2) When $ r_{2}-m_{2}+a_{2}e_{2} > 0 $ and $ r_{2}-m_{2} < 0 $, if $ max \{ 0, I_{1}\} < m < min \{ I_{2}, I_{3}, I_{4}\} $ or $ max \{ 0, I_{1}, I_{2}\} < m < min \{ I_{3}, I_{4}\} $, then the model (2.1) has a unique positive equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $.

    3) When $ r_{2}-m_{2}+a_{2}e_{2} < 0 $, if $ max \{ 0, I_{1}, I_{3}\} < m < min \{ I_{2}, I_{4}\} $ or $ max \{ 0, I_{1}, I_{2}, I_{3}\} < m < I_{4} $, then the model (2.1) has a unique positive equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $.

    Next, we will study the types and stability of the equilibrium points. Firstly, the Jacobian matrix of the model (2.1) is given as follows:

    $ J(x, y) = (r1α1α2(12xk)a1α1ya2α2by(b+α2x)2m1a1α1(xmα2)a2α2xb+α2xe1a1α1y+a2e2α2by(b+α2x)2r2m2+e1a1α1(xmα2)+e2a2α2xb+α2x)
    . $

    Theorem 3. With respect to the origin of this equilibrium point $ E_0 $, we have:

    1) When $ \frac{\alpha_{1}}{\alpha_{2}} < \frac{m_{1}}{r_{1}} $ and $ \frac{\alpha_{1}}{\alpha_{2}} < \frac{m_{2}-r_{2}}{e_{1}a_{1}m} $ hold, the equilibrium point $ E_{0} $ is a stable node.

    2) When $ \frac{\alpha_{1}}{\alpha_{2}} > \frac{m_{1}}{r_{1}} $ and $ \frac{\alpha_{1}}{\alpha_{2}} < \frac{m_{2}-r_{2}}{e_{1}a_{1}m} $ hold, or $ \frac{\alpha_{1}}{\alpha_{2}} < \frac{m_{1}}{r_{1}} $ and $ \frac{\alpha_{1}}{\alpha_{2}} > \frac{m_{2}-r_{2}}{e_{1}a_{1}m} $ hold, the equilibrium point $ E_{0} $ is a saddle.

    Proof. The Jacobian matrix of the equilibrium point $ E_{0} $ is

    $ J_{E_{0}} = (r1α1α2m1α1a1mα20r2m2+α1e1a1mα2)
    . $

    Obviously, $ J_{E_{0}} $ has two characteristic roots $ \lambda_{1} = \frac{r_{1} \alpha_{1}}{\alpha_{2}}-m_{1} $, $ \lambda_{2} = r_{2}-m_{2}+\frac{\alpha_{1}e_{1}a_{1}m}{\alpha_{2}} $. We applied the Routh-Hurwitz criterion to analyze the stability. The equilibrium point $ E_{0} $ is a stable node since $ \lambda_{1} < 0 $ and $ \lambda_{2} < 0 $, while$ \frac{r_{1}\alpha_{1}}{\alpha_{2}} < m_{1} $ and $ r_{2}+\frac{e_{1}a_{1}\alpha_{1}m}{\alpha_{2}} < m_{2} $ hold. The equilibrium point $ E_{0} $ is a saddle since $ \lambda_{1} < 0 $ and $ \lambda_{2} < 0 $, while $ \frac{r_{1}\alpha_{1}}{\alpha_{2}} < m_{1} $ and $ r_{2}+\frac{e_{1}a_{1}\alpha_{1}m}{\alpha_{2}} < m_{2} $ hold. Similarly it can be proved that the equilibrium point $ E_{0} $ is a saddle if $ \frac{r_{1}\alpha_{1}}{\alpha_{2}} < m_{1} $ and $ r_{2}+\frac{e_{1}a_{1}\alpha_{1}m}{\alpha_{2}} > m_{2} $.

    Theorem 4. When $ \frac{\alpha_{1}}{\alpha_{2}} > \frac{m_{1}}{r_{1}} $, then the equilibrium point $ E_{1}(x_1, 0) $ is a stable node if $ r_{2}-m_{2}-\frac{\alpha_{1}e_{1}a_{1}m}{\alpha_{2}}+e_{1}\eta_{1}+e_{2}\eta_{2} < 0 $, and the equilibrium point $ E_{1}(x_1, 0) $ is a saddle if $ r_{2}-m_{2}-\frac{\alpha_{1}e_{1}a_{1}m}{\alpha_{2}}+e_{1}\eta_{1}+e_{2}\eta_{2} > 0 $.

    Proof. The Jacobian matrix of $ E_1 $ is

    $ J_{E_{1}} = (m1α1r1α2α1a1mα2η1η20r2m2α1e1a1mα2+e1η1+e2η2)
    , $

    where

    $ \eta_{1} = \frac{a_{1}k(\alpha_{1}r_{1}-\alpha_{2}m_{1})}{r_{1}}, \eta_{2} = \frac{a_{2}k\alpha_{2}(\alpha_{1}r_{1}-\alpha_{2}m_{1})}{\alpha_{1}r_{1}b-\alpha_{2}k(\alpha_{1}r_{1}-\alpha_{2}m_{1})}. $

    Then the eigenvalues of $ J_{E_{1}} $ are $ \lambda_{3} = m_{1}-\frac{\alpha_{1}r_{1}}{\alpha_{2}} $, $ \lambda_{4} = r_{2}-m_{2}-\frac{\alpha_{1}e_{1}a_{1}m}{\alpha_{2}}+e_{1}\eta_{1}+e_{2}\eta_{2} $. It is obvious that $ \lambda_{3} < 0 $, hence the equilibrium point $ E_{1}(x_1, 0) $ is a stable node if $ \lambda_{4} < 0 $, and the equilibrium point $ E_{1} $ is a saddle if $ \lambda_{4} > 0 $.

    Theorem 5. The conditions for the existence of the internal equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $ are shown in Theorem 2. In this paper, it has been proved that there is and only one internal equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $. If $ Tr(J_{E_1^*}(x_{1}^*, y_{1}^*)) < 0 $, then the internal equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $ is locally asymptotically stable; otherwise, the internal equilibrium point $ E_1^*(x_{1}^*, y_{1}^*) $ is unstable.

    Proof. The Jacobian matrix of the internal equilibrium point $ E_1^* $ is

    $ J_{E_1^*} = (r1α1α2(12x1k)a1α1y1a2α2by1(b+α2x1)2m1a1α1(x1mα2)a2α2x1b+α2x1e1a1α1y1+a2e2α2by1(b+α2x1)20)
    . $

    The determinant of the matrix $ J_{E_1^*} $ is

    $ Det(JE1)=(e1a1α1y1+a2e2α2by1(b+α2x1)2)(a1α1(x1mα2)+a2α2x1b+α2x1).
    $

    The trace of the matrix $ J_{E_1^*} $ is

    $ Tr(JE1)=r1α1α2(12x1k)a1α1y1a2α2by1(b+α2x1)2m1.
    $

    It is obvious that all terms of the product of determinants are positive, hence, we can judge that the determinants of the Jacobian matrix of $ E_1^* $ is always positive. Thus, when $ Tr(J_{E_1^*}) < 0 $, the internal equilibrium point $ E_1^* $ has two negative eigenvalues, which means that the internal equilibrium point $ E_1^* $ is locally asymptotically stable. Otherwise, both characteristic roots are positive, and the internal equilibrium point $ E_1^* $ is unstable.

    In order to probe some possible coexistence evolution modes of cyanobacteria and potential grazers during the outbreak of cyanobacteria bloom, and clarify how the sheltering effect of colonial cyanobacteria affects the coexistence mode of cyanobacteria and potential grazers, we choose $ m $ as the control parameter to study the bifurcation dynamic behavior of the model (2.1), mainly including Transcritical bifurcation and Hopf bifurcation. If the model (2.1) can have Transcritical bifurcation, which implies that the cyanobacteria and potential grazers can coexist in a certain range at least, and there will be no extinction of at least one population with the advance of time. If the model (2.1) can have Hopf bifurcation, which shows that cyanobacteria and potential grazers can coexist in a periodic oscillation mode in a certain range, and the two populations can continue to survive with the advance of time. Therefore, bifurcation dynamics of the model (2.1) has important significance in population ecology and is worthy of our in-depth study.

    The existence of Transcritical bifurcation is studied at the boundary equilibrium point $ E_{1}(x_1, 0) $, we should ensure that $ \frac{\alpha_{1}}{\alpha_{2}} > \frac{m_{1}}{r_{1}} $ is satisfied. Thus, we state the following theorem.

    Theorem 6. The model (2.1) can undergo a Transcritical bifurcation at the boundary equilibrium poin $ E_{1}(x_1, 0) $ when $ m \equiv m_{TC} = \frac{\alpha_{2}}{\alpha_{1}e_{1}a_{1}}(r_{2}-m_{2}+e_{1}\eta_{1}+e_{2}\eta_{2}) $.

    Proof. Now, we use Sotomayor's Theorem to verify the transversality condition for the occurrence of Transcritical bifurcation at $ m \equiv m_{TC} = \frac{\alpha_{2}}{\alpha_{1}e_{1}a_{1}}(r_{2}-m_{2}+e_{1}\eta_{1}+e_{2}\eta_{2}) $, it can be easily verified that $ Det(J_{E_{1}}) = 0 $. Since $ \lambda_{3}\lambda_{4} = Det(J_{E_{1}}) $, which implies that $ J_{E_{1TC}} $ has a zero eigenvalue, that is to say $ \lambda_{4} = 0 $, thus the Jacobian matrix at $ E_{1} $ is

    $ J_{E_{1TC}} = (m1α1r1α2α1a1mα2η1η200)
    . $

    Let $ v $ and $ w $ be the two eigenvectors corresponding to the zero eigenvalue of the matrices $ J_{E_{1TC}} $ and $ J_{E_{1TC}}^T $ respectively, and they are given by

    $ v = (v11)
    , w = (01)
    , $

    where $ v_{1} = \frac{\alpha_{2}(r_{2}-m_{2}-e_{1}\eta_{2}+e_{2}\eta_{2})}{e_{1}(\alpha_{1}r_{1}-\alpha_{2}m_{1})} $.

    Furthermore, we can get

    $ F_{m}(E_{1}, m_{TC}) = (α1a1yα2α1e1a1yα2)
    _{(E_{1}, m_{TC})} = (00)
    , $
    $ DF_{m}(E_{1}, m_{TC})v = (0α1a1α20α1e1a1α2)
    (v11)
    _{(E_{1}, m_{TC})} = (α1a1α2α1e1a1α2)
    , $
    $ D^2F_{m}(E_{1}, m_{TC})(v, v) = (2F1x2v21+22F1xyv1v2+2F1y2v222F2x2v21+22F2xyv1v2+2F2y2v22)
    _{(E_{1}, m_{TC})} $
    $ = (0α2(r2m2e1η2+e2η2)e1(α1(r1+m1)m1)(a1e1α1+a2e2bα2(b+α2(α1r1α2m1)r1α1)2))
    , $

    clearly, they can satisfy the transversality conditions:

    $ w^T F_{m}(E_{1}, m_{TC}) = 0, $

    $ w^T [DF_{m}(E_{1}, m_{TC})v] = -\frac{\alpha_{1}e_{1}a_{1}}{\alpha_{2}} \neq 0, $

    $ w^T [D^2 F_{m}(E_{1}, m_{TC})(v, v)] = 2\frac{\alpha_{2}(r_{2}-m_{2}-e_{1}\eta_{2}+e_{2}\eta_{2})}{e_{1}(\alpha_{1}(r_{1}+m_{1})-m_{1})}(a_{1}e_{1}\alpha_{1}+\frac{a_{2}e_{2}b\alpha_{2}}{(b+\frac{\alpha_{2}(\alpha_{1}r_{1}-\alpha_{2}m_{1})}{r_{1}\alpha_{1}})^2}) \neq 0. $

    Therefore, when $ m \equiv m_{TC} = \frac{\alpha_{2}}{\alpha_{1}e_{1}a_{1}}(r_{2}-m_{2}+e_{1}\eta_{1}+e_{2}\eta_{2}) $, the model (2.1) can induce Transcritical bifurcation at the boundary equilibrium point $ E_{1} $.

    In the view of the proof of Theorem 5, we can know that Hopf bifurcation can occur at the positive equilibrium point $ E_1^* $. At the same time, the positive equilibrium point $ E_1^* $ can change its stability when the value of $ m $ passes through the critical magnitude. Thus, we summarize our findings in the following theorem.

    Theorem 7. Based on Theorem 5, the model (2.1) can induce a Hopf bifurcation when $ m = m_{HP} $.

    Proof. The characteristic equation of matrix $ J_{E_1^*} $ is $ \lambda^2-Tr(J_{E_1^*})\lambda + Det(J_{E_1^*}) = 0 $, and three cross-sectional conditions of Hopf bifurcation occurs, i.e.,

    1) $ Tr(J_{E_1^*})_{m = m_{HP}} = 0 $,

    2) $ Det(J_{E_1^*})_{m = m_{HP}} > 0 $,

    3) $ \frac{ {dTr(J_{E_1^*})}}{ {d\alpha}}\mid_{m = m_{HP}} = 0 $.

    Previously, we can show that $ Det(J_{E_1^*})_{m = m_{HP}} > 0 $. And when $ m = m_{HP} $, we have $ Tr(J_{E_1^*})_{m = m_{HP}} = 0 $. Thus we have $ \frac{ {dTr(J_{E_1^*})}}{ {d\alpha}}\mid_{m = m_{HP}} = 0, $ where

    $ Tr(JE1)m=mHP=r1α1α2(12x1k)a1α1y1a2α2by1(b+α2x1)2m1=0.
    $

    In conclusion, the model (2.1) can induce a Hopf bifurcation when $ m = m_{HP} $. Next, we will discuss the stability of the limit cycle by computing the first Lyapunov number. Translate the equilibrium point $ E_1^* $ to the origin by using the following transformation:

    $ x=xmx1,y=ymy1.
    $

    Hence, we obtain

    $ {˙xm=a10xm+a01ym+a20x2m+a11xmym+a02y2m+a30x3m+a21x2mym+a12xmy2m+a03y3m+Q(xm,ym),˙ym=b10xm+b01ym+b20x2m+b11xmym+b02y2m+b30x3m+b21x2mym+b12xmy2m+b03y3m+P(xm,ym),
    $

    where

    $ \ \ \ \ a_{10} = a_{1}\alpha_{1}y_{1}^*+\frac{a_{2}\alpha_{2}y_{1}^*}{b-x_{1}^*\alpha_{2}}+\frac{r_{1}\alpha_{1}k+2r_{1}\alpha_{1}x_{1}^*}{\alpha_{2}k}+\frac{a_{2}\alpha_2^2x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^2}-m_{1}, \\ \ \ \ \ a_{20} = -\frac{r_{1}\alpha_{1}}{\alpha_{2}k}-\frac{a_{2}\alpha_2^2y_{1}^*}{(b-x_{1}^*\alpha_{2})^2}-\frac{a_{2}\alpha_2^3x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^3}, \ \ \ \ \ \ \ \ a_{30} = \frac{a_{2}\alpha_2^3y_{1}^*}{(b-x_{1}^*\alpha_{2})^3}+\frac{a_{2}\alpha_2^4x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^4}, \\ \ \ \ \ a_{21} = \frac{a_{2}\alpha_2^2}{(b-x_{1}^*\alpha_{2})^2}+\frac{a_{2}\alpha_2^3x_{1}^*}{(b-x_{1}^*\alpha_{2})^3}, \ \ \ \ \ \ \ \ a_{11} = -a_{1}\alpha_{1}-\frac{a_{2}\alpha_{2}}{b-x_{1}^*\alpha_{2}}-\frac{a_{2}\alpha_2^2x_{1}^*}{(b-x_{1}^*\alpha_{2})^2}, \\ \ \ \ \ a_{01} = \frac{\alpha_{1}a_{1}m}{\alpha_{2}}+a_{1}x_{1}^*y_{1}^*+\frac{a_{2}\alpha_{2}x_{1}^*}{b-x_{1}^*\alpha_{2}}, \ \ \ \ \ \ \ \ b_{10} = -y^*a_{1}e_1\alpha_1-\frac{e_2a_{2}\alpha_{2}y_{1}^*}{b-x_{1}^*\alpha_{2}}-\frac{e_2a_{2}\alpha_2^2x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^2}, \\ \ \ \ \ b_{20} = \frac{e_2a_{2}\alpha_2^2y_{1}^*}{(b-x_{1}^*\alpha_{2})^2}+\frac{e_2a_{2}\alpha_2^3x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^3}, \ \ \ \ \ \ \ \ b_{30} = -\frac{e_2a_{2}\alpha_2^3y_{1}^*}{(b-x_{1}^*\alpha_{2})^3}-\frac{e_2a_{2}\alpha_2^4x_{1}^*y_{1}^*}{(b-x_{1}^*\alpha_{2})^4}, \\ \ \ \ \ b_{21} = -\frac{e_2a_{2}\alpha_2^2}{(b-x_{1}^*\alpha_{2})^2}-\frac{e_2a_{2}\alpha_2^3x_{1}^*}{(b-x_{1}^*\alpha_{2})^3}, \ \ \ \ \ \ \ \ b_{11} = \frac{e_2a_{2}\alpha_{2}}{b-x_{1}^*\alpha_{2}}+\frac{e_2a_{2}\alpha_2^2x_{1}^*}{(b-x_{1}^*\alpha_{2})^2}+a_{1}e_{1}\alpha_{1}, \\ \ \ \ \ b_{01} = -\frac{\alpha_{1}a_{1}e_{1}m}{\alpha_{2}}-m_2+r_2-a_{1}e_1\alpha_1x_{1}^*-\frac{e_2a_{2}\alpha_{2}x_{1}^*}{b-x_{1}^*\alpha_{2}}, \ \ \ \ \ \ \ \ a_{02} = a_{12} = a_{03} = b_{02} = b_{12} = b_{03} = 0. $

    And then $ Q(x_{m}, y_{m}) $, $ P(x_{m}, y_{m}) $ are power series in $ (x_{m}, y_{m}) $ with terms $ x_m^iy_m^j $, which can satisfy $ i+j \geq 4 $. Therefore, the first Lyapunov number to determine the stability of limit cycle is given by the formula

    $ l1=3π2a01Δ3/2{[a10b10(a211+a11b02+a02b11)+a10a01(b211+a20b11+a11b02)+b210(a11a02+2a02b02)2a10b10(b202a20a02)2a10a01(a220b20b02)a201(2a20b20+b11b20)+(a01b102a210)(b11b02a11a20)](a210+a01b10)[3(b10b03a01a30)+2a10(a21+b12)+(b10a12a01b21)]}. 
    $

    It is well known that the limit cycle generated by the Hopf bifurcation is unstable if $ l_{1} > 0 $ and stable if $ l_{1} < 0 $. Nevertheless, the expression of the first Lyapunov number is too complex to judge the positive and negative. Therefore, we will give some numerical solutions in Section 6.

    Under the research framework of theoretical derivation, we will not only verify the feasibility and effectiveness of the theoretical results, but also explore how clustering affects the dynamic relationship between cyanobacteria and potential grazers, thus the parameter $ m $ and $ \beta = \frac{\alpha_{1}}{\alpha_{2}} $ are selected as control variables for numerical simulation with $ r_1 = 0.56 $, $ r_2 = 0.2 $, $ k = 5 $, $ a_1 = 0.6 $, $ a_2 = 0.4 $, $ b = 0.15 $, $ m_1 = 0.18 $, $ m_2 = 0.28 $, $ e_1 = 0.6 $, $ e_2 = 0.28 $, it should be noted that all parameter values are estimated according to biological significance and theoretical conditions. Furthermore, $ \beta = \frac{\alpha_{1}}{\alpha_{2}} $ represents the proportion of colonial cyanobacteria and unicellular cyanobacteria, which can also indicate the outbreak degree of cyanobacteria bloom. If $ \beta = \frac{\alpha_{1}}{\alpha_{2}} > 1 $ indicates that the manifestation of cyanobacteria is mainly colonial cyanobacteria, which is not conducive to the capture of potential grazers. If $ \beta = \frac{\alpha_{1}}{\alpha_{2}} < 1 $ indicates that the manifestation of cyanobacteria is mainly unicellular cyanobacteria, which is conducive to the capture of potential grazers. Moreover, it should be pointed out that the predation behavior of potential grazers on cyanobacteria is essentially the implementation of ecological control technology.

    When we take $ \beta = 0.5385 $ with $ \alpha_{1} = 0.35 $ and $ \alpha_{2} = 0.65 $, which represents that the cyanobacteria in the whole ecosystem are in the bloom period, showing the phenomenon of aggregation. In other words, the current situation is not conducive to potential grazers to prey on cyanobacteria. Based on numerical simulation work, the bifurcation diagram of the model (2.1) with control parameter $ m $ has been shown in Figure 2(a). It is easy to find that the model (2.1) will experience a Hopf bifurcation and Transcritical bifurcation with the value of $ m $ changing within $ [0, 1.8] $. When the value of $ m $ gradually decreases from $ 1.8 $ to $ m_{TC} = 1.42282 $, the model $ (2.1) $ will occur a Transcritical bifurcation (Figure 3), the boundary equilibrium point $ E_1 $ will change from a stable state to an unstable state, which suggests that potential grazers will not always be close to extinction.

    Figure 2.  The evolution process of coexistence mode of cyanobacteria and potential grazers. The red line represents the internal equilibrium point $ x_1^* $ changing by $ m $, solid line represents the equilibrium point stable, dashed line represents the equilibrium point unstable and two vertical dashed lines represent the critical value of bifurcation. More detailed, Hp and TC are the critical values for the Hopf bifurcation, Transcritical bifurcation of the model, respectively. (a) $ m_{HP} = 0.2957018 $, $ m_{TC} = 1.42282 $; (b) $ m_{HP} = 0.1412414 $, $ m_{TC} = 1.74941 $.
    Figure 3.  Evolution process of steady-state coexistence mode with $ \alpha_{1} = 0.35 $ and $ \alpha_{2} = 0.65 $.

    And when the value of $ m $ is between $ (m_{HP}, m_{TC}) $, the model $ (2.1) $ has an asymptotically stable internal equilibrium point $ E_1^* $ (seeing Figure 3(b)), which indicates that cyanobacteria and potential grazers can coexist in a stable equilibrium state. If the value of $ m $ gradually decreases to $ m_{HP} = 0.2957018 $, the model $ (2.1) $ will occur Hopf bifurcation, the internal equilibrium point will change from a stable state to an unstable state, and a stable limit cycle will appear, the detailed numerical simulation results are shown in Figures 4 and 5. Furthermore, we can get that the Lyapunov exponent is $ -2.57010 \times 10^{7} $, which can further verify that the limit cycle is stable. Moreover, it should be explained that the dynamic characteristics of the model (2.1) have changed substantially, which means that the coexistence state of cyanobacteria and potential grazers changes from a stable equilibrium state to a stable periodic oscillation state.

    Figure 4.  Evolution process of periodic oscillation coexistence mode with $ \alpha_{1} = 0.35 $ and $ \alpha_{2} = 0.65 $. (a) We have stable periodic orbits bifurcate through Hopf bifurcation around $ E_1^*(0.48682, 0.16138) $ with $ m = m_{HP} = 0.2957018 $; (b) Local amplification of (a); (c) $ E_1^*(0.50822, 0.16493) $ is a locally asymptotically stable with $ m = 0.315 > m_{HP} $; (d) $ E_1^*(0.46966, 0.15840) $ is an unstable focus with $ m = 0.28 < m_{HP} $.
    Figure 5.  Dynamic evolution diagram of Hopf bifurcation with critical value $ m = m_{HP} = 0.2957018 $.

    When we take $ \beta = 1.083 $ with $ \alpha_{1} = 0.52 $ and $ \alpha_{2} = 0.48 $, which represents that the cyanobacteria in the whole ecosystem are still in the bloom period, but the current situation is quite disadvantageous to potential grazers to prey on cyanobacteria because unicellular cyanobacteria account for $ 0.48 $. It is easy to discover from Figure 2(b) that the model $ (2.1) $ has almost the same dynamic evolution process. When the value of $ m $ gradually decreases from $ 2 $ to $ m_{TC} = 1.74941 $, a Transcritical bifurcation will occur (Figure 6). When the value of $ m $ continues to decrease to $ m_{HP} = 0.125 $, the instability of the equilibrium point can produce a stable limit cycle, the detailed dynamic results are shown in Figures 7 and 8.

    Figure 6.  Evolution process of steady-state coexistence mode with $ \alpha_{1} = 0.52 $ and $ \alpha_{2} = 0.48 $.
    Figure 7.  Evolution process of periodic oscillation coexistence mode with $ \alpha_{1} = 0.35 $ and $ \alpha_{2} = 0.65 $. (a) We have stable periodic orbits bifurcate through Hopf bifurcation around $ E_1^*(0.38961, 0.58723) $ with $ m = m_{HP} = 0.14124142 $; (b) Local amplification of (a); (c) $ E_1^*(0.41055, 0.60801) $ is locally asymptotically stable with $ m = 0.155 > m_{HP} $; (d) $ E_1^*(0.36531, 0.56250) $ is unstable focus with $ m = 0.125 < m_{HP} $.
    Figure 8.  Dynamic evolution diagram of Hopf bifurcation with critical value $ m = m_{HP} = 0.1412414 $.

    At the same time, the Lyapunov exponent is $ -5.15158 \times 10^{14} $, which can show the stability of the limit cycle. From the essence of system dynamics evolution, whether the value of $ \beta $ is greater than $ 1 $ or less than $ 1 $, the model $ (2.1) $ has exactly the same dynamic evolution characteristics, that is, Transcritical bifurcation and Hopf bifurcation. That is to say, the ratio of colonial cyanobacteria to unicellular cyanobacteria does not affect the evolution characteristics of the dynamic state of the model $ (2.1) $. Moreover, it is also worth emphasizing that the value of parameter $ m $ seriously affects the dynamic relationship and coexistence mode between cyanobacteria and potential grazers.

    Based on the comparative analysis of numerical simulation results, a Transcritical bifurcation will occur at $ m_{TC} = 1.42282 $ when the value of $ \beta $ is less than $ 1 $, however a Transcritical bifurcation will occur at $ m_{TC} = 1.74941 $ when the value of $ \beta $ is greater than $ 1 $. This result means that during the outbreak of cyanobacteria bloom, the earlier the clustered cyanobacteria population is crushed by physical spray treatment technology, the more opportunities there are to change the dynamic relationship between cyanobacteria and potential grazers, and provide an entry point for the coexistence of cyanobacteria and potential grazers. However, the critical value of inducing Hopf bifurcation is just the opposite. When the value of $ \beta $ is less than $ 1 $, Hopf bifurcation can be induced at $ m = 0.2957018 $, which suggests that potential grazers actively prey on unicellular cyanobacteria to obtain sufficient food, which can lead to a large increase in the biomass of potential grazers. This increased predation ability can lead to essential changes in the dynamic relationship between them, and then play a role in controlling the large-scale outbreak of cyanobacteria bloom. However, when the value of $ \beta $ is greater than $ 1 $, Hopf bifurcation can be induced at $ m = 0.14124142 $, which means that the number of unicellular cyanobacteria can not basically maintain the growth needs of potential grazers, and a large number of colonial cyanobacteria need to be prey, which will cause colonial cyanobacteria to be gradually broken. This result also leads to a stable equilibrium between cyanobacteria and potential grazers in the range of $ m $. Therefore, from the perspective of controlling cyanobacteria bloom by physical spray treatment technology, the colonial cyanobacteria turns into unicellular cyanobacteria after the cyanobacteria bloom is treated by physical spray, which is conducive to potential grazers to prey, and then lead to the coexistence mode of potential grazers and cyanobacteria in a stable state.

    Within the framework of physical and ecological integrated control and management of cyanobacteria bloom, according to the growth characteristics and comprehensive management characteristics of wild cyanobacteria population, we mainly focus on the impact of physical spraying treatment technology on the structural characteristics of cyanobacteria population and the predation effect of potential grazers, and give some ecological modeling assumptions. On this basis, a aquatic ecological dynamic new model is proposed to describe the dynamic relationship between cyanobacteria and potential grazers. At the same time, it is worth emphasizing from Figure 1 that the biggest advantage of this new model is that it integrates three kinds of dynamic action modes, of which the sub model $ (2.1) $ is the most distinctive and has not been studied by predecessors. Furthermore, the biggest advantage of the model $ (2.1) $ is that it can perfectly describe the dynamic relationship between cyanobacteria and potential grazers under the physical and ecological integrated control, most commendably, it can control cyanobacteria bloom from the perspective of minimum value, and the consumption of potential predators is limited.

    Underlying the analysis of mathematical theory, we first explore the existence and stability of the equilibrium point of the model $ (2.1) $, and give some critical conditions of some critical parameters. Then, we obtain some key threshold conditions to ensure that the model $ (2.1) $ has specific dynamic behaviors, such as Transcritical bifurcation and Hopf bifurcation. These theoretical results are the basis for our numerical simulation analysis, which provides theoretical support for in-depth exploration of the dynamic relationship between cyanobacteria and potential grazers.

    Based on the results of numerical analysis, the dynamic behavior type and dynamic evolution process of the model $ (2.1) $ are clarified, focusing on the Transcritical bifurcation and Hopf bifurcation. It is easy to find from Figure 2 that when cyanobacteria bloom break out and gather into clusters, we implement a physical and ecological comprehensive control strategy for cyanobacterial bloom, the dynamic relationship between cyanobacteria and potential grazers suddenly changes by Transcritical bifurcation, which can form coexistence mode in a stable equilibrium state, the detailed dynamic simulation results are shown in Figures 3 and 6. After then, the biomass of cyanobacteria and potential grazers gradually decrease with the reduction of agglomeration, finally, a stable periodic oscillation coexistence mode is formed to maintain a good cycle of ecological algae control strategy, these dynamic simulation results are shown in Figures 4, 5, 7 and 8. Furthermore, it is more meaningful to obtain from Figures 3 and 6 that Transcritical bifurcation represents evolution process of steady-state coexistence mode, and it is the most ecologically significant result form Figures 4 and 7 that Hopf bifurcation represents evolution process of periodic oscillation coexistence mode. In addition, it should be emphasized that the most important role of physical spraying treatment technology is to break up clumps of cyanobacteria in the process of controlling cyanobacteria bloom, which will not change the dynamic essential characteristics of cyanobacteria and potential grazers represented by the model $ (2.1) $.

    Based on theoretical and numerical research results, the main results of this paper are as follows: $ 1) $ The population density of cyanobacteria decreases gradually as the value of key parameter $ m $ decreases, that is to say, the implementation of physical spraying treatment technology can effectively control the population density of cyanobacteria. $ 2) $ Transcritical bifurcation can induce steate-state coexistence mode between cyanobacteria and potential grazers, and Hopf bifurcation can force cyanobacteria and potential grazers to form periodic oscillation coexistence mode. $ 3) $ The evolutionary characteristics of bifurcation dynamics suggest that physical and ecological integrated control technique can effectively inhibit the growth of cyanobacteria, and promote the formation of a stable coexistence mode between cyanobacteria and potential grazers within a controllable range.

    In general, some results have been obtained by proposing a new aquatic ecological model, however, the work of this paper still has some limitations, such as: $ 1) $ the characteristics of toxicity released by cyanobacteria population were not considered in the modeling, $ 2) $ cyanobacteria migration behavior is not mentioned, $ 3) $ modeling assumptions make that the application of the model $ (2.1) $ has certain limitations. Thus there is still much work to be further studied, for example, how does the release of cyanobacteria toxins affects the predatory ability of potential grazers during physical spray treatment, effect of cyanobacteria migration behavior on algae control effect of physical spraying technology, etc. In a word, it is hoped that the research results of this paper will provide some theoretical support for the application of physical and ecological comprehensive management of cyanobacteria bloom.

    This work was supported by the key International Cooperation Projects of China (Grant No: 2018YFE0103700), the National Natural Science Foundation of China (Grant No: 31570364, 61871293, 41876124, 61901303), by the Zhejiang Provincial Natural Science Foundation of China (Grant No: 2Q18C030002), by key research and development projects of Zhejiang Province (Grant No: 2021C03166), by the Science and Technology program of Cangnan, China (Grant No: 2018ZG29), by Science and Technology Major Program of Wenzhou, China (Grant No: 2018ZG002), by the Science and Technology Program of Wenzhou, China (Grant No.X20210097).

    The authors declare there is no conflict of interest.

    [1] Feddema J, Oleson KW, Bonan G, et al. (2005) Importance of land cover change in simulating future climates. Science 310: 1674-1678. doi: 10.1126/science.1118160
    [2] Houghton RA, Hackler JL, Lawrence KT (1999) The U.S. carbon budget: contributions from land-use Change. Science 285: 574-578.
    [3] Caspersen JP, Pacala SW, Jenkins JC, et al. (2000) Contributions of land-use history to carbon accumulation in U.S. forests. Science 290: 1148-1151.
    [4] Vitousek PM, Mooney HA, Lubchenco J, et al. (1997) Human domination of earth's ecosystems. Science 277: 494-499. doi: 10.1126/science.277.5325.494
    [5] Foley JA, Defries R, Asner GP, et al. (2005) Global consequences of land use. Science 309: 570-574. doi: 10.1126/science.1111772
    [6] Pan Y, Birdsey RA, Fang J, et al. (2011) A large and persistent carbon sink in the world's forests. Science 333: 988-993. doi: 10.1126/science.1201609
    [7] U.S. Environmental Protection Agency (2014) Inventory of U.S. greenhouse gas emissions and sinks: 1990 - 2012. Available from: http://www.epa.gov/climatechange/emissions/usinventoryreport.html
    [8] Sleeter BM, Wilson TS, Soulard CE, et al. (2011) Estimation of late twentieth century land-cover change in California. Environ Monit Assess 173: 251-266. doi: 10.1007/s10661-010-1385-8
    [9] Liu J, Vogelmann JE, Zhu Z, et al. (2011) Estimating California ecosystem carbon change using process model and land cover disturbance data: 1951-2000. Ecol Model 222: 2333-2341. doi: 10.1016/j.ecolmodel.2011.03.042
    [10] Winrock International, PIER Energy-Related Environmental Research. California Energy Commission, 500-04-068F. 2014. Availalble from: http://www.energy.ca.gov/pier/project_reports/500-04-068.html
    [11] Birdsey RA, Lewis GM (2003) Carbon in U.S. forests and wood products, 1987-1997: state-by-state estimates. USDA Forest Service, Northeastern Research Station. Gen Tech Rep NE-310: 42.
    [12] State of California, Department of Finance, Report P-1 (Total Population): State and county population projections, 2010-2060. Sacramento CA, December 2014, Available from: http://www.dof.ca.gov/research/demographic/reports/projections/P-1/.
    [13] Nakicenovic N, Swart R, (eds) (2000) Special Report on Emission Scenarios: A special report of Working Group III of the Intergovernmental Panel on Climate Change. Cambridge University Press, UK.
    [14] Millenium Ecosystem Assessment (2005) Ecosystems and Human Well-being: Scenarios, Volume 2. Carpenter SR, Pingali PL, Bennett EM et al., (eds), Island Press, 1718 Connecticut Avenue, Suite 300, NW, Washington, DC 20009.
    [15] Moss RH, Edmonds JA, Hibbard KA, et al. (2010) The next generation of scenarios for climate change research and assessment. Nature 463: 747-756. doi: 10.1038/nature08823
    [16] Sleeter BM, Sohl TL, Bouchard MA, et al. (2012) Scenarios of land use and land cover change in the conterminous United States: Utilizing the special report on emission scenarios at ecoregional scales. Global Environ Chang 22: 896-914. doi: 10.1016/j.gloenvcha.2012.03.008
    [17] Radeloff VC, Nelson E, Plantiga AJ, et al. (2012) Economic-based projections of future land use in the conterminous United States under alternative policy scenarios. Ecol Appl 22: 1036-1049. doi: 10.1890/11-0306.1
    [18] van Vuuren DP, Lucas PL, Hilderink H (2007) Downscaling drivers of global environmental change: Enabling use of global SRES scenarios at the national and grid levels. Global Environ Chang 17: 114-130. doi: 10.1016/j.gloenvcha.2006.04.004
    [19] Verburg PH, Schulp CJE, Witte N, et al. (2006) Downscaling of land use change scenarios to assess the dynamics of European landscapes. Agr Ecosyst Environ 114: 39-56. doi: 10.1016/j.agee.2005.11.024
    [20] Wear D (2011) Forecasts of county-level land uses under three future scenarios: a technical document supporting the Forest Service 2010 RPA Assessment. Department of Agriculture Forest Service, Southern Research Station. Gen Tech Rep SRS-141: 41.
    [21] Zhu Z, et al, (ed) (2010) A method for assessing carbon stocks, carbon sequestration, and greenhouse-gas fluxes in ecosystems of the United States under present conditions and future scenarios. U.S. Geological Survey Scientific Investigations Report 2010-5233, Reston, VA, Available from: http://pubs.usgs.gov/sir/2010/5233/
    [22] IPCC (2007) Climate Change 2007: Synthesis Report. Contribution of Working Groups I, II, and III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Geneva, Switzerland: IPCC. 104.
    [23] Foley JA, Prentice C, Ramankutty N, et al. (1996) An integrated biosphere model of land surface processes, terrestrial carbon balance, and vegetation dynamics. Global Biogeochem Cy 10: 603-628. doi: 10.1029/96GB02692
    [24] Baker WL (1989) A review of models of landscape change. Landscape Ecol 2: 11-133.
    [25] Daniel CJ, Frid L (2012) Predicting Landscape Vegetation Dynamics Using State-and-Transition Simulation Models. In: Kerns BK, Shlisky AJ, Daniel CJ, (eds). U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station, General Technical Report PNW-GTR-869; 5-22.
    [26] APEX Resource Management Solutions,ST-Sim State and Transition Simulator modeling software. 2014. Available from: http://www.apexrms.com/stsm, version 2.4.6.
    [27] IPCC (2006) 2006 IPCC Guidelines for National Greenhouse Gas Inventories, Volume 4: agriculture, forestry and other land use. Paustian K, Ravindranath N.H., Amstel A (eds), Available from: http://www.ipcc-nggip.iges.or.jp/public/2006gl/vol4.html
    [28] Kucharik CJ, Foley JA, Delire C, et al. (2000) Testing the performance of a dynamic global ecosystem model: Water balance, carbon balance, and vegetation structure. Global Biogeochem Cy 14: 795-825. doi: 10.1029/1999GB001138
    [29] Farquhar GD, von Caemmerer S, Berry JA (1980) A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149: 78-90. doi: 10.1007/BF00386231
    [30] Ball JT, Woodrow I, Berry J (1987) A Model Predicting Stomatal Conductance and its Contribution to the Control of Photosynthesis under Different Environmental Conditions. In: Biggins J, editor. Progress in Photosynthesis Research: Springer Netherlands. 221-224.
    [31] Botta A, Viovy N, Ciais P, et al. (2000) A global prognostic scheme of leaf onset using satellite data. Global Change Biol 6: 709-725. doi: 10.1046/j.1365-2486.2000.00362.x
    [32] Parton WJ, Scurlock JMO, Ojima DS, et al. (1993) Observations and modeling of biomass and soil organic matter dynamics for the grassland biome worldwide. Global Biogeochem Cy 7: 785-809. doi: 10.1029/93GB02042
    [33] Verberne E, Hassink J, De Willigen P, et al. (1990) Modelling organic matter dynamics in different soils. Neth J Agr Sci 38: 221-238.
    [34] Parton WJ, Schimel DS, Cole C, et al. (1987) Division s-3-soil microbiology and biochemistry. Soil Sci Soc Am J 51: 1173-1179. doi: 10.2136/sssaj1987.03615995005100050015x
    [35] Liu J, Price DT, Chen JM (2005) Nitrogen controls on ecosystem carbon sequestration: a model implementation and application to Saskatchewan, Canada. Ecol Model 186: 178-195. doi: 10.1016/j.ecolmodel.2005.01.036
    [36] Zhu Q, Liu J, Peng C, et al. (2014) Modelling methane emissions from natural wetlands by development and application of the TRIPLEX-GHG model. Geosci Model Dev 7: 981-999. doi: 10.5194/gmd-7-981-2014
    [37] Anderson JR, Hardy EE, Roach JT, et al. (1976) A Land Use And Land Cover Classification System For Use With Remote Sensor Data. U.S. Geological Survey Professional Paper 964.
    [38] Loveland TR, Sohl TL, Stehman SV, et al. (2002) A strategy for estimating the rates of recent land cover changes in the United States. Photogramm Eng Rem S 68: 1091-1099.
    [39] Sleeter BM, Sohl TL, Loveland TR, et al. (2013) Land-cover change in the conterminous United States from 1973 to 2000. Global Environ Chang 23: 733-748. doi: 10.1016/j.gloenvcha.2013.03.006
    [40] Omernik JM (1987) Ecoregions of the Conterminous United States. Ann Assoc Am Geogr 77: 118-125. doi: 10.1111/j.1467-8306.1987.tb00149.x
    [41] Kim SJ, Flato GM, Boer GJ (2003) A coupled climate model simulation of the Last Glacial Maximum, Part 2: approach to equilibrium. Clim Dynam 20: 635-661.
    [42] Kim SJ, Flato GM, Boer GJ, et al. (2002) A coupled climate model simulation of the Last Glacial Maximum, Part 1: transient multi-decadal response. Clim Dynam 19: 515-537. doi: 10.1007/s00382-002-0243-y
    [43] Flato GM, Boer GJ (2001) Warming symmetry in climate change simulations. Geophys Res Lett 28: 195-198. doi: 10.1029/2000GL012121
    [44] Price DT, McKenney DW, Papadopol P, et al. (2004) High resolution future scenario climate data for North America; Proceedings of the American Meteorological Society 26th Conference on Agricultural and Forest Meteorology, Vancouver, BC, Canada, 23-26 August 2004. 13. http://ams.confex.com/ams/pdfpapers/78202.pdf
    [45] Cohen WB, Spies TA (1992) Estimating structural attributes of Douglas Fir/Western Hemlock forest stands from Landsat and SPOT imagery. Remote Sens Environ 41: 1-17. doi: 10.1016/0034-4257(92)90056-P
    [46] Kellndorfer J, Walker W, LaPoint E, et al. (2012) NACP Aboveground Biomass and Carbon Baseline Data (NBCD 2000), U.S.A., 2000. Available on-line at http://daac.ornl.gov: ORNL DAAC Oak Ridge, Tennessee, U.S.A.
    [47] Sleeter R, Acevedo W, Soulard CE, et al. (2015) Methods used to parameterize the spatially explicit components of a state and transition simulation model. AIMS Environ Sci [submitted].
    [48] Wilson BT, Woodall CW, Griffith DM (2013) Imputing forest carbon stock estimates from inventory plots to a nationally continuous coverage. Carbon Balance Manage 8: 1. doi: 10.1186/1750-0680-8-1
    [49] Zhu Z (ed) (2012) Baseline and Projected Future Carbon Storage and Greenhouse-Gas Fluxes in Ecosystems of the Western United States, U.S. Geological Survey Professional Paper 1797, 192 . Available from: http://pubs.usgs.gov/pp/1797/
    [50] Blackard J, Finco M, Helmer E, et al. (2008) Mapping U.S. forest biomass using nationwide forest inventory data and moderate resolution information. Remote Sens Environ 112: 1658-1677.
    [51] Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture. Soil Survey Geographic (SSURGO) Database. Available from: http://sdmdataaccess.nrcs.usda.gov/
    [52] Daly C, Halbleib M, Smith JI, et al. (2008) Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int J Climatol 28: 2031-2064. doi: 10.1002/joc.1688
    [53] Kurz WA, Apps MJ (1999) A 70-year retrospective analysis of carbon fluxes in the Canadian forest sector. Ecol Appl 9: 526-547. doi: 10.1890/1051-0761(1999)009[0526:AYRAOC]2.0.CO;2
    [54] Kurz WA, Dymond CC, White TM, et al. (2009) CBM-CFS3: A model of carbon-dynamics in forestry and land-use change implementing IPCC standards. Ecol Model 220: 480-504. doi: 10.1016/j.ecolmodel.2008.10.018
    [55] Bachelet D, Neilson RP, Lenihan JM, et al. (2001) Climate Change Effects on Vegetation Distribution and Carbon Budget in the United States. Ecosystems 4: 164-185. doi: 10.1007/s10021-001-0002-7
    [56] Westerling AL, Bryant BP (2007) Climate change and wildfire in California. Climatic Change 87: 231-249.
  • This article has been cited by:

    1. Mudasir Younis, Haroon Ahmad, Mahpeyker Ozturk, Deepak Singh, A novel approach to the convergence analysis of chaotic dynamics in fractional order Chua’s attractor model employing fixed points, 2025, 110, 11100168, 363, 10.1016/j.aej.2024.10.001
  • Reader Comments
  • © 2015 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(8984) PDF downloads(1574) Cited by(22)

Figures and Tables

Figures(12)  /  Tables(7)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog