Abstract
Plankton blooms are complex nonlinear phenomena whose occurrence can be described by the two-timescale (fast-slow) phytoplankton-zooplankton model introduced by Truscott and Brindley (Bulletin of Mathematical Biology 56(5):981–998, 1994). In their work, they observed that a sufficiently fast rise of the water temperature causes a critical transition from a low phytoplankton concentration to a single outburst: a so-called plankton bloom. However, the dynamical mechanism responsible for the observed transition has not been identified to the present day. Using techniques from geometric singular perturbation theory, we uncover the formerly overlooked rate-sensitive quasithreshold which is given by special trajectories called canards. The transition from low to high concentrations occurs when this rate-sensitive quasithreshold moves past the current state of the plankton system at some narrow critical range of warming rates. In this way, we identify rate-induced tipping as the underlying dynamical mechanism of largely unpredictable plankton blooms such as red tides, or more general, harmful algal blooms. Our findings explain the previously reported transitions to a single plankton bloom, and allow us to predict a new type of transition to a sequence of blooms for higher rates of warming. This could provide a possible mechanism of the observed increased frequency of harmful algal blooms.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
Introduction
Marine phytoplankton do not only form the basis of marine food webs and provide approximately half of the global primary production, they also contribute to essential biogeochemical processes in the ocean (Field et al. 1998; Falkowski 2012). A striking feature of many phytoplankton populations is the occurrence of so-called plankton blooms — a phenomenon characterized by a rapid increase in plankton concentrations followed by an almost equally rapid decline. Plankton blooms can occur seasonally (e.g., spring blooms) or in a more erratic fashion (e.g., red tides). A potentially catastrophic situation occurs when some of the bloom-forming species release toxic chemicals which can paralyze or even kill some affected marine species (Amaya et al. 2018; Griffith et al. 2019). The corresponding blooms are called harmful algal blooms (HABs).
In the past, various drivers of both types of blooms have been identified. Regarding spring blooms, the role of light and mixing processes has been emphasized for many decades (Gran and Braarud 1935; Riley 1946; Riley et al. 1949; Sverdrup 1953). The certainly most prominent hypothesis regarding the driving factor is the so-called critical depth hypothesis which has already been suggested in the middle of the last century (Gran and Braarud 1935; Sverdrup 1953) and which can still be found in more recently published textbooks (Simpson and Sharples 2012). Its key element constitutes the critical mixing depth above which improved upper ocean growth conditions allow for the accumulation of phytoplankton and thus ultimately for the formation of plankton blooms. Later, the role of biotic factors such as grazing pressure (Behrenfeld 2010; Behrenfeld et al. 2013), competition (Chakraborty and Feudel 2014; Busch et al. 2019), viral infection and parasitism (Chambouvet et al. 2008; Velo-Suárez et al. 2013; Richards 2017) or community composition (Lewandowska et al. 2015) attracted increasing attention. These studies include the disturbance-recovery hypothesis proposed by Behrenfeld (2010). The hypothesis presumes that blooms are initiated by recurring physical processes that disrupt the balance between phytoplankton reproduction and grazer consumption. According to the work of Behrenfeld and colleagues, this imbalance is caused by the annual deepening of the mixed layer which “dilutes” the grazing pressure on the phytoplankton which, in turn, allows for its significant accumulation (Behrenfeld 2010; Behrenfeld and Boss 2014).
In contrast to spring blooms, red tides or harmful algal blooms are more localized outbreaks which are not strongly correlated with a particular seasonal change in the upper ocean growth conditions. However similar to spring blooms, it is assumed that red tides occur due to an imbalance that affects the major contributing factors to phytoplankton growth: increasing light availability, rising temperature, nutrient enrichment and the reduction of the grazing pressure (Chakraborty and Feudel 2014; Truscott and Brindley 1994; Winder et al. 2012; Lewandowska et al. 2014; Hjerne et al. 2019; Trombetta et al. 2019; Sommer and Lewandowska 2011).
In order to fully understand the mechanisms underlying the emergence of annually occurring spring blooms and less predictable HABs such as red tides, it is not only important to identify the environmental conditions that induce the imbalance, but also to describe how this imbalance emerges. Interestingly, the magnitude of environmental change may not be the one decisive factor for the onset of a plankton bloom, which could disclose why the full explanation of this phenomenon has been elusive. The rate of change of environmental conditions also appears to play an important role (Pinek et al. 2020). The indications of rate dependence come from a study in which the Helgoland Roads long-term data series are analyzed (Freund et al. 2006). The study reveals that the spring bloom onset is correlated with a rapid increase of temperature, and not so much with the temperature itself.
The speed of temperature increase also plays a decisive role in the formation of red tides in the theoretical work of Truscott and Brindley (1994). Inspired by the idea of an excitable medium capable of mimicking fast growing harmful plankton blooms, they formulate a plankton model which consists of a fast evolving phytoplankton population that is controlled by a much slower reproducing zooplankton population. Depending on the speed of temperature increase, their time-scaled model shows two disparate transient behaviors: (i) balance of phytoplankton and zooplankton at low plankton densities (Fig. 1a) or (ii) the formation of a single plankton bloom (Fig. 1b). They further ascertain that both behaviors are delimited by a specific value of the speed of temperature increase (\(\upsilon \approx 0.003\) in Fig. 1). However, they were not able to identify the rate-sensitive quasithreshold that separates these two very different transient responses to temperature increase. In other words, the trigger mechanism of the harmful algal bloom has not been identified.
The phytoplankton (P)-zooplankton (Z) model developed by Truscott and Brindley (1994) is sensitive to rates of environmental change \(\upsilon\). It changes its dynamics from the stationary coexistence of P and Z at low densities. a to the sudden formation of a phytoplankton bloom b if solely the growth rate of the phytoplankton increases faster than the rate \(\upsilon \approx 0.003\) day\(^{-1}\) (as response to increasing temperatures). The phytoplankton density is displayed in green and red while the zooplankton density is shown in gray
Mathematically, the model of Truscott and Brindley (1994) consists of one fast variable, the phytoplankton, and two slow variables, the zooplankton and the time-varying environmental condition which changes at a given rate. In such two-timescale systems, exceptional solutions called canards are typical (Benoît 1981, 1983; Dumortier and Roussarie 1996; Szmolyan and Wechselberger 2001), and form boundaries between different dynamical regimes (Wieczorek et al. 2011; Wechselberger et al. 2013; Perryman and Wieczorek 2014). For instance, canards can separate small-amplitude oscillations from relaxation oscillations (Brøns et al. 2008; Desroches et al. 2012), different types of spiking behavior in neuronal models (Izhikevich 2007; Mitry et al. 2013), different states of CO\(_2\)-concentration in the atmosphere (Wieczorek et al. 2011; O’Sullivan et al. 2022), or disparate transient (Vanselow et al. 2019) or asymptotic (O’Keeffe and Wieczorek 2020) dynamics in biological systems characterized by different dominating species. In the latter four examples, similar to the observations of Truscott and Brindley (1994), a variation of the rate of environmental change alone can cause a transition from one dynamical regime to another by crossing canard solutions. These critical transitions have been classified as rate-induced tipping points by Ashwin et al. (2012). Others identified so-called folded-saddle canards as non-obvious thresholds for the onset of rate-induced tipping (Wieczorek et al. 2011; Mitry et al. 2013; Perryman and Wieczorek 2014; O’Sullivan et al. 2022). Specifically, as the rate of environmental change increases, the position of the rate-dependent threshold shifts in state space and might pass the current state of the system, giving rise to a large nonlinear response. It is important to note that in these cases there is no classical bifurcation (no loss of stability in the classical autonomous sense), such as the dangerous saddle-node bifurcation, separating different dynamical regimes. Instead, the different dynamical regimes are solely separated by rate-sensitive canard solutions.
In this work, we uncover the rate-sensitive boundary, that is the singular canard solution, in the model of Truscott and Brindley (1994). Furthermore, we demonstrate that this boundary gives rise to a quasithreshold for rate-induced tipping, which is the dynamical mechanism responsible for the observed harmful algal blooms (Fig. 1). To this end, we start by introducing the phytoplankton-zooplankton model of Truscott and Brindley (1994) (Sec. “The model”). Afterwards, we reproduce their simulations in which the phytoplankton population is exposed to rapidly increasing growth rates induced by increasing temperatures. Then, we study the resulting rate-sensitive dynamics using geometric singular perturbation theory, and reveal the singular canard solution (Sec. “Rate-induced tipping triggers plankton bloom”). Moreover, we show that, depending on the location of the singular canard, the model can show more than one sequentially occurring plankton bloom (Sec. “R-tipping can trigger two sequentially occurring plankton blooms”). Next, we relate this singular canard to a family of canards that form a quasithreshold. Finally, we discuss our results (Sec. “Discussion”).
The model
The Truscott-Brindley model (Truscott and Brindley 1994) (in the following the TB-model) combines a logistic growth of a phytoplankton population P, a zooplankton population Z growing by the Holling-type III functional response (Holling 1959a, b), and dying with standard linear mortality.
In the following, we vary exclusively the growth rate r of the phytoplankton population while the other parameters are kept constant. The constant parameters K, \(R_m\), \(\alpha\), \(\gamma\) and \(\mu\) represent the carrying capacity of the phytoplankton population, the attack rate, the half-saturation phytoplankton concentration, the conversion efficiency and the mortality rate of zooplankton population (see Table 1 for the parameter values).
We will analyze the response of the plankton model to a gradual increase of the phytoplankton growth rate r, which could be caused by a gradual increase in temperature. Specifically, we consider the simplest case in which r increases linearly over time at a rate \(\eta\) as long as \(r\in (r_{min},r_{max})\), and remains constant when \(r=r_{min}\) or \(r=r_{max}\) (Fig. 2a):
In the extended TB-model (1)–(3), the phytoplankton population P evolves on a different time scale than its grazer Z. This can be demonstrated by transforming the model into a non-dimensional form, which we show in detail in Appendix A.1. Considering the non-dimensional TB-model Eqs. (A5), (A6) and (A9), it becomes apparent that the time scale separation is determined by the conversion efficiency \(\gamma\) of the zooplankton which quantifies its turnover of phytoplankton biomass: the lower the turnover, the slower the zooplankton density changes in time. To emphasize that the change of the growth rate r is as slow as the development of zooplankton we make the ansatz, that \(\eta =\gamma \upsilon\). Now we discuss the representation of the model taking this time scale separation into account. Following Truscott and Brindley (1994), we fix the conversion efficiency to \(\gamma = 0.05\). Therefore, the phytoplankton population evolves on the fast time scale \(\tau\), while the zooplankton density and r change on the slow time scale t. Using the relationship between these two timescales \(t = \gamma \tau\) we can write the extended dimensional TB-model (1)–(3) in slow time t:
with \(\kappa =\frac{\mu }{\gamma }\). From now on we will analyze (4), (5) and (6) and fix all parameters K, \(R_m\), \(\alpha\), \(\gamma\) and \(\mu\) according to Table 1, which is based on the original parameterization in Truscott and Brindley (1994) following the reasonable values formulated in (Uye 1986; Wake 1991). To allow for an easy comparison between our results and the results obtained by Truscott and Brindley (1994), we again employ the original formulation of the phytoplankton-zooplankton model (4)–(6) for the remainder of this work. Furthermore, we assume that the phytoplankton and zooplankton populations coexist at low densities in a stable equilibrium at \(t=0\) before the growth rate of the phytoplankton r starts to increase linearly in time at the rate \(\upsilon\).
Regular and quasi-static equilibria
When growth rate r is fixed, the autonomous TB-model (1)–(2) has three stationary solutions (equilibria). In addition to the extinction and the phytoplankton-only equilibria, both of which are unstable, there is a stable coexistence equilibrium in the (P, Z) phase plane (see Appendix A.2 for more details):
This equilibrium corresponds to the balance between the phytoplankton P and zooplankton Z populations, and is the starting point for our analysis of harmful algal blooms. It is always taken as the initial condition for the simulations. Most importantly, this equilibrium is linearly stable for the chosen parameter settings (Table 1) and all \(r \in (r_{min}, \; r_{max})\). In other words, e(r) never bifurcates in the classical autonomous model within the range of phytoplankton growth rates considered here.
The growth rate r changes over time, e.g., with increasing temperatures, and we need to consider the extended TB-model (4)–(6), where r becomes an additional dynamical variable. When \(r = r_{max}\) or \(r=r_{min}\), the long-term behavior of the TB-model is determined by the unique stable equilibrium \(e(r=r_{max})\) or \(e(r=r_{min})\), respectively. However, when \(r \in ( r_{min}, r_{max})\), the system has no equilibrium solutions because \(\frac{\textrm{d}r}{\textrm{d}t} =\upsilon >0\). In this range, we will be interested in how the position of the state of stable coexistence \(e(\upsilon t)\) changes with time, and how the system evolves relative to this changing \(e(\upsilon t)\). Note that \(e(\upsilon t)\) is referred to as a quasi-static equilibrium (Ashwin et al. 2012; Wieczorek et al. 2011), or a moving equilibrium (Vanselow et al. 2019; O’Keeffe and Wieczorek 2020; Wieczorek et al. 2021), which is computed by inserting \(r=r(t)=\upsilon t\) into ((7)). The coexistence of the phytoplankton and zooplankton population at low densities is guaranteed when a solution stays near \(e(\upsilon t)\) for small enough growth rates \(\upsilon > 0\). In this case, the enhancing growth conditions of the phytoplankton are compensated by an increasing grazing pressure (increasing r causes an increase of \(e_Z(r)\), Fig. 1a). However, an interesting instability occurs for larger rates of change \(\upsilon\). In this case, the top-down control by the zooplankton relaxes since it is no longer able to balance the faster growing phytoplankton. This provokes an imbalance which results in a deviation of the solution of the extended TB-Model (1)–(3) from the quasi-static equilibrium \(e(\upsilon t)\), manifesting in the formation of a phytoplankton bloom (Fig. 1b). This instability cannot be explained by the classical autonomous bifurcation theory, and requires an alternative approach. To fully understand this instability, we recall some basic concepts from the geometric singular perturbation theory for fast-slow systems.
The critical manifold
As shown in the previous Section, the system under investigation possesses a timescale separation. In fact, the dynamics of this 1-fast 2-slow system is slow for most of the time (Fig. 2b, single-headed arrows) with short time intervals of fast motion (double-headed arrows). Accordingly, it is reasonable to begin the study of our system by approximating the slow dynamics by taking \(\gamma \rightarrow 0\). Note that we take \(\gamma \rightarrow 0\) while keeping \(\kappa\) constant, which requires that \(\mu \rightarrow 0\).
a The growth rate of the phytoplankton r increases linearly in time at the rate \(\upsilon\) between \(r(0)=r_0>r_{min}\) and \(r_{max}\) (6). b The stable (red), unstable parts (blue) and folds \(F_1(r)\), \(F_2(r)\) (black lines) of the critical manifold S (11) organize the slow-fast motion in phase space. Notice that S approximates the dynamics in the limit \(\gamma \rightarrow 0\). Parameters of red trajectory: a \(r(0) = r_0 = 0.4\), \(\upsilon = 0.004\), b, c \(P(0) = P_0 = e_P\), \(Z(0) = Z_0 = e_{Z}\), \(r(0) =r_0=0.4\), \(\upsilon = 0.004\), \(Rm = 0.4\), \(\gamma = 0.05\). Other parameters see Table 1
Taking the limit \(\gamma \rightarrow 0\) for the slow time scale t gives the reduced model
for the evolution of the slow variables Z and r in slow time t on the so-called critical manifold S:
in the (P, Z, r) phase space. Rearranging the condition above with respect to the Z-coordinate gives the cubic critical manifold formula
Besides its approximation of the slow flow (one-headed arrow, Fig. 2b), the critical manifold S organizes the fast dynamics as follows: the (red) stable parts of S attract the fast flow, whereas the (blue) unstable part of S repels the fast flow. The stable and unstable parts connect at the two folds \(F_1(r)\) and \(F_2(r)\) (black lines) of S (see Appendix A.3 for the formulation of the folds). In the following, we display the critical manifold and the trajectories of the TB-model together in several figures. Notice that the critical manifold approximates the dynamics in the limit \(\gamma = 0\), whereas the trajectories represent solutions of the TB-model for \(0 < \gamma \ll 1\).
Rate-induced tipping triggers plankton bloom
If \(\upsilon =0.004\), the extended TB-model (4)–(6) in slow time t reveals a plankton bloom (Fig. 1b, red), but it does not if \(\upsilon =0.002\) (Fig. 1a, green). Since the difference in \(\upsilon\) is the only difference between the two setups (see Table 1 for the remaining parameters and the initial conditions), this indicates that, between the rates \(\upsilon =0.002\) and \(\upsilon =0.004\), the system crosses some quasithreshold that separates both dynamical regimes. To comprehend what creates the qualitative different behavior of the green and red trajectories (Fig. 1), we examine both trajectories in the (P, Z, r) phase space (Fig. 3). Since the stable (red), unstable parts (blue) and the folds (black solid lines) of the critical manifold organize the flow in phase space, we add the critical manifold S of the extended TB-model (4)–(6) to the phase portraits.
Starting at \(e(\upsilon t)\) (7) at \(t=0\), the green trajectory (\(\upsilon = 0.002\)) reveals the expected behavior: it slowly follows the pathway of \(e(\upsilon t)\) (gray dotted line, Fig. 3a) until it settles on the stable long-term state \(e(r = r_{max})\). In this case, the slight increase in the zooplankton density, which implies an increasing grazing pressure, is sufficient to balance the growth of the phytoplankton. Consequently, we observe no plankton bloom. Since the green trajectory remains close to the quasi-static state \(e(\upsilon t)\) (7) at all times, we say that the system tracks the stable quasi-static state \(e(\upsilon t)\).
On the contrary, the red trajectory (\(\upsilon = 0.004\)) leaves the vicinity of the stable quasi-static state \(e(\upsilon t)\) and performs a large excursion before converging to \(e(r_{max})\) (Fig. 3b). At first, the red trajectory slowly proceeds towards the fold \(F_1(r)\) (black solid line). In the vicinity of \(F_1(r)\), its motion changes from slow to fast and remains fast until it reaches the vicinity of the right stable part of S. This fast motion away from the quasi-static state \(e(\upsilon t)\) entails a rapid increase in the phytoplankton density — the plankton bloom arises. When the trajectory reaches the vicinity of the right stable part of S, the plankton bloom possesses its maximal concentration. From here, the red trajectory slowly proceeds towards lower phytoplankton concentrations along the right stable part. Hence, the bloom slowly starts to recede due to an increasing zooplankton concentration controlling it. When the red trajectory approaches the second fold \(F_2(r)\), it starts to move fast towards the left stable part: the bloom collapses rapidly due to the still persisting high zooplankton concentration. Along the left stable part, it finally converges to the stable long-term state \(e(r_{max})\). Since the zooplankton feeds on the phytoplankton on a slower time scale than the phytoplankton reproduces, its evolution in time always lags behind the development of the phytoplankton. For instance, the zooplankton bloom reaches its peak when the phytoplankton density is close to zero (Figs. 1b and Fig. 3b).
Dynamics of the three-dimensional TB-model (4)–(6) for two different rates \(\upsilon = 0.002\) and \(\upsilon = 0.004\). Stable, unstable part and folds \(F_1(r)\), \(F_2(r)\) of the critical manifold S (11) are shown in red, blue and as solid black line. a The green trajectory remains close to the quasi-static equilibrium \(e(\upsilon t)\) (7) (gray dotted line). b The red trajectory leaves the vicinity of \(e(\upsilon t)\) close to the fold \(F_1(r)\) towards the fast P-direction causing the formation of the plankton bloom. \(P_0 = e_P\), \(Z_0 = e_Z\) and \(r_0 = 0.4\). Parameters see Table 1
Hence, in order to better understand the initiating mechanism, we examine the evolution of the fast variable P on the critical manifold S. To this end, we consider the reduced system, which describes the evolution of the slow variables Z and r on the critical manifold S. On the basis of the reduced system, we find a three-dimensional system that also captures the evolution of the fast variable P on S. Finally, for the sake of simplicity, we lower the complexity of the resulting system by projecting its three-dimensional flow on the two-dimensional critical manifold S using \(Z = h(P,r)\) (see Appendix B for the transformation of the extended TB-model (4)–(6) to the projected-reduced system (13)–(14) and also for the detailed formulation of the numerator \(\Lambda (P,r,\upsilon )\) and the denominator \(\frac{\partial f^S(P,r)}{\partial P}\Big \vert _S\) in Eq. (13)). The resulting projected-reduced system is given by:
Since the plankton bloom arises when the dynamics of the TB-model (4)–(6) changes from slow to fast motion close to the fold \(F_1(r)\), we start our analysis of the projected-reduced system (13)–(14) on the stable part of S near the fold \(F_1(r)\) and vary the phytoplankton growth rate r. Depending on the sign of the numerator \(\Lambda (P_{F1},r,\upsilon )\) which changes with r, we can distinguish three different types of trajectories (see Fig. 4 and Appendix B for more details): (i) For growth rates where \(\Lambda (P_{F1},r,\upsilon )<0\), (red) trajectories are attracted to the fold \(F_1(r)\). However, at the fold \(F_1(r)\), the denominator \(\frac{\partial f(P^S_{F1},r,\upsilon )}{\partial P}\Big \vert _S\) becomes zero. Hence, solutions of the projected-reduced system P(t) blow up (go to infinity in finite time t) when they reach typical points on the fold \(F_1(r)\). In other words, solutions cease to exist within S when they reach typical points on \(F_1(r)\) (Fig. 4b, d). (ii) For values of r where \(\Lambda (P_{F1},r,\upsilon ) > 0\), the (green) trajectories never reach \(F_1(r)\) because they are repelled from the fold (Fig. 4a, c). (iii) There can be special points (special values of r) along the fold \(F_1(r)\) at which both the numerator \(\Lambda (P_{F1},r,\upsilon )= 0\) and the denominator become zero such that \(\frac{\textrm{d}P}{\textrm{d}t}\) remains finite (see Appendix B for more details). The corresponding (blue) trajectory approaches such a point on the fold \(F_1(r)\) slowly and is able to cross it with finite speed (Fig. 4). Afterwards it proceeds slowly along the unstable part of S without being repelled in the fast P-direction (Fig 4c). Hence, this special trajectory — called singular canard — combines aspects of both dynamical regimes, i.e., of the green and red trajectory: moving away from the quasi-static state (red, bloom) but slowly (green, no bloom). The special fold point \(FS = (P_{FS}, r_{FS})\) (Fig. 4c), at which \(\Lambda (P_{F1},r,\upsilon )=0\), is called the folded saddle singularity (Szmolyan and Wechselberger 2001) (see Appendix B for the computation of the folded saddle and the canard). Since the term \(\Lambda (P_{F1},r, \upsilon )\) depends on the rate \(\upsilon\), the position of the folded saddle and thus the position of the singular canard on the critical manifold S also depends on the rate \(\upsilon\) (see Appendix B). For this reason, the location of the singular canard changes from the scenario without bloom (\(\upsilon = 0.002\), green trajectory) to the scenario with plankton bloom (\(\upsilon = 0.004\), red trajectory). For \(\upsilon = 0.002\), the singular canard is located below the initial condition of the green trajectory. Fold points above the singular canard, respectively above the folded saddle FS (\(P_F > P_{FS}\) and \(r_F >r_{FS}\)), are repelling fold points (\(\Lambda > 0\), Fig. 4c). Consequently, the green trajectory is repelled by the part of the fold \(F_1(r)\) above the singular canard and tracks the stable quasi-static state \(e(\upsilon t)\): no bloom is formed. For \(\upsilon = 0.004\), the singular canard is located at higher values of the growth rate r. As a consequence, the singular canard is now located above the initial condition of the red trajectory. The fold points below the singular canard are attracting \((\Lambda < 0)\) and therefore, the red trajectory crosses the fold and runs away in the fast direction. As a result, a plankton bloom emerges (\(\Lambda < 0\), Fig. 4d). Thus, in the limit \(\gamma \rightarrow 0\), there exists an isolated critical rate \(\upsilon =\upsilon _{crit}\) for each initial condition. This critical rate is the value of \(\upsilon\) for which the given initial condition lies exactly on the \(\upsilon\)-dependent singular canard. This singular canard can be thought of as a singular threshold for rate-induced tipping.
a, c The green trajectory starts above the canard (blue) where it is repelled by the fold points of \(F_1(r)\) (\(\Lambda >0\), black line) and returns to \(e(\upsilon t)\) — no bloom emerges. b, d The red trajectory starts below the canard where it is attracted by the fold points (\(\Lambda <0\)), reaches them and moves fast away resulting in the large excursion away from \(e(\upsilon t)\) — the bloom occurs. The canard (blue trajectory) crosses the fold \(F_1(r)\) at the folded saddle FS where \(\Lambda = 0\). Parameter: \(P_0 = e_{P}\), \(Z_0 = e_{Z}\), \(r_0=0.4\). Other parameters see Tab. ()
The critical manifold S only approximates the slow dynamics in the limit when the time scale separation is perfect, i.e., \(\gamma = 0\) (Wechselberger et al. 2013). However, in the full extended TB-model (4)–(6), we consider a finite time scale separation between phytoplankton and zooplankton, and, hence, \(0 < \gamma \ll 1\). To evaluate if the singular canard persists in the full extended TB-model (neither reduced nor projected), we need to translate the dynamics in the singular limit (\(\gamma = 0\)) to the full dynamics (\(0 < \gamma \ll 1\)). In the full system, the critical manifold S is replaced with a nearby slow manifold \(S_\gamma\), and the stable and unstable parts of \(S_\gamma\) typically split along the fold (Fig. 5b). They only intersect near the point FS (white circle in Fig. 5a) which represents the folded saddle in the singular limit system. This intersection point gives rise to a maximal canard that crosses from the stable part of the slow manifold into the unstable part of the slow manifold where it stays for as long as the unstable part exists (Fig. 5b). On the lower-r side of the maximal canard, the unstable part is located above the stable part. Hence, (red) trajectories can proceed beneath the unstable part and then run away in the fast direction — a bloom arises. Conversely, on the higher-r side of the maximal canard, the unstable part is located below the stable part whereby (green) trajectories proceed above the unstable part and are thus repelled back towards the left stable part — P and Z remain at low densities.
a Stable (red) and unstable (blue) part of the critical manifold S merge along the fold \(F_1(r)\) (black solid line). (Red) trajectories cease to exist at the fold whereas a singular canard (blue trajectory) is able to cross the fold via the point FS. (Green) trajectories are repelled by the fold and approach \(e(\upsilon t)\) (gray dashed line). b Stable and unstable part are perturbed along the position of the fold \(F_1(r)\) for \(\gamma = 0\). At their intersection point (FS when \(\gamma \rightarrow 0\)), the maximal canard crosses from stable to unstable part separating (green) trajectories that are repelled by the unstable part from (red) trajectories that leave the stable part in the fast P-direction
Therefore, also in the full system (\(0<\gamma \ll 1\)), there is a boundary separating (green) trajectories, reflecting the maintenance of the balance between P and Z, from (red) trajectories, that show the formation of a plankton bloom. The difference is that, in contrast to the singular case (\(\gamma = 0\)), this boundary is not clear cut. This can be understood as follows. In addition to the maximal canard, which is in a one-to-one correspondence with the singular canard (Wechselberger et al. 2013) and follows the unstable part of the slow manifold for the longest time, there are additional canards. These additional canards stay close to the maximal canard for some shorter time, after which they leave the unstable part of the slow manifold towards the left stable part. This whole family of canards is responsible for the transition to a plankton bloom, and is referred to as a quasithreshold for rate-induced tipping (Wieczorek et al. 2021; O’Sullivan et al. 2022). Crossing such a quasithreshold occurs for a narrow critical range of \(\upsilon\) rather than at one isolated critical rate \(\upsilon =\upsilon _{crit}\).
In summary, we have demonstrated that the harmful algal bloom or red tide is solely triggered by rate-induced tipping (Ashwin et al. 2012; Wieczorek et al. 2011). An increase of the rate of change \(\upsilon\), which describes the speed at which the growth rate of the phytoplankton r increases, changes the position of the maximal canard on the slow manifold \(S_\gamma\) in a way that causes the (red) trajectory to leave the region nearby \(e(\upsilon t)\), resulting in the formation of a plankton bloom in the full system (Fig. 4b).
R-tipping can trigger two sequentially occurring plankton blooms
Interestingly, the extended TB-model (4)–(6) can show more than one plankton bloom (Fig. 6). In the previous section, we have demonstrated that if initial conditions of the extended TB-model (4)–(6) are located on \(S_\gamma\) on the lower-r side of the maximal canard, the model shows a rate-induced plankton bloom. Clearly, when the red trajectory still finds itself on the lower-r side of the maximal canard following the bloom, it can exhibit another bloom before it settles on the long-term state \(e(r_{max})\) (Fig. 6). The number of possibly recurring blooms naturally increases with parameter changes that shift the maximal canard towards higher values of r accompanying a higher maximum growth rate \(r_{max}\). For instance, increasing the maximum growth rate \(r_{max}\) from 0.6 to 0.8 and simultaneously increasing the rate of change \(\upsilon\) to 0.006 triggers two recurring plankton blooms (Fig. 6). Other parameter changes that promote the occurrence of multiple blooms are an increase of the zooplankton’s mortality \(\omega\) or a decrease of the zooplankton’s attack rate \(R_m\) (see Appendix C for more details). Obviously, parameter changes reducing the grazing pressure and therefore relaxing the top-down control by the zooplankton encourage the formation of more than one bloom.
Two sequentially recurring plankton blooms occur if the red trajectory reaches the left stable part of the critical manifold S (red) beneath the canard (blue) after performing the first excursion. Parameter: \(P_0 = e_{P}\), \(Z_0 = e_{Z}\), \(r_0=0.4\), \(r_{max} = 0.8\). Other parameters see Table 1
Discussion
In ecology, several factors which are responsible for the emergence of regularly and irregularly occurring plankton blooms have already been identified (Behrenfeld 2014; Sommer et al. 2012; Lewandowska et al. 2015; Freund et al. 2006; Chakraborty and Feudel 2014) but still, some mechanisms are not fully understood. Hence, there exists some uncertainty how future climate change will impact the severity and frequency of plankton blooms (Doney 2006; Hillebrand et al. 2018; Winder et al. 2012). We demonstrate that rate-induced tipping constitutes a possible mechanism explaining the occurrence of erratically occurring blooms such as red tides. We proceed from the work of Truscott and Brindley (1994) in which they developed a time-scaled phytoplankton (P, fast variable)-zooplankton(Z, slow variable) model using ideas from the theory of excitable media. In their work, they studied the response of the model to fast environmental changes that provoke an increase of the phytoplankton’s growth rate. Depending on its speed of change, the model reveals two disparate behaviors: conservation of the balance of P and Z (no bloom), or an imbalance of both that manifests in a plankton bloom. Using fast-slow system theory, we uncover the quasithreshold phenomenon which separates both behaviors in state space: a special trajectory called maximal canard together with a family of shorter canards. The position of the quasithreshold determines whether solutions remain close to the balance state corresponding to low plankton concentrations or leave its vicinity leading to an imbalance that manifests in a large plankton bloom. Since the location of the quasithreshold depends on the speed, or the rate, at which the phytoplankton growth rate increases, rate-induced tipping constitutes the mechanism being responsible for the emergence of the harmful algal bloom. We further demonstrate that decreasing the grazing pressure or broadening the interval in which the growth rate increases allows for multiple recurring plankton blooms.
In accordance with Truscott and Brindley (1994), we assume that the rate-induced HAB is triggered by fast enhancing growth conditions of the phytoplankton due to temperatures increasing at a certain speed. However, other scenarios in which the growth conditions improve gradually at a certain speed are just as conceivable, e.g., due to an increase of the light intensity (Rumyantseva et al. 2019; Winder et al. 2012) or the nutrient supply (Largier 2020; Guseva and Feudel 2020). Moreover, an imbalance between phytoplankton growth and predator control can be obtained if traits of the grazers, such as attack rate or mortality rate, are affected by rapid environmental changes (Behrenfeld 2010; Busch et al. 2019). Consequently, the notion that rate-induced tipping is able to cause irregularly occurring plankton blooms does not depend on any particular environmental driver but rather on its dynamics: As long as any environmental condition changes at a certain speed, rate-induced tipping is a potential trigger mechanism of plankton blooms. The same argument holds also for regular spring blooms as long as their grazing function is described with the same functional response as for harmful algal species.
A second precondition for the occurrence of rate-induced blooms is that phytoplankton reproduces faster than their grazers. This condition is often met since the hourly to daily cell division of phytoplankton (Franks 2001) is typically much faster than the reproduction of their grazers which mainly possess generation times from hours to a month (Hirsche 2013; Klais et al. 2016). Empirical evidence for the importance of time scale separation in the formation of plankton blooms comes from iron enrichment experiments carried out in the 1990s (Cavender-Bares et al. 1999; Cullen 1995; Morel et al. 1991). A particularly vivid example was given by the IronEx II experiment (Coale et al. 1996), in which iron enrichment only led to about a doubling of the picoplankton biomass, but up to an 85-fold increase in some diatom species (phytoplankton). Picoplankton are usually grazed by protists which have often short response time scales — similar to the generation times of the picoplankton. Accordingly, protists can keep the picoplankton at low densities. On the contrary, phytoplankton like diatoms tend to be grazed by metazoan zooplankton, which have relatively long generation times and therefore long response times which allows for the formation of temporary imbalances and thus for the formation of plankton blooms (Franks 2001).
Finally, we want to point out two shortcomings of the current state of the art in ecological modeling concerning the impact of environmental change. We have seen that the observation of rate-induced phenomena, like the formation of plankton blooms, depends on (i) the evolution of the time horizon of environmental change and (ii) the presence of multiple time scales in species’ development. (i) Regarding the time scale of environmental change, most studies employ the concept called “tipping point” which underlies the idea of a catastrophic bifurcation of the involved long-term states if a critical threshold is exceeded (Scheffer et al. 2001; Lenton et al. 2008; Lenton 2020; Krönke et al. 2020; Dakos et al. 2019). Bifurcation analysis implies that environmental conditions change infinitely slow. Hence, tipping points are only relevant in cases where it can be assumed that the speed of environmental change is much slower than the intrinsic ecological dynamics. By contrast, rate-induced tipping occurs on time scales comparable to the dynamics of ecosystems and, hence becomes apparent only in the transient dynamics — the dynamics prior to the long-term dynamics (Hastings et al. 2018; Morozov et al. 2020). Hence, the plankton bloom described in this work would be overlooked if we would focus exclusively on the long-term dynamics.
Concerning the omission of the presence of multiple time scales: if the presence of multiple time scales is neglected, no canard solution can exist and thus, without the existence of any other instability (e.g., a tipping point), no rate-induced phenomena can be observed in the system. Such phenomena are, however, crucial when examining the effect of the speed of environmental change as they reveal critical states which would be missed otherwise. Since present climate change accelerates environmental changes, corresponding studies will be of paramount importance in the future (Walther et al. 2002; Parmesan 2006; Smith et al. 2015).
In summary, to meet the accelerating speed of environmental change, we suggest three key elements for consideration in future ecological studies. First of all, we believe that it will be of paramount importance to take into account the temporal evolution of environmental disturbances or changes (which is usually not infinitely slow). Secondly, in order to capture the full impact of fast rates of environmental change, it is essential to explicitly consider the different intrinsic time scales which are present in an ecological system. Thirdly, since critical phenomena — like rate-induced tipping — act on time scales comparable to the system dynamics, the evaluation of environmental impacts should not be exclusively based on the long-term response but on the transient dynamics as well. In addition we would like to emphasize that an experimental test for rate-induced transitions in the laboratory requires explicit measurement of the rate of environmental change.
For certain, the Truscott-Brindley model misses essential physical processes, such as vertical mixing and sinking, as well as other ecological processes, such as viral infection and the diverse composition of interacting communities. However, just like any theoretical model, due to its delightful simplicity it enables the understanding of otherwise unsolvable relationships between ecological actors and their environment. Our findings suggest that plankton communities, which typically involve species evolving on multiple time scales, are potentially prone to environmental disturbances evolving at a certain speed. Hence, we propose to consider multiple time scales in theoretical models of ecosystems, such as food webs, and to examine their transient dynamics.
Availability of data and materials
Not applicable.
Code availability
Not applicable.
References
Amaya O, Quintanilla R, Stacy BA et al (2018) Large-scale sea turtle mortality events in El Salvador attributed to paralytic shellfish toxin-producing algae blooms. Front Mar Sci 5:411
Ashwin P, Wieczorek S, Vitolo R et al (2012) Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system. Phil Trans R Soc A 370:1166–1184
Behrenfeld MJ (2010) Abandoning sverdrup’s critical depth hypothesis on phytoplankton blooms. Ecology 91(4):977–989. https://doi.org/10.1890/09-1207.1
Behrenfeld MJ (2014) Climate-mediated dance of the plankton. Nat Clim Chang 4(10):880–887
Behrenfeld MJ, Boss ES (2014) Resurrecting the ecological underpinnings of ocean plankton blooms. Ann Rev Mar Sci 6:167–194
Behrenfeld MJ, Doney SC, Lima I et al (2013) Annual cycles of ecological disturbance and recovery underlying the subarctic Atlantic spring plankton bloom. Global Biogeochem Cycles 27(2):526–540
Benoît E (1981) Chasse au canard. Collect Math 31–32:37–119
Benoît E (1983) Systemes lents-rapides dans \(\mathbb{R} ^3\) et leur canards. Asterisque 109–110:159–191
Brøns M, Kaper T, Rotstein H (2008) Introduction to focus issue: mixed mode oscillations: experiement, computation, and analysis. Chaos 18(015):101
Busch M, Caron D, Moorthi S (2019) Growth and grazing control of the dinoflagellate Lingulodinium polyedrum in a natural plankton community. Mar Ecol Prog Ser 611:45–58
Cavender-Bares KK, Mann EL, Chisholm SW et al (1999) Differential response of equatorial Pacific phytoplankton to iron fertilization. Limnol Oceanogr 44(2):237–246
Chakraborty S, Feudel U (2014) Harmful algal blooms: combining excitability and competition. Thyroid Res 7:221–237
Chambouvet A, Morin P, Marie D et al (2008) Control of toxic marine dinoflagellate blooms by serial parasitic killers. Science 322(5905):1164,387
Coale KH, Johnson KS, Fitzwater SE et al (1996) A massive phytoplankton bloom induced by an ecosystem-scale fertilization experiment in the equatorial Pacific Ocean. Nature 383:495–501
Cullen J (1995) Status of the iron hypothesis after the open-ocean enrichment. Limnol Oceanogr 40:1336–1343
Dakos V, Matthews B, Hendry AP et al (2019) Ecosystem tipping points in an evolving world. Nat Ecol Evol 3:355–362
Desroches M, Guckenheimer J, Krauskopf B et al (2012) Mixed-mode oscillations with multiple time scales. SIAM Rev 54(2):211–288
Doney SC (2006) Plankton in a warmer world. Nature 444:695–696
Dumortier F, Roussarie R (1996) Canard cycles and Center Manifolds, 577th edn. Memories of the American Mathematical Society
Falkowski P (2012) The power of plankton. Nature 483:7–10
Field CB, Behrenfeld MJ, Randerson JT et al (1998) Primary production of the biosphere: Integrating terrestrial and oceanic components. Science 281(5374):237–240
Franks PJS (2001) Phytoplankton blooms in a fluctuating environment: the roles of plankton response time scales and grazing. J Plankton Res 23(12):1433–1441
Freund JA, Mieruch S, Scholze B et al (2006) Bloom dynamics in a seasonally forced phytoplankton-zooplankton model: Trigger mechanisms and timing effects. Ecol Complex 3(2):129–139. https://doi.org/10.1016/j.ecocom.2005.11.001
Gran H, Braarud T (1935) A quantitative study on the phytoplankton on the Bay of Fundy and Gulf of Maine (including observations on hydrography, chemistry and morbidity). J Biol Board Can 1:219–467
Griffith AW, Shumway SE, Gobler CJ (2019) Differential mortality of North Atlantic bivalve molluscs during harmful algal blooms caused by the dinoflagellate, Cochlodinium (a.k.a Margalefidinium) polykrikoides. Estuarines and Coasts 42:190–203
Guseva K, Feudel U (2020) Numerical modelling of the effect of intermittent upwelling events on plankton blooms. J R Soc Interface 14(20190):880
Hastings A, Abbott KC, Cuddington K et al (2018) Transient phenomena in ecology. Science 361:990
Hillebrand H, Brey T, Gutt J et al (2018) Climate change: warming impacts on marine Biodiversity. In: Salomon M, Markus T (eds) Handbook on marine environment protection. Springer, Cham
Hirsche HJ (2013) Long-term experiments on lifespan, reproductive activity and timing of reproduction in the Arctic copepod Calanus hyperboreus. Mar Biol 160:2469–2481
Hjerne O, Hajdu S, Larsson U et al (2019) Climate driven changes in timing, composition and magnitude of the Baltic Sea phytoplankton spring bloom. Front Mar Sci 6:482
Holling CS (1959a) Some characteristics of simple types of predation and parasitism. Can Entomol 91:385–398
Holling CS (1959b) The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Can Entomol 91(5):293–320
Izhikevich EM (2007) Dynamical systems in neuroscience. MIT Press, Cambridge
Klais R, Lehtiniemi M, Rubene G et al (2016) Spatial and temporal variability of zooplankton in a temperate semi-enclosed sea: implications for monitoring design and long-term studies. J Plankton Res 38(3):652–661
Krönke J, Wunderling N, Winkelmann R et al (2020) Dynamics of tipping cascades on complex networks. Phys Rev E 101(042):311
Largier JL (2020) Upwelling bays: how coastal upwelling controls circulation, habitat, and productivity in bays. Ann Rev Mar Sci 12:415–447
Lenton TM (2020) Tipping positive change. Proc R Soc B 375(20190):123
Lenton TM, Held H, Kriegler E et al (2008) Tipping elements in the Earth ’ s climate system. Proc Natl Acad Sci USA 105(6):1786–1793
Lewandowska AM, Boyce DG, Hofmann M et al (2014) Effects of sea surface warming on marine plankton. Ecol Lett 17(5):614–623
Lewandowska AM, Striebel M, Feudel U et al (2015) The importance of phytoplankton trait variability in spring bloom formation. ICES J Mar Sci 72:1908–1915
Mitry J, McCarthy M, Kopell N et al (2013) Excitable neurons, firing threshold manifolds and canards. J Math Neurosci 3(12)
Morel F, Rueter J, Price N (1991) Iron nutrition of phytoplankton and its possible importance in the ecology of ocean regions with high nutrient and low biomass. Oceanography 4:56–61
Morozov A, Abbott KC, Cuddington K et al (2020) Long transients in ecology: Theory and applications. Phys Life Rev 32:1–40
O’Keeffe PE, Wieczorek S (2020) Tipping phenomena and points of no return in ecosystems: Beyond Classical Bifurcations. SIAM J Appl Dyn Syst 19(4):2371–2402
O’Sullivan E, Mulchrone K, Wieczorek S (2022) Rate-induced tipping to metastable zombie fires. arXiv preprint arXiv:2210.02376
Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37:637–669
Perryman C, Wieczorek S (2014) Adapting to a changing environment: non-obvious thresholds in multi-scale systems. Proc R Soc A Math Phys Eng Sci 470(2170):20140,226
Pinek L, Mansour I, Lakovic M et al (2020) Rate of environmental change across scales in ecology. Biol Rev 95:1798–1811
Richards KJ (2017) Viral infections of oceanic plankton blooms. J Theor Biol 412:27–35
Riley G (1946) Factors controlling phytoplankton populations on Georges Bank. J Mar Res 6:54–73
Riley G, Stommel H, Bumpus D (1949) Quantitative ecology of the plankton of the western North Atlantic. Bulletin of the Bingham Oceanographic Collection 12:1–169
Rumyantseva A, Henson S, Martin A et al (2019) Phytoplankton spring bloom initiation: The impact of atmospheric forcing and light in the temperate North Atlantic Ocean. Prog Oceanogr 178(102):202
Scheffer M, Carpenter S, Foley JA et al (2001) Catastrophic shifts in ecosystems. Nature 413:591–596
Simpson JH, Sharples J (2012) Introduction to the physical and biological oceanography of shelf seas. Cambridge University Press, New York
Smith SJ, Edmonds J, Hartin CA et al (2015) Near-term acceleration in the rate of temperature change. Nat Clim Chang 5:333–336
Sommer U, Lewandowska AM (2011) Climate change and the phytoplankton spring bloom: warming and overwintering zooplankton have similar effects on phytoplankton. Glob Chang Biol 17(1):154–162
Sommer U, Adrian R, De Senerpont Domis L et al (2012) Beyond the plankton ecology group (PEG) model: Mechanisms driving plankton succession. Annu Rev Ecol Evol Syst 43:429–448
Sverdrup H (1953) On conditions for the vernal blooming of phytoplankton. ICES J Mar Sci 18(3):287–295
Szmolyan P, Wechselberger M (2001) Canards in \(\mathbb{R} ^3\). J Differ Equ 177:419–453
Trombetta T, Vidussi F, Mas S et al (2019) Water temperature drives phytoplakton blooms in coastal waters. PLoS ONE 14(4):e0214933
Truscott J, Brindley J (1994) Ocean plankton populations as excitable media. Bull Math Biol 56(5):981–998
Uye S (1986) Impact of copepod grazing on the red tide flagellate. Mar Biol 92:35–43
Vanselow A, Wieczorek S, Feudel U (2019) When very slow is too fast - collapse of a predator-prey system. J Theor Biol 479:64–72
Velo-Suárez L, Brosnahan ML, Anderson DM et al (2013) A quantitative assessment of the role of the parasite Amoebophrya in the termination of Alexandrium fundyense blooms with a small coastal embayment. PLoS ONE 8(12):e81,150
Wake G (1991) Picoplankton and marine food chain dynamics in a variable mixed layer: a reaction-diffusion model. Ecol Model 57:193
Walther GR, Post E, Convey P et al (2002) Ecological responses to recent climate change. Nature 416:389–395
Wechselberger M, Mitry J, Rinzel J (2013) Canard theory and excitability. In: Kloeden PE, Pötzsche C (eds) Nonautonomous Dynamical Systems in the Life Sciences. Springer, pp 89–133
Wieczorek S, Ashwin P, Luke C et al (2011) Excitability in ramped systems: the compost-bomb instability. Proc R Soc Lond Ser A 467:1243–1269
Wieczorek S, Xie C, Ashwin P (2021) Rate-induced tipping: Thresholds, edge states and connecting orbits. arXiv preprint arXiv:2111.15497
Winder M, Berger SA, Lewandowska A et al (2012) Spring phenological responses of marine and freshwater plankton to changing temperature and light conditions. Mar Biol 159:2491–2501
Funding
Open Access funding enabled and organized by Projekt DEAL. This work was supported by the DAAD-DST grant No. 57452783. U.F. acknowledges support by the European Union’s Horizon 2020 Research and Innovation program under the Marie Sklodowska-Curie Action Innovative Training Networks Grant Agreement No. 956170 (CriticalEarth). P.P. acknowledges support from the DST-DAAD project (Ref. No.: DST/INT/DAAD/P-15/2019).
Author information
Authors and Affiliations
Contributions
AV, SW, UF conceived and designed the research, AV, SW, PP set up the model and the simulations. At last, AV analyzed the model results. AV, LH, SW and UF designed the storyline and the concept of the manuscript. All authors read and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval
Not applicable.
Consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interest
The authors declare no competing interests.
Appendices
Appendix A. The Truscott-Brindley model
A.1 Non-dimensionalization
The phytoplankton(P)-zooplankton(Z) model developed by Truscott and Brindley (1994) is given by the following two equations:
with the growth rate of the phytoplankton population r and its carrying capacity K. The attack rate of the zooplankton \(R_m\), its half-saturation constant \(\alpha\), conversion efficiency \(\gamma\) and mortality rate \(\mu\) complete the model.
In the following, we introduce the non-dimensional variables \(\hat{P} = \frac{P}{K}\), \(\hat{Z} = \frac{Z}{K}\) and the non-dimensional time \(\hat{t} = \tau \cdot R_m\). Replacing P, Z and \(\tau\) by their non-dimensional equivalent, we can rewrite the TB-model (A1)–(A2) as follows:
Setting \(\delta = \frac{\alpha }{K}\), \(\beta = \frac{r}{R_m}\) and \(\omega = \frac{\mu }{\gamma R_m}\), the non-dimensional TB-model (A3)–(A4) becomes:
The parameter \(\gamma\) quantifies the separation between the time scale of the phytoplankton and the time scale of the zooplankton. For this reason, we call \(\gamma\) the time scale parameter. If \(\gamma < 1\), the zooplankton population reproduces slower than the phytoplankton and vice versa. For \(\gamma = 1\), there exists no time scale separation.
We further assume that the growth rate of the phytoplankton r becomes time-dependent: \(r=r(\tau )\) and increases at a constant rate \(\eta\):
Using the dimensionless time \(\hat{t} = \tau \cdot R_m\), Eq. (A7) becomes:
with \(\beta = \frac{r}{R_m}\), we can write:
Equations (A5), (A6) and (A9) show the three-dimensional plankton model in fast time \(\hat{t}\). We now assume that the temporal evolution of r and \(\beta\) happens on the same time scale as the temporal evolution of zooplankton. Therefore we use the time scale separation parameter to scale the rate of change of r: \(\eta =\gamma \upsilon\). We can write the whole model system in slow time \(\tilde{t} = \gamma \hat{t}\) using \(\frac{\textrm{d}}{\textrm{d}\hat{t}} = \frac{\textrm{d}}{\textrm{d}\tilde{t}} \frac{\textrm{d}\tilde{t}}{\textrm{d}\hat{t}} = \frac{\textrm{d}}{\textrm{d}\tilde{t}}\gamma\):
A.2 Long-term states and linear stability
To find the long-term states of the TB-model (A1)–(A2), we set the dynamics of the phytoplankton and zooplankton population equal to zero:
Solving Eqs. (A13)–(A14) with respect to P and Z leads to the following three long-term states \(e_1\)–\(e_3(r)\) (A15)–(A17)
The equilibrium \(e_1\) describes the situation in which phytoplankton and zooplankton are extinct. In \(e_2\), the phytoplankton possesses the maximum density that the environment can “carry” (K) while the zooplankton is extinct. In the third equilibrium \(e_3(r)\), phytoplankton and zooplankton coexist at densities unequal zero. Notice that only the equilibrium \(e_3(r)\) depends on the growth rate of the phytoplankton r. Within the interval \(r \in [0.2\; 2]\), it represents the unique stable state of the TB-model (Fig. 7). Hence for values of the growth rate between \(r_{min}\) and \(r_{max}\), the phytoplankton and zooplankton coexist in a stable long-term equilibrium when \(t \rightarrow \infty\).
Linear stability of the three equilibria \(e_1\)–\(e_3(r)\) (A15)–(A17) of the TB-model (A1)–(A2) depending on the growth rate \(r \in [0.2\; 2]\). Solid/dotted line denotes the real parts of the eigenvalues \(\lambda _{1,2} < 0\) (stable) respectively \(\lambda _{1,2} > 0\) (unstable). For the other parameters see Table 1
A.3 The critical manifold
The extended Truscott-Brindley model with time-dependent growth rate r can be written in slow time \(t = \gamma \tau\) (\(\kappa = \mu /\gamma\), \(\upsilon = \eta /\gamma\)) (Truscott and Brindley 1994):
The TB-model with time-dependent growth rate r (A18)–(A20) is determined, as the original TB-model, by a fast and a slow time scale. In fact, the dynamics of such slow-fast systems is primarily slow with short interruption by fast motion. The critical manifold S approximates the slow motion and is therefore a useful tool to obtain a first impression of a part of the full dynamics.
The critical manifold S (A21) can be further written as:
with the two folds \(F_{1,2}(r) = (P_{F_{1,2}}, Z_{F_{1,2}})\) given by the following equations
Appendix B. The canard trajectory
When studying the slow-fast dynamics of the extended TB-model (B25)–(B27) with time-dependent growth rate r for different rates \(\upsilon\), we find a plankton bloom when the growth rate r increases faster than \(\upsilon \approx 0.003\) (\(\eta \approx 0.0015\)) (Fig. 1, see Sec. “Rate-induced tipping triggers plankton bloom” for more details)
To simplify notations, we use from now on the following abbreviations. Note that q is denoting f and g and k stands for P, Z and r.
The bloom forms when (red) trajectories cross the fold \(F_1(r)\) of the critical manifold S and run away in the fast P-direction (Fig. 4). For this reason, we study the fast flow close to the fold \(F_1(r)\) on the critical manifold S. To study the dynamics on S, we set the fast dynamics \(\frac{\textrm{d}P}{\textrm{d}t}\) equal to zero (\(\gamma \rightarrow 0\)) which gives the reduced system:
Differentiating the algebraic constraint (B30) with respect to the slow time t leads to:
The Eq. (B35) describes the fast flow P on the critical manifold S. Replacing the constrain Eq. (B30) by Eq. (B35), we obtain:
Using \(Z = h(P,r)\) (A22), we project the flow onto the critical manifold S. The so-called projected-reduced system is given by:
With \(\Lambda (P,r,\upsilon ) = f^S_Z(P,r)\; g^S(P) + f^S_{r}(P,r) \upsilon\) and \(f^S_P(P,r) = \frac{\partial f(P,r)}{\partial P}\bigg \vert _{S}\) we can write the projected-reduced system as:
For completeness, we write the projected-reduced system in full terms
If \(f^S_P(P,r) = 0\) (B39), the flow on the critical manifold S (11) is not defined — it goes to infinity in finite time. Unfortunately, these points for which \(f^S_P(P,r) = 0\) are the fold points of the critical manifold S (11). Since the bloom formation occurs close to the fold, we have to enable the analysis near the fold \(F_1(r)\). This can be achieved by using the scaling known as desingularization:
which preserves the direction of time on the stable part of the critical manifold S (red, Fig. 9), but reverses it on the unstable part (blue, Fig. 9). The desingularized system is given by:
respectively in full terms,
Interestingly, setting \(\frac{\textrm{d}r}{\textrm{d}\hat{t}}=0\):
and rearranging Eq. (B50) with respect to P gives the \(P_{F_1}\)-values of the fold \(F_1(r)\). When we further set \(P_{F_1}\) into Eq. (B48), we can analyze the fast flow \(\frac{\textrm{d}P}{\textrm{d}\hat{t}}\) at the fold \(F_1(r)\) (Fig. 8). For \(\upsilon = 0.002\), there exist a special solution at which \(\frac{\textrm{d}P}{\textrm{d}\hat{t}} = 0\) (white circle) within \(r \in [0.2\; 0.6]\). This special solution is an equilibrium of the desingularized system (B48)–(B49) and it marks the boundary between fold points that attract trajectories (\(\frac{\textrm{d}P}{\textrm{d}\hat{t}} < 0\), red) and fold points that repel them (\(\frac{\textrm{d}P}{\textrm{d}\hat{t}} > 0\), blue). Such special points are called folded singularities (Szmolyan and Wechselberger 2001). For \(\upsilon = 0.004\), we found no folded singularity in the interval \(r \in [0.2\; 0.6]\).
The fast flow \(\frac{\textrm{d}P}{\textrm{d}\hat{t}}\) (B48) (black solid line) on the fold \(F_1(r)\) (A23)–(A24) is shown for two different rates \(\upsilon = 0.002\) and \(\upsilon = 0.004\). Fold Points with \(\frac{\textrm{d}P}{\textrm{d}\hat{t}} <0\) attract trajectories (red) while fold points with \(\frac{\textrm{d}P}{\textrm{d}\hat{t}} >0\) are repelling (blue). The fast flow \(\frac{\textrm{d}P}{\textrm{d}\hat{t}}\) remains only finite at the point FS. For other parameters see Table 1
Studying the linear stability of the folded singularity at \(F_1(r)\) gives one negative eigenvalue (\(\lambda <0\)) and one positive eigenvalue (\(\lambda > 0\)). Hence, it is called a folded saddle singularity FS.
And indeed, solutions that cross the fold via the folded saddle singularity show a kind of boundary behavior: they can cross the fold with finite speed and move away from the quasi-static state (bloom) but slowly (no bloom). Such trajectories are called (singular) canards (Fig. 9). In the desingularized system (B48)–(B49) exist besides the folded saddle FS (\(\lambda _1 < 0,\lambda _2 >0\)) folded focus FF (Re(\(\lambda _{1,2}\)) \(< 0\)) at the second fold \(F_2(r)\). In the case of a folded saddle FS, the singular canard trajectory is given by its stable manifolds.
Dynamics of the desingularized system (B48)–(B49) for two different rates \(\upsilon = 0.002\) (a, c, e) and \(\upsilon = 0.004\). (b, d, f). The singular canard (blue trajectory) is given by the stable manifolds of the folded saddle FS (white filled circle). It represents the threshold separating (green) tracking from (red) tipping trajectories. a, c The green trajectory starts above the singular canard and tracks — no bloom emerges (e). b, d The red trajectory starts below the canard, is attracted by the fold \(F_1(p,r)\), crosses it and approaches the stable folded focus FF (black filled circle) — a plankton bloom occurs (f)
Appendix C. More than one recurring bloom — parameter studies
The occurrence of more than one plankton bloom while thes growth rate r increases in time depends i.a. on the attack rate \(R_m\) and the mortality rate \(\mu\) of the zooplankton as well as the maximum growth rate \(r_{max}\) of the phytoplankton. Of course, \(r_{max} = 20\) is far from any realistic approach; nevertheless, we choose this value to evaluate how many recurring blooms can be theoretically observed for extreme high maximum growth rates.
In the following, we outline our procedure for finding the number of blooms exemplary for \(r_{max} = 0.8\) (Fig. 10b, d) and different values of the attack rate \(R_m\) (Fig. 10a). (i) At first, we compute the parameter interval of \(R_m\) for which \(e(\upsilon t)\) represents the unique stable state. (ii) Then, we evaluate the rate \(\upsilon\) for which the r-coordinate of the folded saddle FS is equal to \(r_{max}\). (iii) We simulate the three-dimensional TB-model (B25)–(B27) for the corresponding attack rates \(R_m\) and values of the rate \(\upsilon\) (\(P(0) = e_{3P}\), \(Z(0) = e_{3Z}\) (see Table 1 for other parameters). (vi) If the maximum phytoplankton density \(P_{max}\) exceeds the threshold \(P_b = 60\), we increase the parameter \(N_b\) which displays the number of blooms, by one. (v) We check randomly if we have detected the correct number of blooms by visible examination of the corresponding trajectory.
Obviously, decreasing the attack rate \(R_m\) and increasing the mortality rate \(\mu\) of the zooplankton cause an increase of the number of recurring blooms (Fig. 10). Hence, decreasing the predation pressure on the phytoplankton allows for more and more blooms while r increases linearly in time.
The number of recurring blooms \(P_{max} > P_b = 60\) depending on the attack rate \(R_m\) (a, c) and the mortality rate \(\mu\) of the zooplankton (b, d) as well as the maximum growth rate \(r_{max}\). Decreasing attack rate \(R_m\) and increasing mortality rate \(\mu\) promote the occurrence of more and more plankton blooms during the growth rate r increases linearly in time. Parameter: \(P(0) = e_P\), \(Z(0) = e_Z\); others see Table 1
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Vanselow, A., Halekotte, L., Pal, P. et al. Rate-induced tipping can trigger plankton blooms. Theor Ecol 17, 89–105 (2024). https://doi.org/10.1007/s12080-024-00577-5
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s12080-024-00577-5