EFFECT OF SEASONAL CHANGING TEMPERATURE ON THE GROWTH OF PHYTOPLANKTON

. An non-autonomous nutrient-phytoplankton interacting model incorporating the eﬀect of time-varying temperature is established. The impacts of temperature on metabolism of phytoplankton such as nutrient uptake, death rate, and nutrient releasing from particulate nutrient are investigated. The ecological reproductive index is formulated to present a threshold criteria and to characterize the dynamics of phytoplankton. The positive invariance, dis- sipativity, and the existence and stability of boundary and positive periodic solution are established. The analyses rely on the comparison principle, the coincidence degree theory and Lyapunov direct method. The eﬀect of seasonal temperature and daily temperature on phytoplankton biomass are simulated numerically. Numerical simulation shows that the phytoplankton biomass is very robust to the variation of water temperature. The dynamics of the model and model predictions agree with the experimental data. Our model and analysis provide a possible explanation of triggering mechanism of phytoplankton blooms.


1.
Introduction. Phytoplankton constitutes the bottom level of the aquatic food web and plays a fundamental role in the aquatic ecosystem. Phytoplankton bloom (PB) can assert severe negative consequences on ecosystems and social economics. In the past decades, there have been numerous field or experimental and also some mathematical modeling studies on the impact of some specific factors contributing to the growth of phytoplankton. Many researchers focused on the role of nutrients played on PB. Many studies show that eutrophication is the main reason leading to PB [9,18]. Johnson et al. [24] and Phillips et al. [30] stressed the importance of sediment phosphorus release in the restoration of very shallow lakes. Shukla et al. [32] modeled and analyzed PB in a lake caused by discharge of nutrients. Those studies can help to understand the contribution to PB from the different specific points. Nevertheless, most of the studies ignored the impact of temperature, which is thought to be a possible key triggering factor of PB.
There have been some studies emphasizing on the effect of constant temperature on the growth of phytoplankton [8,40]. Especially, Chen et al [5] presented a sound frame work on the effect of constant temperature on the growth of phytoplankton. They defined a helpful index called basic reproductive index to fully characterize the relationship between the dynamic evolution of phytoplankton and temperature.
The model considering the effect of constant temperature on the growth of phytoplankton is unsuitable for modeling PB that recurs seasonally. Then some aquatic ecologists have been shifting their attentions to seasonal variations of phytoplankton composition and has been attracted by concepts of periodicity, succession and response to environmental changes [33,38]. Some field or experimental studies (such as our field study in Fig. 2 and 3) suggested that the dominant variations in both phytoplankton composition and environmental variables are seasonal over a consistent annal cycle in lake and strongly related to the annual cycles of temperature [7,21,39]. Our goal of this paper is to develop modeling approaches for understanding the effect of seasonal changing temperature on the growth of phytoplankton and model the complex dynamics of seasonally recurring phytoplankton blooms.
The temperature plays important roles in many aspects during the process of phytoplankton metabolism. The impacts of temperature on the metabolism of phytoplankton are traditionally described by the Q 10 rule [11,14,42], which asserts that a change of the temperature by 10 o C will multiply the rate at mean temperature by a factor Q 10 , that is to say, the Q 10 temperature coefficient is a measure of the rate of change of a biological or chemical system as a consequence of increasing the temperature by 10 o C. The nonlinear nature of biological responses to temperature are indeed topics of high priority [18]. Aquatic ecologists focused on the relationship between temperature and the growth rate of phytoplankton [2,31,40], and the nutrients releasing rate from sediments [26], and the death rate of phytoplankton [8,25].
Located in the southern part of the Yangtze River delta, Lake Tai is the third largest freshwater lake in China with a surface area of about 2,338 square kilometers (about 902 square miles), an average depth of 2 meters, and maximum depth 3 meters. Lake Tai has been used as the important freshwater sources in the yangtze river delta on edge of Chinese eastern coast. Its water quality and environmental conditions have an important impact on the development of this area. In recent years, with the fast development of industry and agriculture process, waste water and domestic sewage without thorough treatment were discharged into Lake Tai Basin, then those pollutions gradually worsen the water quality and assert considerably negative impact on the ecological environment system of Lake Tai Basin. The water pollution of cyanobacteria blooms has become a severe problem in this area. At present, most of the studies were focused on Wuxi Meiliang bay of Lake Tai and the pollution in other typical areas around Lake Tai were less explored.
From March 2008 to April 2009, we have made an detailed field investigation of the water qualities and temperature variations around Lake Tai. The study was conducted in 12 selected monitor stations (i.e., MS) set in the typical areas around Lake Tai (Fig 1). We sampled the water, collected data, monitored and assessed  Table 1. the water quality once a month. The specimen water was sampled from 0.5m deep below the surface of water. The water temperature (W T ) was field determined simultaneously on the sampling spot by water thermometer (instrument: Waterchecker 2010, HACH). Then the water samples were taken back to the laboratory immediately in seal evades the light preservation and the other index of performance were mensurated and analyzed. Particularly, Chlorophyll-a is determined by the ultraviolet spectrophotometry (instrument: UV2450 Spectrophotometer, Shimadzu) and the total soluble reactive phosphorus (PO 4 -P) is monitored by Ion Chromatography method (instrument: DIONEX ICS-90). The annual mean values of PO 4 -P (P in1 ) are listed in Table 3. The density of phytoplankton is expressed by the density of Chlorophyll-a. The time series and observed patterns of water temperature and Chlorophyll-a in the targeted monitor stations T i (i = 1, 2, .., 12) are illustrated in Fig. 2 and 3, respectively. Our experiments, long-term observations, and our compilation of collected data suggest obviously seasonal patterns of the water temperature and Chlorophylla. The water temperature and the Chlorophyll-a increase first, then decrease, and their maxima and minima occur around August and January, respectively. It is well known that, since 1990s, the PB or algae bloom occurs annually in Lake Tai  [9]. The annual cyclic dynamics of phytoplankton in shallow lake is far from being well understood.
The remainder of this paper is organized as follows. In section 2, we propose a time-varying nutrient-phytoplankton dynamic model including the time-varying effect of temperature on metabolism of phytoplankton. In section 3, an ecological reproductive index is proposed and the existence and stability of periodic solutions are proved by coincidence degree theory and Lyapunov direct method. Section 4 carries out numerical simulations to expound both the long term cyclic dynamics and the short term effect of daily temperature on the growth of phytoplankton. The filed application is then carried out for Lake Tai to validate our model. The paper ends with conclusion and some interesting discussion.
2. Model formulation. In this section, following the clues and framework in [5], we formulate a dynamical model governed by non-autonomous ordinary differential equations for a lake with inflow and outflow of nutrient and study the impact of time-varying temperature on the dynamic evolution of both the phytoplankton and the nutrient. We assume that the lake is a physically homogeneous shallow lake or the well-mixed surface layer of a stratified water column, which thoroughly mixed at all times, so that there are no spatial gradient of concentrations.
The forced phytoplankton-nutrient compartment model consists of two state variables. Let A(t) and P (t) be the densities of phytoplankton and soluble mineral nutrient (say phosphorus) at time t in units of milligram of Chlorophyll-a per stere and in units of milligram of nutrient per stere, respectively. The system is assumed  to be governed by three basic processes: the phytoplankton nutrient uptake and growth, the death of phytoplankton, and the dilution of the nutrient.
The loss of phytoplankton due to natural mortality is assumed to be proportional to its density A and is determined by m, a parameter taking into account relatively constant grazing or predation by higher tropic levels. The loss due to the intraspecific competition is modeled by the quadratic mortality term bA 2 . The exchange of the nutrient in lake with the external water is denoted by a(P s − P ), where P s is the density of total soluble nutrient input. The uptake of phytoplankton for nutrient is assumed to be U p f (P )A with U p being the maximum nutrient uptake rate of phytoplankton and f (P ) being the nutrient-limited functional response. The growth of phytoplankton is assumed to be maximum growth rate µ opt limited by the nutrient uptake. Then the frame model is regulated by the following equations Next, we discuss the mathematical formulation of the components of (1), in particular, we will expound their dependence on the temperature.
Nutrient-limited functional response. The nutrient-limited functional response f (P ) usually satisfies the following properties [13,27] f (P ) ∈ C 2 , f (0) = 0, f (P ) > 0, lim p→∞ f (P ) = 1. (2) The most typical and common used is the so-called Hill function f n (P ) = P n H n P + P n , which incorporates the Holling family of functional responses [22] as special cases, where H P is the half-saturation constant for nutrient uptake. For example, if n = 1, then (3) reduces to the classical Holling type II functional response and depicts a decreasing rate of change of ingestion as food becomes more abundant such that the relationship is a concave downward curve. If n = 2, then (3) is the classic sigmoidal Holling type III functional response, which can be expressed in terms of Michaelis-Menten parameters and is concave downward only once the food surpasses a certain level (the inflection point). Motivated by the Holling family, a more general one can be extended by the following formulation to model the nutrient uptake of phytoplankton where r 1 > r 2 ≥ 0 and H i > 0 for i = 1, 2. The four parameters govern the shape of the response curve, no matter it is a concave downward curve or a sigmoidal curve. Temperature is a key environmental factor for the lake system and asserts essential impacts on the nutrient uptake, growth and death of phytoplankton, the nutrient releasing from particulate nutrients, and so on. Next, we shall deliberately discuss these metabolism progresses respectively. Temperature-dependent uptake and growth. The nutrient uptake and growth of phytoplankton are heavily affected by the variation of temperature. There are several model candidates depicting the effect of temperature such as linear model [1,34], Michaelis-Menten-Monod saturating model [34], optimal exponential model [3,40] and Arrhenius model [17,19]. In this study, we adopt the so-called cardinal temperature model with inflexion (CTMI) proposed by Rosso et al [31], which is validated for a large range of species. As we all known, phytoplankton grows much better under more suitable temperature condition and halts growth under extremely low or high temperature. Bernard [2] shows that the CTMI can accurately represent the growth of microalgae. The CTMI includes four key parameters and three of them are cardinal temperatures with a biological significance, which makes the model rather straightforward to calibrate and facilitates the estimation with experimental data. It predicts the growth rate for temperatures between T min and T max , that is, µ opt · g(T ), where with Here, g(T ) is the impact factor of water temperature modeling the direct effect of temperature on phytoplankton growth. The maximum growth rate µ opt occurs when T = T opt , where T opt is the optimal (most appropriate) temperature for the phytoplankton's production (growth). Fig. 4 illustrates the shape of CTMI. The CTMI depicts the relationship between growth rate and temperature and is more credible and acceptable [2,4,20]. Temperature-dependent mortality. Water temperature essentially affects the natural mortality of phytoplankton. Higher temperature may lead to higher death rate [25]. Ding [8] adopted an exponential model to address the death rate with respect to water temperature, which reads where m 0 is the minimum death rate of phytoplankton, K 1 is the respiration rate of phytoplankton, θ 1 is the temperature coefficient, and T ref is the reference water temperature.
Temperature-dependent nutrient release. We assume that there is an external source of nutrients (say phosphorus) flowing into the lake and the densities write as P in1 and P in2 , where P in1 denotes the soluble nutrient and P in2 represents the particulate nutrient. Denote by P b the particulate nutrient in bottom sediments. Moreover, the particulate nutrient in the water column and the nutrient in the sediments release soluble nutrient at rate γ. The releasing rate γ is affected by temperature and then it can be written as [26] , where K 2 and θ 2 are positive coefficients. Therefore, the total input soluble nutrient in the water column readsP s = aP in1 + γ(T )(P in2 + P b ) := aP s .
Time-varying temperature. In real world application, the environmental factors and resources vary or fluctuate in time. The environmental fluctuations are likely to have a significant effect on the dynamic evolution of aquatic ecosystems. In order to study the effect of time-varying temperature forcing phytoplankton, T must be time-dependent (say T (t)). T (t) can be fitted by some simple functions or their combination to illustrate the varying temperature such as the constant function, step function, quadratic function, and so on.
Consider the specific case that the temperature is assumed to be varying periodically, i.e., T (t) represents a periodic function. A convenient and commonly applied scheme is that of sinusoidal forcing [14,23,29] T (t) = T mean + α sin(Ωt + ϕ), where T mean represents the mean temperature and Ω governs the period. When the forcing is annual, and t is in units of months, then Ω = π/6. The parameter α controls the amplitude of forcing. It is obvious that, if α = 0, the model collapses to an unforced system. This case has been well explored by Chen et al [5], where the temperature is simplified to be constant or on an average, i.e., T = T mean . Chen et al [5] have proved that such a model can not admit any periodically oscillating dynamics. That is to say, the time-independent phytoplankton-nutrient model proposed in [5] can not catch or characterize the seasonal pattern of phytoplankton bloom. Although (5) is also a crude simplification, it serves as a first approximation for investigating how seasonal factors influences the model's dynamics.
To sum up, the objective phytoplankton-nutrient system takes the following form The description, unit, default value (when available), and reference resource for each parameter are summarized in Table 2. 3. Oscillatory dynamics of model (6). In this section, we explore the dynamics of (6) and present some results on the positive invariance, the existence and globally asymptotic stability of boundary and internal cyclic dynamics. In order to facilitate the following discussion, write which are assumed to be continuous and nonnegative.

Positive invariance and dissipativity. For a bounded continuous function
is positively invariant with respect to (6), where Proof. Let (A(t), P (t)) be the solution of (6) through (A(t 0 ), P (t 0 )) ∈ Γ. From the phytoplankton equation in (6), it follows that A standard comparison argument shows that which, together with the nutrient equation in (6), produces and hence Similarly, the nutrient equation in (6) yields whence Moreover, by the phytoplankton equation of system (6), we have which implies Therefore, Γ is positively invariant with respect to (6). The proof is complete. (6) is ultimately bounded with respect to Γ , i.e., (6) is dissipative.

Ecological reproductive index.
In the rest of Section 3, T (t) is assumed to be ω-periodic, i.e., there exists an ω > 0 such that T (t) = T (t + ω), and hence In order to characterize the population dynamics of phytoplankton, motivated by the computational formula of a threshold parameter for periodic compartmental epidemic models given by [37], here we introduce the ecological reproductive index R 0 for the periodic nutrient-phytoplankton model (6) . It is not difficult to show that (6) always admit a phytoplankton-free state E 0 = (0, P * (t)), where is the unique positive periodic solution of the following system Let E = (A, P ) T , then we rewrite the model (6) as where .
(16) By carrying out arguments similar to those in [10, Lemma 1] , one has where F (t) = µ opt g(t)f (P * (t)) and V (t) = m(t). The reproduction rate of phytoplankton at E 0 is F (t) and the removal rate of phytoplankton at E 0 is V (t). Then we get the following cooperative system Denote by Φ F (·)−V (·) (t) the monodromy matrix of (18), by ρ(Φ F (·)−V (·) (ω)) the spectral radius of Φ F (·)−V (·) (ω), and by Φ V (·) (t) the monodromy matrix of the linear ω-periodic differential system Hence, the monodromy matrix Assume that φ(s) is ω-periodic in s and denotes the initial distribution of phytoplankton individuals. Then F (s)φ(s) is the rate of the new population generated by initial fertile phytoplankton individuals being introduced at time s. Given t > s, then Y (t, s)F (s)φ(s) gives the distribution of those fertile individuals that were newly produced at time s and remain in the fertile compartments at time t [36]. It follows that is the distribution of the accumulative new individuals at time t produced by all those fertile individuals φ(s) introduced at time previous to t.
Let C ω denote the ordered Banach space of all ω-periodic functions from R to R and is equipped with the maximum norm · . Define a linear operator L on C ω by that is From [36] and [37], it follows that the ecological reproductive index R 0 is given by ρ(L), where ρ denotes the spectral radius of L.
In order to formulate the explicit expression of R 0 , we need to introduce the linear ω-periodic system with σ ∈ (0, ∞). Since F (t) is nonnegative and −V (t) is cooperative, it follows that ρ(Φ(ω, σ)) is continuous and nonincreasing in σ ∈ (0, ∞) and lim λ→∞ ρ(Φ(ω, σ)) < 1. It is not difficult to verify that (6) satisfies the assumptions (A 1 ) − (A 7 ) in [37]. Now, we list below two main conclusions from [37], which will be essential for the formulation of R 0 and our main findings in this study. 2) The following statements hold.
Let W (t, s, λ), t ≥ s be the evolution operator of (22) on R. Note that By Theorem 2.1 and 2.2 in [37], the ecological reproductive index can be also defined as σ 0 such that ρ(Φ (F/σ0)−V (ω)) = 1, which can be straightforward to calculate as follows The ecological reproductive index R 0 is a function of varying temperature and nutrient input. R 0 is a key indicator of phytoplankton fecundity measuring the average new born phytoplankton produced by one unit of phytoplankton during the average phytoplankton life span in one temperature period.
3.3. Boundary cyclic dynamics. In this subsection, we explore the boundary cyclic dynamics of (6) by establishing the existence and globally asymptotic stability of the boundary ω-periodic solution of (6).

3.4.
Internal cyclic dynamics. In the following, we investigate the internal cyclic dynamics of (6) by establishing some sufficient criteria for the existence and GAS of positive periodic solution of (6). For the reader's convenience, we shall summarize below a few concepts and results from [16] that will be essential for the following discussion. Let X, Z be normed vector spaces, L : Dom L ⊂ X → Z be a linear mapping, and N : X → Z be a continuous mapping. The mapping L will be called a Fredholm mapping of index zero if dim Ker L = codim Im L < +∞ and Im L is closed in Z. If L is a Fredholm mapping of index zero, then there exist continuous projectors P : X → X and Q : Z → Z such that Im P = Ker L, Im L = Ker Q = Im(I − Q). It follows that L| Dom L∩Ker P : (I −P )X → Im L is invertible. We denote the inverse of that map by K P . If Ω is an open bounded subset of X, the mapping N will be called L-compact on Ω if QN (Ω) is bounded and K P (I − Q)N : Ω → X is compact. Since Im Q is isomorphic to Ker L, there exist isomorphisms J : Im Q → Ker L. Theorem 3.9. Assume that f (P ) is strictly increasing with 0 ≤ f (P ) < 1 for P ≥ 0 and f (0) = 0. IfR then (6) has at least one strictly positive ω-periodic solution, where Proof. Considering the biological significance of (6), we specify A(0) > 0 and P (0) > 0. To facilitate the discussion below, we introduce the following notations Let A(t) = exp{x 1 (t)} and P (t) = x 2 (t), then (6) rewrites Define and x = (x 1 , x 2 ) T = max t∈[0,ω] |x 1 (t)| + max t∈[0,ω] |x 2 (t)| for any x ∈ X(orZ). Then X and Z are Banach spaces endowed with the norm · . Let Then, it is not difficult to show that Ker L = R 2 , Im L = {z ∈ Z : ω 0 z(t)dt = 0} is closed in Z, dim Ker L = 2 = codim Im L, and P , Q are continuous projectors such that Im P = Ker L, Ker Q = Im L = Im(I − Q). Therefore, L is a Fredholm mapping of index zero. Furthermore, the generalized inverse (to L)K P : Im L → Ker P ∩ Dom L is given by K P (z) = Let x = x(t) ∈ X be a solution of (27) for some certain λ ∈ (0, 1). Integrating (27) from 0 to ω produces From (27)- (31), it follows that By (28), one has bω exp( By (29), one has x 2 (ξ 2 ) ≤ P s , then ).
Next, we go ahead with investigating the globally asymptotic stability of the positive ω-periodic solution of (6). Theorem 3.11. Let (A * (t), P * (t)) be the positive periodic solution of (6). If where m i , M i , i = 1, 2, are defined in (8). Then (A * (t), P * (t)) is GAS.
Consider the Lyapunov function defined by A direct calculation of the right derivative of V (t) along the solutions of (6) leads to Note that there are some terms containing f (P (t))A(t) − f (P * (t))A * (t)) in the right-hand side of the above inequality and from (32), it follows that there exists a positive constant φ > 0 such that Integrating both sides of (33) from t 0 + T 1 to t (≥ t 0 + T 1 ) produces The boundedness of A * (t) and P * (t) and the ultimate boundedness of A(t) and P (t) imply that A(t), P (t), A * (t) and P * (t) all have bounded derivatives for t > t 0 + T 1 (from the equations satisfied by them). Then it follows that |A(t) − A * (t)| + |P (t) − P * (t)| is uniformly continuous on [t 0 + T 1 , +∞). By lemma 3.10, one has The proof is complete.
In summary, we have established some sufficient criteria for the boundary and internal cyclic dynamics of (6). Fig. 5 (a) and (b) illustrate such scenarios, where the parameter values are deliberately specified such that the criteria in Theorem  (6). It reveals that the conditions in Theorem 3.7 and Theorem 3.11 are sufficient for the cyclic dynamics and those two theorems can be further improved. (a) R 0 < 1 and the boundary ω-periodic solution is GAS (Theorems 3.7 is valid), here T (t) = 5−5 sin(π/6t+8), P in1 = 0.05. (b) R 0 > 1 and the internal ω-periodic solution is GAS (Theorem 3.11 holds) here T (t) = 15 − 15 sin(8πt + 8), P in1 = 5, µ opt = 5. (c) The conditions in Theorem 3.11 are not satisfied (i.e., R 1 < 1 < R 0 ), but (6) still admits a positive periodic solution being GAS, here T (t) = 15 − 15 sin(π/6t + 8), P in1 = 0.5. All the other parameters take the default values listed in Table 2. 3.7 and Theorem 3.11 are satisfied, respectively. In addition, it is not difficult to show that The criteria established above are of sufficient type and have some room for further improvement. Fig. 5 (c) illustrates such an example, where the conditions in Theorems 3.11 are not satisfied (i.e., R 1 < 1 < R 0 ) but the numerical simulation shows that (6) still admits a positive periodic solution being GAS. Moreover, in Fig. 6, the bifurcation diagrams of (6) with R 0 being the bifurcation parameter suggests that (6) admits a dichotomous cyclic dynamics, that is, the global attractor of (6) is the boundary ω-periodic solution if R 0 < 1, and is the internal positive ωperiodic solution if R 0 > 1. Now we shall propose a challenging conjecture that the ecological reproductive index R 0 = 1 is the threshold that completely characterize the global dynamics of (6). From Fig. 6, it is also observed that the amplitude of cyclic oscillation of phytoplankton increases with the increasing of R 0 and it is quite the opposite for nutrient,  Figure 6. Bifurcation diagram for (6) with R 0 being the bifurcation parameter. R 0 = 1 is the threshold that characterizes the global dynamics of (6), that is, the boundary ω-periodic solution is GAS (phytoplankton tends to perish) if R 0 < 1 and the internal ω-periodic solution is GAS (phytoplankton persists) if R 0 > 1.

Numerical analysis and applications.
In this section, based on the above theoretical analysis, by using the numerical simulations and field applications, we further investigate the ecological impact and consequence of temperature variation and the system's corresponding response.

4.1.
Effect of seasonal temperature. The general form of (6) and the definition of R 0 imply that the periodic dynamics produced by (6) are very robust to a wide array of periodic functions in that the qualitative dynamics can be reproduced independent of the qualitative form of the forcing. The seasonal patterns of phytoplankton follow characteristic climatology. This temporal structure of seasonal variability arises because the annual climate cycle is the primary driver of biomass variability. This finding biologically implies that both phytoplankton and nutrient variables displayed strong temporal variations strongly related to cycles of water temperature (Fig. 7) and their annual cycles can be tuned to match the climate cycle. In Fig. 7, some numerical experiments are performed to (6) in order to observe the seasonal periodic patterns of phytoplankton growth and the robustness to the temperature variation. Fig. 8 elaborates the relationship between the amplitude of phytoplankton cyclic oscillation and the cyclic forcing of water temperature. The figure shows the amplitude of the cyclic oscillation of phytoplankton increases first and then decreases with the increasing of the strength of cyclic forcing of water temperature. The relationship between the amplitude of temperature and the phytoplankton amplitude (b) Time-series plot of phytoplankton (A(t)) (solid line in red) and limiting nutrient (P (t)) (dotted line in black) in the reference model (6). It reveals that the annual cycle of temperature is the primary driver of the biomass variability. Here T = T (t) = 15−15 sin(π/6t+ 8) and other parameters take the default values in Table 2. is not simply monotonously consistent. The amplitude of phytoplankton oscillation achieves a maximum value at some intensity of forcing temperature.   (6) with respect to the strength or intensity (α) of external temperature forcing, which reveals a strong correlation between the water temperature and the cyclic dynamics of phytoplankton. Dashed line represents the bifurcation diagram for the amplitude (distance between two solid lines) of phytoplankton's periodic oscillation with respect to α. Here T (t) = 15 − α sin(π/6 · t + 8).

4.2.
Filed application in Lake Tai. Based on the analysis of 47 years observation Figure 9. Cyclic fluctuation pattern of monthly mean temperature of Lake Taihu (from [35]).
data of 4 weather stations (i.e., Dongshan, Wuzhong, Wuxi, and Yixing) around Lake Tai [35], Fig. 9 illustrates the trends of the monthly mean temperature of Lake Tai from 1961-2007, which develops a typical cyclic fluctuation pattern. Then, the modulation of water temperature T can be represented as a periodic solution T (t+ω) = T (t) with maxima at the optimal season, where ω is the period of forcing, which is taken to be monthly. A convenient and commonly applied scheme is that of sinusoidal forcing T (t) = T mean + α sin(Ωt + ϕ).
(34) Since the forcing is assumed to be monthly and t is in units of months, then Ω = π/6. We choose T 7 as a case study and we also simulate the whole lake by averaging the collected data. We fitted model parameters in (34) to the data sets of water temperature collected from the monitor stations T i (i = 1, 2, ..., 12) using the least square method. The fitted values of parameters in (34) are listed in Table 3 and the fitted curves are presented by solid lines in Fig. 10 (a) and Fig. 11 (a). Except the water temperature T (t), the nutrient P in1 takes the values in Table 3 and all the other parameters in (6) are held fixed at the values from published sources in the available experimental evidence of Lake Tai (see Table 2 for details) and are assumed to be constant across monitor stations and constant over time. The approach in this way has the advantage of providing additional check on the performance and adaptability of the model: does it produce patterns similar to the data, even though it was not fitted to those data?
According to studies on zooplankton composition and feeding rates of crustacean in Lake Tai [6,7], the zooplankton there could never reach values high enough to control Microcystis production. Moreover, the dominant fish population in Lake Tai was zooplankton-feeding fish and phytophagous fish abundance was insignificant to control the Microcystis bloom. That is to say, phytoplankton in Lake Tai should therefore not be controlled by the grazing of either zooplankton or fish. Therefore, the top-down control can be reasonably ignored and the bottom-up formulation is adequate as a general description in mathematical modeling. Therefore, the nutrient-phytoplankton model established in this study can be applied to model and to predict phytoplankton dynamics in Lake Tai.
In Lake Tai, temperature plays an important role in the phytoplankton composition and the dominant variations in phytoplankton composition are seasonal and are strongly related to the annual cycles of temperature [6,7]. It is reasonable to assume that the annual cycles of growth, reproduction, and senescence of phytoplankton are finely tuned to the annual climate cycle having a period of one year. Cyclic behavior of phytoplankton has been shown to occur in both the field and the laboratory cultures and reconciling the predictions of mathematical models with the nonequilibrium dynamics of experimental communities has proved to be challenging [15].  Figure 10. (a) Seasonal patterns of water temperature at monitor site T 7 . In the figures, dot-solid curve represents field data while solid curve comes from fitting. (b) Comparison of model predictions (solid line) with real experimental data (stars) for the monitor station T 7 , which shows the fundamental agreement of oscillatory population behavior and cycle phases between experiment and model prediction of Chlorophyll-a in Lake Tai.
As example, we only fit our model to monitor set T 7 and the average Chlorophylla of 12 monitor stations. From Fig. 10 (b), it is observed that the predictions of long-term dynamics of the parameterized model are in good agreement with the observed dynamics in both the oscillatory population behavior and cycle phases. Both the observed dynamics and the prediction of our model reveals a strong correlation between the water temperature and the cyclic dynamics of phytoplankton. With the rising of water temperature, the phytoplankton growth rare speeds up and the phytoplankton (represented by the density of Chlorophyll-a) increases accordingly. The water temperature not only affects the physicochemical process but also induces the variations of other environmental factors such as the PH value et al. The seasonal patterns of phytoplankton in Lake Tai follow characteristic climatology. This temporal structure of seasonal variability arises because the annual climate cycle is the primary driver of biomass variability. Our model reproduces this structure with input from the temperature and provides further insight into the mechanisms that led to the annually cyclic oscillation of phytoplankton population.
We treat Lake Tai as a whole and consider the average of the collected data in each monitoring station. Fig. 11 shows the average temperature and chlorophyll-a. The numerical simulation shows that the prediction curve well matches the changes of chlorophyll-a biomass. It is not difficult to observe that the fitness in Fig. 10 is better than that in Fig. 11, which is due to the fact that the whole lake contains more uncertain impact factors. It can be more accurate to predict a portion than a whole. Despite this, our model can simulate the trend of changes of phytoplankton biomass in the whole lake.  Figure 11. (a) Seasonal patterns of mean water temperature among 12 stations in Lake Tai. In the figures, dot-solid curve represents field data while solid curve comes from fitting. (b) Comparison of model predictions (solid line) with real experimental data (stars) for the average chlorophyll-a of 12 monitor stations, which shows the fundamental agreement of oscillatory population behavior and cycle phases between experiment and model prediction of Chlorophyll-a in Lake Tai.

4.3.
Short term response to the changing temperature. In order to better assess the model performance, adaptability, and to analyze the impact and consequences of temperature on model dynamics, below we present numerical simulations to show short time responses to the changing temperature.
We consider the water temperature being constant with an outbreak. In general, the temperature can be simulated by the following convenient and common formulation 15], where T mean is the mean temperature in the corresponding season, T high denotes the highest temperature which occurs in the 5th day, and δ is the parameter governing the width of 'open mouth' of the temperature variation. The dashed lines in green in Fig. 12 figure out the trend of temperature.
In all the four cases in Fig. 12, the density of phytoplankton also shows an outbreak under the corresponding temperature forcing. One can study Fig. 12 by comparing (a) with (b), (c), and (d). The wider 'open mouth' of temperature variation leads to the wider 'open mouth' of the phytoplankton outbreak ( Fig. 12 (a) and  (b)). In addition to this finding, we also find that the peak value of phytoplankton in (b) is higher than that in (a). That is to say, longer period high temperature causes the phytoplankton stay in higher value for long time. The magnitude of the temperature variation is also important. The bigger the temperature changes, the stronger the phytoplankton fluctuates ( Fig. 12 (a) and (c)). Even though with the same magnitude of the temperature variation, the mean temperature also plays an important role. One can consider the case (a) as summer and case (d) as winter.
The temperature with the same magnitude causes the different results ( Fig. 12 (a) and (d)). That is to say, phytoplankton bloom is more likely to happen during the summer. The four cases in Fig. 12 also illustrate that the peak of the density of phytoplankton occurs later than the time maximum temperature occurs ( Fig. 12  (a)). That is to say, temperature may have a lag effect on growth of phytoplankton. The time of maximum of phytoplankton falls in the range where T (t) decreases with time. It provides us warning the phytoplankton bloom after the high temperature. We would like to point out that the phytoplankton with different optimal temperatures can have different responses to even the same temperature forcing. The time series of phytoplankton with optimal temperature T opt = 23 o C is illustrated by the solid line in black in Fig. 13. For the phytoplankton with lower optimal temperature, something shows different. There will present a valley between two peaks when encountering the same temperature condition (Fig. 12 (a) and Fig. 13). It is not difficult to explain this phenomenon. Around the 5th day, the temperature is higher than the optimal temperature for growth of phytoplankton. The excessive temperature may have a negative effect on growth of phytoplankton. This result warns us that, even though the density of phytoplankton decreases, we should also be on alert of the outbreak of phytoplankton when the temperature is high. 5. Conclusion and discussion. Because of ignoring the time-varying nature of growth rate and environmental biomass of phytoplankton, the traditional autonomous model could not well explain why phytoplankton usually occur in high temperature seasons [28]. In this study, in order to investigate the effect of seasonal temperature on the growth of phytoplankton, by incorporating the time-varying feature of temperature, we propose a new time-dependent bottom-up nutrientphytoplankton interacting model. Our model provides a reasonable explanation to understand the mechanism of phytoplankton bloom. Our model has made some improvements of the traditional models. Effect of temperature on metabolism of phytoplankton like growth, nutrient uptake, natural mortality, and nutrient releasing from sediments are analyzed systematically and comprehensively. Especially, a more realistic model named CTMI is incorporated into our model. CTMI is determined by three parameters (T opt , T min and T max ) with biological significance depicting temperature effect on growth rate. An ecological reproductive index is formulated to establish the threshold result of (6). Through the analysis of the non-autonomous periodic model, the existence of positive periodic solution is proved by using coincidence degree theory and its global stability is achieved by constructing suitable Lyapunov function. Our theoretical findings work for very general functional response. Through the results of our simulation, the phytoplankton biomass oscillates periodically with the seasonal temperature and the amplitude of phytoplankton biomass is very robust to the variation of water temperature. But the oscillation amplitude of phytoplankton seems not positively correlated to the amplitude of temperature. The amplitude of phytoplankton cyclic oscillation increases first and then decreases with the increasing of the cyclic forcing of water temperature. In addition, the effect of daily temperature on the growth of phytoplankton is also discussed by some short-term numerical simulations. The results can explain some mechanisms of the phytoplankton bloom. The outbreak of phytoplankton biomass is related to the daily average temperature, the variation of temperature and how long it remains in high temperature. Our model without delay term can even simulate that the peak of phytoplankton biomass occurs a few days later than the peak of temperature. This gives a better understanding of the phytoplankton bloom often happens following the high temperature period. The numerical simulation and filed application are carried out for Lake Tai to assess the performance and the adaptability of our model. The dynamics of the model predictions agree well with our experimental data.
Certainly, in order to obtain a desired model to well expound the mechanism of temperature working on the growth of phytoplankton or even the triggering mechanism of PB, we still has a very long way to go. Freund [14] found that an average blooms are correlated with rapid upward temperature fluctuations and speculate on their possible role as trigger mechanisms. Daily average temperature may have more immediate effect on the phytoplankton growth and bloom. The effect of daily temperature on phytoplankton bloom should be studied intensively. To some extend, our model is still an idealized one, for example, here we view temperature as a periodic function of time, yet in the real world, the water temperature changes in some random way [28]. Hence, it is more realistic to take into account the stochastic noise in the equation predicting water temperature in our model and to implement further theoretical or numerical analysis.
Our study provides a preliminary explanation of how temperature contributes to the recurrent outbreaks of PB. To further understand the mechanism of PB, one should consider other abiotic and biotic factors that are responsible for PB. The light intensity, precipitation, ecological stoichiometry of growth-limiting nutrient (e.g., the ratio of nitrogen to phosphorus, expressed as N:P) etc., may also have significant effect on the metabolism of phytoplankton. Although our climate-driven bottom-up nutrient-phytoplankton model sheds some new light on the mechanism of phytoplankton dynamics, the alternative regulation of phytoplankton abundance can come from other biotic processes such as the important top-down control. The food web interactions also have a profound effect on the population dynamics of phytoplankton bloom species and it is left for our future effort.