Data-driven dynamics of phytoplankton blooms in a reaction–diffusion NPZ model

Abstract

The causes of seasonal phytoplankton blooms in the ocean are a debated topic. The disturbance-recovery hypothesis posits that blooms are initiated when seasonally changing light and mixing conditions disrupt attractive equilibrium states in the predator–prey relationship between zooplankton and phytoplankton, leading to an accumulation in phytoplankton biomass. This paper follows up on this notion with a nutrient-phytoplankton-zooplankton (NPZ) model incorporating diffusion and depth-dependent coefficients in which we explore the nature of these attractive states. The reaction–diffusion NPZ model, incorporated with seasonal solar radiation and mixed layer depth data, simulates blooms with better timing than a related ordinary differential equations model but still leaves room for further modeling and improvement. The results of the data-driven, reaction–diffusion model aid in understanding the driving forces of the simulated phytoplankton blooms. The results also reveal a significant influence on the vertical structure of plankton concentration by the attenuation of light with depth in the ocean. Models that accurately simulate blooms tend to share characteristics of the reaction–diffusion model in this paper. The dynamics described in this paper could be a key mechanism that seasonally disrupts the ecological balance between phytoplankton and zooplankton and facilitates a high-latitude marine phytoplankton bloom.

INTRODUCTION

A variety of hypotheses about the initiation of spring phytoplankton blooms in the North Atlantic have been published since Gran and Braarud’s Bay of Fundy Project (Gran and Braarud, 1935), the conclusion of which eventually became known as the critical depth hypothesis (CDH). Revisiting Gran and Braarud’s study, Sverdrup (1953) proposed that, assuming a well-mixed surface layer and constant background losses due to respiration, phytoplankton blooms in the North Atlantic are caused primarily by alleviating light limitation through the shoaling of the surface mixed layer depth above a critical depth, where the net phytoplankton production is equal to net respiration. Since then, alternative theories have been offered. The critical turbulence hypothesis (CTH) (Huisman et al., 1999) suggests that if turbulent mixing rates in the surface layer are less than growth rates, then growth will overcome losses due to mixing and a bloom will form. Behrenfeld put forward the disturbance-recovery hypothesis (DRH) (Behrenfeld, 2010; Behrenfeld and Boss, 2014; Behrenfeld and Boss, 2018) that suggests that background losses to phytoplankton populations are not constant. A major feature of DRH is that a bloom initiates in early winter as the deepening mixed layer reduces contact rates between the predators (zooplankton) and the prey (phytoplankton) through dilution, allowing the phytoplankton population to bloom by escaping grazing pressure. During this period, however, phytoplankton concentrations may still be decreasing due to dilution, whereas depth-integrated stocks are increasing. Behrenfeld and Boss (2014) argue that a paradigm shift is needed as to the way blooms are explained. Their analysis suggests that the primary cause of marine phytoplankton blooms is not by processes that affect the rates of cell division, as in the CDH and CTH, but rather by processes that disrupt the ecological equilibria between predator (zooplankton) and prey (phytoplankton). Given the numerous complexities of the world’s oceans, it is plausible that many, if not all, of these theories have some merit. It is not necessarily true that one of these theories is entirely correct, whereas the others are completely wrong. However, in this manuscript, we examine the perspective of DRH, as it relates to the formation of phytoplankton blooms. Although a bloom is understood as an excessive accumulation of phytoplankton biomass, there is no consensus on what defines a bloom beginning or ending. When we refer to bloom initiation, we do so in a manner consistent with a “rate of change” method of identifying the start date of a phytoplankton bloom (Brody et al., 2013). The rate of change methods estimate the time of bloom initiation at the point of most rapid increase on a chlorophyll time series. In our model, we refer to bloom initiation as the point of most rapid increase of our simulated phytoplankton concentration time series, although we do not calculate this exact point.

We consider a nutrient-phytoplankton-zooplankton (NPZ) model, a system of coupled, nonlinear differential equations, to simulate phytoplankton blooms as transitions between ecological equilibria. Since the quantities of solar radiation, mixed layer thickness and turbulent mixing have been used to understand theoretical causes of phytoplankton blooms, these quantities are incorporated into the model. An autonomous form of an NPZ model (with no time dependence in its coefficients) can have multiple equilibrium states. Time-dependent coefficients can cause the attracting equilibrium state of the system to periodically change from one state to another. Analogously, environmental conditions can influence the balance of a marine ecosystem. Seasonally changing conditions sometimes cause drastic transitions between one balanced state and another. A phytoplankton bloom may be such a transition.

To begin our narrative of blooms, we incorporate a light response coefficient |$\rho$|⁠, which depends on solar radiation |$I$|⁠, into the ordinary differential equations model of Franks et al. (1986). If |$I$|⁠, and consequently |$\rho$|⁠, is independent in time, the resulting autonomous model has three equilibrium states determined by Busenberg et al. (1990): plankton free, phytoplankton only and co-existing phytoplankton and zooplankton. Busenberg et al. (1990) also show that only one of these states can attract at one time, and if a state attracts, the other two must repel. In this paper, the effect of changing the value of |$\rho$| on the dynamics of the system is examined, with attention to conditions that facilitate a model phytoplankton bloom. Based on the value of |$\rho$|⁠, it is determined analytically which equilibrium state of the system is attracting and how that can change to another state, which is a key tenant of the DRH hypothesis. Using North Atlantic solar radiation data to describe the quantity |$I(t)$|⁠, which changes |$\rho$|⁠, the resulting phytoplankton predictions are compared with North Atlantic chlorophyll data. This produces a solution that exhibits blooming behavior with the peaks of the blooms occurring shortly after the new year. Blooms in the North Atlantic peak much later in the spring; to address this, the model is enhanced to a reaction diffusion system. However, analysis of the dynamics of the autonomous ordinary differential equations model illuminates the process of a phytoplankton bloom in the reaction–diffusion model.

Other studies have analyzed the dynamics of ordinary differential equation plankton ecology models to give explanations of phytoplankton blooms. Some use state variables for nutrients and phytoplankton only to study blooms that are primarily bottom-up controlled by nutrients rather than top-down controlled by higher trophic levels. Huppert et al. (2002) construct a nutrient-phytoplankton model where blooms are triggered when nutrients exceed a threshold and blooms are further controlled by initial values of nutrients and phytoplankton. Many studies follow up on the autonomous phytoplankton-zooplankton model of Truscott and Brindley (1994), commonly referred to as the T-B model, in which the phytoplankton grow logistically. Blooms in the T-B model occur following a sufficiently low initial concentration of zooplankton. Freund et al. (2006), Gao et al. (2009), and Zhao and Yan (2016) analyze the effect of seasonal forcing on the blooming behavior of the T-B model. Luo (2013) defines the phytoplankton maximum growth rate and carrying capacity of the T-B model as functions of time and nutrient concentration. Luo (2013) observes that blooms are sensitive to phytoplankton carrying capacity, which depends on nutrient concentration. Chatterjee and Pal (2011) and Zhang and Wang (2012) each describe annual cycles of phytoplankton populations, including blooms, as oscillations around a coexisting equilibrium in an NPZ model. The present paper identifies the effect of the phytoplankton-only equilibrium on a bloom in an ordinary differential equations NPZ model when sunlight intensity increases. This paper also identifies an analogous effect that the phytoplankton-only equilibrium has on a depth-dependent reaction–diffusion (partial differential equation) NPZ model as sunlight intensity increases and vertical diffusion decreases.

To obtain a model which simulates blooms that reflect actual phytoplankton blooms in the ocean, we incorporate diffusion (used as a proxy for turbulent mixing), mixed layer depth and attenuation of light with depth into the ordinary differential equations model, resulting in a reaction–diffusion system. Using the techniques of Edwards et al. (2000), the equilibrium states of the reaction–diffusion model and the conditions under which they attract are determined computationally. Solar radiation and mixed layer depth data from the North Atlantic Ocean are incorporated into the reaction–diffusion model. The effect of the seasonal data on the dynamics of the NPZ reaction–diffusion system is analyzed graphically. The improvements of bloom timing resulting from the incorporation of diffusion, mixed layer depth and attenuation of light with depth are discussed. This analysis suggests a mechanism for facilitating marine phytoplankton blooms that is congruent with the concept of understanding phytoplankton blooms as disruptions in predator–prey equilibria suggested by the DRH (Behrenfeld and Boss, 2014). Limitations of the reaction–diffusion model are also discussed. Finally, similarities between the results of this paper and others are discussed that suggest similar dynamical explanations of phytoplankton blooms.

METHODS

Construction of the ordinary differential equations model

The ordinary differential equations model from Cowall et al. (2019), which is similar to the model of Franks et al. (1986) with the addition of the parameter |$\rho$|⁠, is as follows:

$$\begin{eqnarray} {\displaystyle \begin{array}{l}\frac{dP}{dt}\!=\!\rho (I){\mu}_{max}\frac{N}{k_N\,+\,N}P-{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_{P}P\\[4pt] {}\frac{dZ}{dt}\!=\!\alpha{h}_{max}\left(1\!-\!{e}^{- aP}\right)Z\!-\!{m}_ZZ\\[4pt] {}\frac{dN}{dt}\!=\!-\rho (I){\mu}_{max}\frac{N}{k_N\,+\,N}P\!+\!\left(1\!-\!\alpha \right){h}_{max}\left(1\!-\!{e}^{- aP}\right)Z\\ \qquad\quad +\,{m}_PP+{m}_ZZ.\end{array}} \end{eqnarray}$$

(1)

There is an equation for each state variable: phytoplankton (⁠|$P$|⁠), zooplankton (⁠|$Z$|⁠) and nutrients (⁠|$N$|⁠). Assuming nitrogen is the limiting nutrient of system, the units of each state variable is mmol N m⁻³. Starting with the model of Franks et al. (1986), we incorporate a phytoplankton light response coefficient

$$\begin{equation} \rho (I)=\frac{\Gamma I}{\sqrt{\mu_{max}^2+{\Gamma}^2{I}^2}} \end{equation}$$

(2)

where |$I(t)$| is the intensity of the solar radiation at the surface and |$0<\rho <1$|⁠. Table I gives a description of the parameters and the choices of values indicated by Gentleman and Neuheimer (2008) with the exception of |${h}_{max}$|⁠, which is reduced from 1.5 d⁻¹ to 0.5 d⁻¹, the value chosen by Busenberg et al. (1990). This reduction is chosen to increase the period of the periodic solutions to 1 year, which we would expect from an ecosystem that is dependent upon annual cycles of solar radiation intensity. These values will be used throughout the paper unless otherwise specified. Since the right-hand side of ordinary differential equation system (1) sums to zero, then |$N+P+Z={N}_T$| for constant |${N}_T$|⁠. The system (1) can be reduced to the following Equation (3).

$$\begin{equation} {\displaystyle \begin{array}{l}\frac{dP}{dt}=\rho{\mu}_{max}\frac{N_T\,-\,P\,-\,Z}{k_N\,+\,{N}_T\,-\,P\,-\,Z}P-{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_PP\\{}\frac{dZ}{dt}=\alpha{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_ZZ\end{array}} \end{equation}$$

(3)

Table I

Parameter descriptions and values for Equation ( 1)

Parameter	Description	Value	Units
|$\alpha$|	Assimilation efficiency	0.7	–
|${\mu}_{max}$|	Phytoplankton maximum growth rate	2	d−1
|$a$|	Ivlev decay constant	1	(mmol N m−3)−1
|${h}_{max}$|	Maximum grazing rate	0.5	d−1
|${k}_N$|	Michaelis–Menten half-saturation constant	1	mmol N m−3
|${m}_P$|	Phytoplankton mortality rate	0.1	d−1
|${m}_Z$|	Zooplankton mortality rate	0.2	d−1
|${N}_T$|	Total nitrogen in the system	2	mmol N m−3
|$\Gamma$|	Initial photosynthetic slope	0.055	(E m−2)−1d−1

Values from Gentleman and Neuheimer (2008) except for the value of |${h}_{max}$|⁠, which is chosen by Busenberg et al. (1990), and the value of Γ, which is chosen to remain in the range indicated by Kuhn et al. (2015), Fennel et al. (2006) and Schartau and Oschlies (2003).

Data-driven blooms in the ordinary differential equations model

In system (3), solar radiation data are incorporated to describe the quantity |$I(t)$|⁠. The satellite data photosynthetically active radiation (PAR) and chl-ocx (chlorophyll) collected by the NASA MODIS-Aqua satellite are used. The PAR and chlorophyll data are available at https://oceancolor.gsfc.nasa.gov/. The region 25–35°W, 40–45°N in the North Atlantic Ocean is investigated. Time series of the horizontally averaged PAR and chlorophyll data from 2009 through 2015 are used for this region. The predictions of model concentrations of phytoplankton from Equation (3) using PAR data for |$I(t)$| are compared with a time series of the chlorophyll data. Figure 1 shows the predictions of the system (3) on the North Atlantic region. Sinusoidal best-fit approximations for PAR, shown in Fig. 1, are used to describe |$I(t)$| to avoid noise in the solution. The resulting phytoplankton time series is similar to that of a result from Behrenfeld and Boss (2018).

Fig. 1

Time series of solution of ordinary differential equation model system (3) with I(t) described by PAR data. All data averaged over the region 40–45°N, 25–35°W in the North Atlantic. Parameter values from Table I. Top: time series of PAR data with sinusoidal best-fit curve. Middle: time series of model phytoplankton concentration predictions compared with chlorophyll data. Bottom: time series of surface model phytoplankton and zooplankton predictions.

Seasonal variation of solar radiation is the only driving force behind the simulated phytoplankton bloom described by Equation (3) with |$\rho$| defined as in (2). The bloom cycle begins in the summer when solar radiation is at its maximum. Decreasing sunlight in autumn and early winter diminish the phytoplankton. This then depletes the food supply for the zooplankton, causing their population to decrease. Low grazing pressure and increasing sunlight after the winter solstice facilitate rapid phytoplankton growth. It is during this time that the phytoplankton concentration achieves its point of most rapid growth rate: the beginning of the bloom according to the “rate of change” method of identifying bloom initiation. The bloom is essentially the phytoplankton overcompensating from their low winter levels. This overcompensation is aided by the fact that the phytoplankton can grow faster than the zooplankton. The delayed recovery by the zooplankton brings the phytoplankton back down to summer levels and removes the bloom.

The description of this cycle is neither a complete nor an accurate narrative of a bloom cycle in the ocean. This is due in part to sunlight being the only time-dependent driver of the model. One observable problem in the results in Fig. 1 is that the model phytoplankton begin their rapid increase in concentration before the winter solstice. The bottom plot of Fig. 1 shows the beginning of the increase in phytoplankton concentration slightly preceding the lowest zooplankton concentration of the year. This suggests that the dominant initiator of the model bloom is a critically low level of grazing pressure. This is not evident in the data, as the chlorophyll concentration begins its rapid increase later in the spring. Evidently, this is a limitation of the model Equation (3) because the system is driven only by varying seasonal sunlight intensity. Also, since the magnitude of the annual conditions influencing blooms is not the same every year, then the magnitude of the observed chorophyll peaks varies from year to year. This feature is not captured in our model results because the input from sunlight intensity is periodic. Despite these limitations, a deeper mathematical description of the model bloom from Equation (3) is helpful for interpreting our yet to be described reaction–diffusion NPZ model that incorporates depth-dependent factors necessary for describing the bloom cycle.

Mathematical underpinnings of the bloom

Busenberg et al. (1990) identified the three equilibria of the coupled Equation (1) without the parameter ρ. The equivalent equilibria when the parameter ρ is included are:

• a “plankton-free” equilibrium |$(P,Z,N)=(0,0,{N}_T),$|

• a “phytoplankton-only” equilibrium |$(P,Z,N)=(\overline{P},0,{N}_T-\overline{P})$|⁠, and

• a “coexisting” equilibrium |$(P,Z,N)=({P}^{\ast },{Z}^{\ast },{N}_T-{P}^{\ast }-{Z}^{\ast}),$|

where

$$\begin{eqnarray} \overline{P}\!\!\!\!&={N}_T-\frac{m_P{k}_N}{\rho{\mu}_{max}-{m}_P}>0,\nonumber\\{P}^{\ast}\!\!\!\!&=-\frac{1}{a}\ln \left(1-\frac{m_Z}{\alpha{h}_{max}}\right)>0, \end{eqnarray}$$

(4)

and |${Z}^{\ast }$| being the unique solution, satisfying |$0<{Z}^{\ast }<{N}_T-{P}^{\ast }$|⁠, of the quadratic equation

$$\begin{equation} {\displaystyle \begin{array}{l}{Z}^2-Z\left[\frac{\alpha \left(\rho{\mu}_{max}-{m}_P\right)}{m_Z}{P}^{\ast }+{k}_N+{N}_T-{P}^{\ast}\right]\\{}\quad+{P}^{\ast}\frac{\alpha }{m_Z}\left[\left(\rho{\mu}_{max}-{m}_P\right)\left({N}_T-{P}^{\ast}\right)-{k}_N{m}_P\right]=0.\end{array}} \end{equation}$$

(5)

|${P}^{\ast }$| and |${Z}^{\ast }$| are the phytoplankton and zooplankton concentrations, respectively, when the ecosystem is balanced. |$\overline{P}$| is the expected phytoplankton concentration in the absence of zooplankton. Although we do not expect a marine ecosystem to be devoid of zooplankton, the phytoplankton-only equilibrium |$(\overline{P},0)$| significantly influences the dynamics of the model bloom.

We use three results from Busenberg et al. (1990) to give an alternative description of the bloom process in the previous section. These results hold whenever the coexisting equilibrium mathematically exists, which depends on the values of the parameters in Table I. First, the coexisting equilibrium point is asymptotically stable (solutions that start near the point move toward the point in time) and |$\overline{P}>{P}^{\ast }$|⁠, i.e. the concentration of phytoplankton would be higher if they do not have to share the ocean with zooplankton. Also, the phytoplankton-only equilibrium is a saddle point with stable manifold |$\{(P,Z):Z=0,0<P\le{N}_T\}$|⁠. This means trajectories that start close to the stable manifold (horizontal axis, i.e. phytoplankton-only populations) will move toward the phytoplankton-only equilibrium but then turn away from that toward the coexisting equilibrium. Figure 2 illustrates a saddle point in the context of our model. When sunlight decreases from its summer maximum, the phytoplankton population also decreases. With their food source diminishing, the zooplankton population declines more rapidly than the phytoplankton due to their higher mortality rate. Thus, the system temporarily approaches a phytoplankton-only equilibrium, a saddle point, (implying |$P$| approaches |$\overline{P}$|⁠). Since |$\overline{P}>{P}^{\ast }$| then the phytoplankton concentration will be higher than |${P}^{\ast }$| when the system is close to the phytoplankton-only equilibrium. At this point, the bloom is at its peak. Then, as the zooplankton recover, the system moves away from the phytoplankton-only equilibrium, ending the bloom (see Fig. 2). Similar behavior is observed in the reaction–diffusion model. However, in addition to sunlight, the stability properties of the reaction–diffusion system are also controlled by mixed layer depth and diffusion. Thus, all three of these physical drivers will be used to describe the bloom cycle simulated by the reaction–diffusion model.

Construction of the reaction–diffusion model

Predictions of Equation (1) using |$\rho$| defined as in (2) produce a rapid increase in phytoplankton biomass that starts, approximately, just before the winter solstice, see Fig. 1. This can be interpreted as the beginning of the bloom by a rate-of-change method, although the DRH posits that the blooming process starts much sooner than the winter solstice. This biomass increase generally happens 2–4 months later in the North Atlantic. This issue is addressed by incorporating physical driving forces, drivers that have previously been used to explain phytoplankton blooms: diffusion and mixed layer depth. Despite the shortcomings of the ordinary differential equations model (1), the study of its dynamical properties does inform the study of the expansion of (1) to a spatial, reaction–diffusion system (6) with depth |$\xi$|⁠.

$$\begin{eqnarray} {\displaystyle \begin{array}{l}{\partial}_tP=\rho{e}^{-{k}_W\xi }{\mu}_{max}\frac{N}{k_N\,+\,N}P-{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_PP\\ \qquad +\,{\partial}_{\xi}\left(D{\partial}_{\xi }P\right)\\{}{\partial}_tZ=\alpha{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_ZZ+{\partial}_{\xi}\left(D{\partial}_{\xi }Z\right)\\{}{\partial}_t{N}=-\rho{e}^{-{k}_W\xi }{\mu}_{\mathrm{m}\mathrm{ax}}\frac{{N}}{{{k}}_{{N}}+{N}}{P}\\ \qquad +\,\left(1-\alpha \right){{h}}_{\mathrm{m}\mathrm{ax}}\left(1-{{e}}^{-{aP}}\right){Z}+{{m}}_{{P}}{P}\\ \qquad +\,{{m}}_{{Z}}{Z}+{\partial}_{\xi}\left(D{\partial}_{\xi }N\right)\\{}\rho (I)=\frac{\Gamma I}{\sqrt{\mu_{max}^2\,+\,{\Gamma}^2{I}^2}}\\{}D\left(\xi, t\right)={D}_{min}-\frac{D_{min}-{D}_{max}(t)}{2}\left(1+\tanh \left(\frac{M(t)-\xi }{\psi}\right)\right)\\{}{D}_{max}(t)={10}^{-2}M{(t)}^2\end{array}}\!\!\!\!\!\!\!\!\!\!\!\! \end{eqnarray}$$

(6)

The phytoplankton response to light is expressed as |$\rho{e}^{-{k}_W\xi }$|^, which accounts for light attenuating exponentially with depth. In (6), the expression for the parameter |$\rho$| is used as in Kuhn et al. (2015) and Smith (1936). The diffusion coefficient |$D(\xi, t)$|⁠, a simplified version of the diffusion coefficient from Gibson et al. (2005), accounts for low diffusion in the lower ocean and higher diffusion in the upper mixed layer. The upper mixed layer has diffusion |$\mathcal{O}({D}_{max})$| and the lower ocean has diffusion |$\mathcal{O}({D}_{min})$|⁠. |$M(t)$| is the mixed layer depth (i.e. the depth of the bottom of the upper mixed layer). The approximation for |${D}_{max}$| is like that of Kuhn et al. (2015) but we include the coefficient 10⁻² to keep |${D}_{max}$| in the range of diffusion levels suggested in Edwards et al. (2000). When |${D}_{max}$| is of order 10⁻⁵ m²s⁻¹ or smaller and all other parameters are fixed, diffusion is weak and model predictions behave similarly to those governed solely by the reaction terms, that is the ordinary differential equation system (3). When |${D}_{max}$| is significantly larger than 10⁻³ m²s⁻¹ diffusion overwhelms the system. The parameter |$\psi$| controls the thickness of the transition layer between the mixed layer and the lower ocean. We choose |$\psi =25$|m for a smooth transition from the mixed layer to the lower ocean. This is computationally suitable for our vertical grid size of |$\Delta \xi =5$|m. Smaller values of |$\psi$| necessitates a smaller grid size and consequently, longer running time for the algorithm. Figure 3 shows a profile of the depth-dependent solar radiation intensity and of the diffusion coefficient. The diffusion coefficient is designed to separate the mixed layer from the lower ocean. Table II gives the descriptions of and the values chosen for the additional parameters we have introduced since expanding Equation (3) into the reaction–diffusion system (6). No-flux boundary conditions, |${\partial}_{\xi }P={\partial}_{\xi }Z={\partial}_{\xi }N=0$|⁠, are imposed at the surface, |$\xi =0$|⁠, and at the bottom, |$\xi =H=350$| m. The system is solved using a Crank–Nicolson scheme (LeVeque, 2007) with an evenly partitioned vertical grid size of Δξ = 5 m.

Fig. 2

Phase plane illustration of the solutions of Equation (3). The phytoplankton-only equilibrium (green cross) is a saddle point, and the coexisting equilibrium (red “X”) is a stable fixed point. Trajectories that start close to the stable manifold (horizontal axis, i.e. phytoplankton-only populations) will move toward the phytoplankton-only equilibrium but then turn away from that toward the coexisting equilibrium.

Fig. 3

Depth profile of solar radiation intensity (solid black) and diffusion coefficient (dashed blue) with surface solar radiation intensity I = 50 Em⁻²d⁻¹, |${k}_W=0.06$| m⁻¹, |${D}_{max}={10}^{-3}$| m²s⁻¹ and mixed layer depth M = 40 m for purposes of illustration. Light attenuates exponentially with depth. The mixed layer depth M is at the inflection point of the diffusion profile and at the center of the transition layer. M varies seasonally and geographically and typically ranges from 20 to 300 meters. Redrawn after Cowall et al. (2019).

Table II

Parameter descriptions and values introduced in Equations (6)

Parameter	Description	Value	Units
|${k}_W$|	Light attenuation constanta	0.02–0.06	m−1
|$\psi$|	Transition layer parameterb	25	m
|${D}_{min}$|	Diffusion level in the lower oceanc	10−6	m2 s−1
|$H$|	Water column depthd	350	m

a

Range as indicated from https://oceancolor.gsfc.nasa.gov/.

b

Value chosen for a smooth transition from mixed layer to lower ocean.

c

Value chosen from the lower range of diffusion values from Edwards et al. (2000).

d

Value chosen to accommodate most values of mixed layer depth M.

Bloom mechanics in a reaction–diffusion system

The reaction–diffusion system has bloom mechanics similar to that of the system of ordinary differential Equations (3). To show this, first, the time dependence of the right-hand side of (6) is removed and the system is reduced to two (partial differential) equations; simplifying the system saves computational time when computing the steady states, giving the system (7). We will observe the effect of the steady states on the time-dependent solutions of (7).

$$\begin{equation} {\displaystyle \begin{array}{l}{\partial}_tP=\rho{e}^{-{k}_W\xi }{\mu}_{max}\frac{N_T-P-Z}{k_N+{N}_T-P-Z}P-{h}_{max}\left(1-{e}^{- aP}\right)Z\\ \qquad -\,{m}_PP+{\partial}_{\xi}\left(D{\partial}_{\xi }P\right),\\{}{\partial}_tZ=\alpha{h}_{max}\left(1-{e}^{- aP}\right)Z-{m}_ZZ+{\partial}_{\xi}\left(D{\partial}_{\xi }Z\right),\\{}D\left(\xi \right)={D}_{min}-\frac{D_{min}-{D}_{max}}{2}\left(1+\tanh \left(\frac{M-\xi }{\psi}\right)\right).\end{array}} \end{equation}$$

(7)

System (7) has three depth-dependent equilibrium states whose stability properties depend on the parameters in Table I and Table II: a plankton-free equilibrium, a phytoplankton-only equilibrium, and a coexisting phytoplankton and zooplankton equilibrium (Cowall et al., 2019). To demonstrate the bloom mechanics in the reaction–diffusion system, a coexisting equilibrium state of (7) subject to low (winter-level) sunlight is used as initial conditions. Such equilibrium has a low concentration of zooplankton. Then the autonomous system (7) is solved subject to higher (summer-level) sunlight. The solution is expected to include a phytoplankton bloom like in Fig. 1.

The techniques of Edwards et al. (2000) are utilized to compute the equilibrium to be used for initial conditions. For illustrative purposes, |${h}_{max}$| is set to 0.75 d⁻¹ and |${k}_W$| is set to 0.06 m⁻¹. Mixed layer depth |$M$| is set to 130 m, which is a common winter value in the North Atlantic, and |${D}_{max}={10}^{-3}$| m² s⁻¹, which is in the upper range of diffusion values in Edwards et al. (2000). We write system (7) as

$$\begin{equation} {\displaystyle \begin{array}{c}{\partial}_tP=F\left(P,Z\right)+{\partial}_{\xi}\left(D{\partial}_{\xi }P\right)\\{}{\partial}_tZ=G\left(P,Z\right)+{\partial}_{\xi}\left(D{\partial}_{\xi }Z\right).\end{array}} \end{equation}$$

(8)

After discretizing depth, we create an initial guess by calculating the equilibrium |$\Big({\overset{\sim }{P}}_j,{\overset{\sim }{Z}}_j\Big)$| at each depth |${\xi}_j,j=1,\dots, m$|⁠, with |$D\equiv 0$| (i.e. we analytically calculate the equilibrium for Equation (9)).

$$\begin{equation} {\displaystyle \begin{array}{l}\frac{d{P}_j}{dt}=\rho{e}^{-{k}_W{\xi}_j}{\mu}_{max}\frac{N_T-{P}_j-{Z}_j}{k_N+{N}_T-{P}_j-{Z}_j}{P}_j\\ \qquad-{h}_{max}\left(1-{e}^{-a{P}_j}\right){Z}_j -\,{m}_P{P}_j\\{}\frac{d{Z}_j}{dt}=\alpha{h}_{max}\left(1-{e}^{-a{P}_j}\right){Z}_j-{m}_Z{Z}_j\end{array}} \end{equation}$$

(9)

We define |$\overrightarrow{P}$| and |$\overrightarrow{Z}$| as the vectors with entries |${\overset{\sim }{P}}_j$| and |${\overset{\sim }{Z}}_j$| (⁠|$j=1,\dots, m$|⁠), repspecively. Using (⁠|$\overrightarrow{P},\overrightarrow{Z}$|⁠) as an initial guess, we use a nonlinear equation solver to solve the system (10)

$$\begin{equation} {\displaystyle \begin{array}{c}F\left(P,Z\right)+{\partial}_{\xi}\left(D{\partial}_{\xi }P\right)=0\\{}G\left(P,Z\right)+{\partial}_{\xi}\left(D{\partial}_{\xi }Z\right)=0.\end{array}} \end{equation}$$

(10)

Figure 4a shows depth profiles of the computed coexisting equilibrium of system (7) with |$\rho =0.35$|⁠, |${h}_{max}=0.75$| d⁻¹, |$M=130$| m and |${D}_{max}={10}^{-3}$| m²s⁻¹. The zooplankton concentration is low at all depths in the water column, a prime condition for a bloom to begin. This equilibrium is used as the initial condition for system (7) with the same parameterization except that |$\rho$| is chosen to have the fixed value of 1. Similarly, Fig. 4b shows the depth profile of the coexisting equilibrium with the same parameterization as Fig. 4a except that |$\rho$| is chosen to have the fixed value of 1. In this case, the solution of (7) is expected to converge to the equilibrium in Fig. 4b.

Fig. 4

Depth profiles of the coexisting equilibrium for system (7) with |${D}_{max}={10}^{-3}$| m²s⁻¹, M = 130 m, |${k}_W=0.06$| m⁻¹ and |${h}_{max}=0.75$| d⁻¹. In (a), ρ = 0.35 and in (b), ρ = 1.

One would expect the plankton concentrations in the mixed layer to be homogeneous but this is not the case in the depth profiles of Fig. 4. This is in large part due to the light attenuation parameter |${k}_W$| having a greater influence on the vertical structure of plankton concentration than does diffusion or mixed layer depth. This is evident in Fig. 5, where |$\rho$| is set to 1 and the vertical axis spans the entire depth of the water column down to 350 m. In Fig. 5a, |${k}_W$| is reduced from 0.06 to 0.02 m⁻¹, which as expected, allows for existence of plankton deeper in the water column. In Fig. 5b, |${k}_W$| is reduced to 0.01 m⁻¹, which is below the range indicated from https://oceancolor.gsfc.nasa.gov/ for the North Atlantic region studied in this paper. Here, the phytoplankton concentration is essentially homogeneous from the surface to a depth of ~250 m, well below the mixed layer depth of 130 m. Although it is not expected for light attenuation to be this low in the North Atlantic, this outcome highlights the model’s sensitivity to the light attenuation coefficient |${k}_w$|⁠.Figure 6 shows the predicted surface concentrations as a function of time as predicted by the autonomous (no explicit time dependence) system (7) with |$\rho \equiv 1$| and initial conditions as shown in Fig. 4a. As expected, the phytoplankton bloom before the system subsides to the coexisting equilibrium state in Fig. 4b.

Fig. 5

Depth profiles for the coexisting equilibrium for system (7) with |${D}_{max}={10}^{-3}$| m²s⁻¹, M = 130 m, ρ = 1 and |${h}_{max}=0.75$| d⁻¹. In (a), |${k}_W=0.02$| m⁻¹ and in (b), |${k}_W=0.01$| m⁻¹.

Fig. 6

Phytoplankton and zooplankton surface concentrations computed from the depth-dependent Equation (7) with |${D}_{max}={10}^{-3}$| m²s⁻¹, M = 130 m, |${h}_{max}=0.75$| d⁻¹, |${k}_W=0.06$| m⁻¹ and ρ = 1, and initial conditions as shown in Fig. 4a and conditions on Day 90 shown in Fig. 4b.

Figure 7a–d shows the phytoplankton depth profiles from days 1, 10, 30 and 50 of the solution shown in Fig. 6 along with the phytoplankton depth profiles from the phytoplankton-only equilibrium and coexisting equilibrium with |$M=130$| m, and |${D}_{max}={10}^{-3}$| m²s⁻¹, |${h}_{max}=0.75$| d⁻¹ and |$\rho =1$|⁠. Figure 7e is a model temporal cross-section that shows phytoplankton concentration from this solution as a function of depth and time. This cross-section shows the temporal evolution of the depth-dependent phytoplankton concentration but the depth profiles in Fig. 7a–d are necessary to show how the solution compares to the steady states. Note that |$M$| and |${D}_{max}$| are kept constant. This is done so that the reader does not have to keep track of more than one moving curve in Fig. 7 or in the supplemental movie. One can see how the bloom works in the reaction–diffusion model when only |$\rho$| varies. However, the bloom will be affected by varying |$M$| and |${D}_{max}$| as will be seen in the next section when data will be used to parameterize these two quantities. On Day 1 (Fig. 7a), the phytoplankton concentration is relatively low at all depths and close to the phytoplankton depth profile of the coexisting equilibrium. On Day 10 (Fig. 7b), the phytoplankton concentration at all depths has risen relatively close to the levels of the phytoplankton depth profile of the phytoplankton-only equilibrium. With a dearth of grazing pressure, increasing sunlight allows the phytoplankton biomass to approach its maximum level that the water column can support. Figures 6 and 7e show that on Day 10, the surface phytoplankton is near its peak concentration. On Day 30 (Fig. 7c), the phytoplankton concentration at all depths has fallen to levels close to the corresponding levels of the phytoplankton depth profile of the coexisting equilibrium. As grazing pressure increases from rebounding zooplankton, the phytoplankton biomass declines toward a level such that the ecosystem is balanced. In Figures 6 and 7e, the surface phytoplankton concentration on Day 30 has fallen just below what appears to be its level in a coexisting equilibrium state; the system appears to essentially achieve this coexisting equilibrium by Day 90. On Day 50 (Fig. 7d), the phytoplankton concentration at all depths almost matches the corresponding levels of the phytoplankton depth profile of the coexisting equilibrium. Accordingly, the surface phytoplankton concentration in Fig. 6 and Fig. 7e also appears to be very close to its level in the coexisting equilibrium. Evidently, when the coexisting equilibrium of the NPZ reaction–diffusion system (7) is stable, the phytoplankton-only equilibrium behaves like its counterpart in system (1). Phytoplankton blooms in system (7) materialize in a manner analogous to the way they do in the system of ordinary differential Equation (3). A key difference is that the attracting states of the reaction–diffusion NPZ model are influenced not only by seasonally changing solar radiation, but also by diffusion, and by mixed layer depth (Cowall et al., 2019).

Fig. 7

Depth profiles of the phytoplankton component of the phytoplankton-only equilibrium and coexisting equilibrium of the system (7) with M = 130 m, |${D}_{max}={10}^{-3}$| m²s⁻¹, |${h}_{max}=0.75$| d⁻¹ and ρ = 1 along with depth profiles of the solution of (7): Day 1 (a), Day 10 (b), Day 30 (c) and Day 50 (d). Temporal cross-section (e) of the solution.

Fig. 8

Time series of solution of system (6) with I(t) described by PAR data and M(t) described by mixed layer depth data. All data averaged over the region 40–45°N, 25–35°W in the North Atlantic. Parameter values from Table I and Table II except for |${k}_W=0.02$| m⁻¹ and |${h}_{max}=0.5$| d⁻¹. Top: time series of PAR data with sinusoidal best-fit curve. Second from top: time series of mixed layer depth data with sinusoidal best-fit curve. Second from bottom: time series of surface model phytoplankton concentration predictions compared with chlorophyll data. Bottom: time series of surface model phytoplankton and zooplankton predictions.

The supplemental movie provides a depth profile animation of the phytoplankton predictions of Equation (7) with initial conditions |$(P(\xi ),Z(\xi ))$| from Fig. 4a and |$\rho =1$|⁠. The low initial zooplankton concentration and high sunlight intensity cause the phytoplankton depth profile to move toward the phytoplankton-only equilibrium depth profile. Then the phytoplankton depth profile turns toward the coexisting equilibrium depth profile and converges to it.

RESULTS AND DISCUSSION

In system (6), solar radiation data are incorporated to describe the quantity |$I(t)$| and mixed layer depth data are used to describe the quantity |$M(t)$|⁠. The satellite data PAR and chl-ocx (chlorophyll) collected by the NASA MODIS-Aqua satellite is used. The PAR and chlorophyll data are available at https://oceancolor.gsfc.nasa.gov/. The mixed layer depth data are retrieved from HYbrid Coordinate Ocean Model at https://hycom.org/. The region 25–35°W, 40–45°N in the North Atlantic Ocean is investigated. Time series of the horizontally averaged PAR and chlorophyll data from 2009 through 2015 are used for this region. The predictions of model surface concentrations of phytoplankton from Equation (6), using PAR data for |$I(t)$| and mixed layer depth data for |$M(t)$|⁠, are compared with a time series of the chlorophyll data. Figure 8 shows the predictions of the system (6) on the North Atlantic region, using |${h}_{max}=0.5$| d⁻¹ and |${k}_W=0.02$| m⁻¹. Sinusoidal best-fit approximations for PAR and mixed layer depth, shown in Fig. 8, are used to describe the quantities |$I(t)$| and |$M(t)$|⁠, respectively, to avoid noise in the solution. In the second plot from the bottom, the annual pattern of model phytoplankton concentration more closely resembles the chlorophyll data than it does in Fig. 1. The rapid increase in phytoplankton concentration begins after the winter solstice, instead of before as shown in Fig. 1. However, the peaks of the annual spring blooms are still earlier for the model than for the data. This is likely due to our simplified model’s lack of inclusion of effects from nutrient influx, thermal convection, ocean temperature, and sinking of phytoplankton among others.

As mentioned earlier, during the winter deepening of the mixed layer depth, phytoplankton concentrations may be decreasing, whereas the depth-integrated stocks are increasing. As shown in Fig. 9, this feature is captured in model (6) with parameter values set as shown in Fig. 8. In Fig. 9, there are winter periods when the model phytoplankton concentration averaged over the mixed layer (middle plot) decreases, whereas the integrated mixed layer phytoplankton biomass (bottom plot) increases. It is also seen in Fig. 9, unlike Fig. 1, that phytoplankton concentration begins to increase before the winter solstice and while the mixed layer depth is deepening. These features are consistent with the DRH.

Fig. 9

Comparison of time series of average model phytoplankton concentration in mixed layer (middle) and depth-integrated model phytoplankton biomass in mixed layer (bottom). Parameter values from Table I and Table II except for |${k}_W=0.02$| m⁻¹ and |${h}_{max}=0.5$| d⁻¹, same as shown in Fig. 8. Top plot is mixed layer depth. All data averaged over the region 40–45°N, 25–35°W in the North Atlantic.

Fig. 10

Stability diagram for equation system (7) with M = 130 m, |${k}_W=0.02$| m⁻¹, a = 1 (mmol N m⁻³)⁻¹ and |${h}_{max}=0.5$| d⁻¹. The white region (upper left) is where the plankton-free equilibrium (no phytoplankton or zooplankton survives) is stable. The aqua region is where the phytoplankton-only equilibrium is stable. The yellow region (right) is where the coexisting equilibrium is stable. The dashed curve is the seasonal trajectory of sinusoidal approximations for light response ρ and upper mixed layer diffusion |${D}_{max}$| based on the corresponding data from the North Atlantic region 40–45°N, 25–35°W. Four distinct symbols indicate points on the trajectory representing the vernal and autumnal equinoxes and the winter and summer solstices.

To further study the effect of PAR, diffusion and mixed layer depth data on the dynamics of the model Equation (7), we construct stability plots as is done in Cowall et al. (2019). Figure 10 is a stability diagram constructed for system (7) with the same parameter values used to produce the results in Fig. 8. The parameter |$\rho$| is plotted against |${D}_{max}$| to show the regions of stability for the plankton-free equilibrium, the phytoplankton-only equilibrium and the coexisting equilibrium. The mixed layer depth |$M$| is fixed at 130 m although it is known to change with the seasons. Higher values for |$M$| have an effect similar to that of higher values of |${D}_{max}$| as demonstrated in Cowall et al. (2019). The white region (upper left) is where the plankton-free equilibrium is stable. The aqua region is where the phytoplankton-only equilibrium is stable. The yellow region (right) is where the coexisting equilibrium is stable. No regions were found where the coexisting equilibrium exists and is unstable. The dashed curve is the seasonal trajectory of sinusoidal approximations for light response |$\rho$| and upper mixed-layer diffusion |${D}_{max}$|⁠. This diagram and Fig. 12 are constructed using the techniques from Edwards et al. (2000) and Cowall et al. (2019). Figure 10 shows that the |${D}_{max}$| verses |$\rho$| trajectory lands close to the border between the region of stability for the phytoplankton-only equilibrium and the region of stability for the coexisting equilibrium on December 20 but turns directly away from the border immediately afterward. Under these circumstances, the phytoplankton biomass is expected to increase almost immediately after the winter solstice in our model. This is what happens in the second plot from the bottom of Fig. 8 as the phytoplankton only equilibrium is first an attractor, that is the phytoplankton population grows with little change in zooplankton, but then the zooplankton population begins to grow with the result that the solution trajectory pulls away from the phytoplankton-only equilibrium. As the winter solstice approaches, sunlight is decreasing but mixed layer diffusion is increasing. After the winter solstice, the increase in phytoplankton accelerates as the sunlight increases.

Fig. 11

Time series of solution of system (6) with I(t) described by PAR data and M(t) described by mixed layer depth data. All data averaged over the region 40–45°N, 25–35°W in the North Atlantic. Parameter values from Table I and Table II except for |${k}_W=0.04$| m⁻¹, a = 0.5 (mmol N m⁻³)⁻¹ and |${h}_{max}=0.75$| d⁻¹. Top: time series of PAR data with sinusoidal best-fit curve. Second from top: time series of mixed layer depth data with sinusoidal best-fit curve. Second from bottom: time series of surface model phytoplankton concentration predictions compared with chlorophyll data. Bottom: time series of surface model phytoplankton and zooplankton predictions.

The model phytoplankton biomass increases occur later in the year if the |$\log{D}_{max}$| verses |$\rho$| trajectory turns away from the border between the domain of stability of the phytoplankton-only equilibrium and the domain of stability of the coexisting equilibrium later than in Fig. 10. This border can be moved to the right by increasing |${k}_W$| or decreasing |$a$| as shown in Cowall et al. (2019). Figure 11 shows the predictions of the system (6) in the North Atlantic region using |${h}_{max}=0.75$| d⁻¹, |$a=0.5$| (mmol N m⁻³)⁻¹ and |${k}_W=0.04$| m⁻¹. The value of |${h}_{max}$| is chosen because, if that change were not made, then the new graph with the indicated changes in |$a$| and |${k}_W$| would result in blooms lasting too long (not shown). The second plot from the bottom of Fig. 11 illustrates that the timing of the phytoplankton decline has improved over that of Fig. 8. However, two problems materialize with this change of parameterization. First, there is a late-year decline in the model phytoplankton concentration that is not indicated in the chlorophyll data. Second, the bottom plot of Fig. 11 shows the zooplankton concentration declining to |$4\times{10}^{-4}$| mmol N m⁻³, an unrealistically low level. Consequently, with negligible grazing pressure and a higher maximum grazing rate, the phytoplankton concentration increases more rapidly, remains at its peak longer and declines faster. Hence the phytoplankton peaks in Fig. 11 are broader than those in Fig. 8. These two problems can be explained by the stability diagram in Fig. 12.

Fig. 12

Stability diagram for equation system (7) with M = 130 m, |${k}_W=0.04$| m⁻¹, a = 0.5 (mmol N m⁻³)⁻¹ and |${h}_{max}=0.75$| d⁻¹. The white region (upper left) is where the plankton-free equilibrium is stable. The aqua region (between the white and the yellow) is where the phytoplankton-only equilibrium is stable. The yellow region (right) is where the coexisting equilibrium is stable. The dashed curve is the seasonal trajectory of sinusoidal approximations for light response ρ and upper mixed layer diffusion |${D}_{max}$| based on the corresponding data from the North Atlantic region 40–45°N, 25–35°W. Four distinct symbols indicate points on the trajectory representing the vernal and autumnal equinoxes and the winter and summer solstices.

The stability diagram of system (7) with parameters as shown in Fig. 11 along with a trajectory of |${D}_{max}$| versus |$\rho$| based on the sinusoidal approximation of the data from the North Atlantic region is shown in Fig. 12. The phytoplankton bloom of the model occurs well after the winter solstice. This can be understood because the |${D}_{max}$| versus |$\rho$| trajectory moves away from the domain of stability for the phytoplankton-only equilibrium and into the domain of stability for the coexisting equilibrium about half-way between the winter solstice and the vernal equinox. The late year decline in phytoplankton concentration can be explained by the |${D}_{max}$| verses |$\rho$| trajectory entering the domain of stability of the phytoplankton-only equilibrium and getting close to the domain of stability of the plankton-free equilibrium at the winter solstice. The trajectory entering the domain of stability for the phytoplankton-only equilibrium causes the zooplankton concentration to rapidly decline to unrealistically low levels. The zooplankton loss rate exceeds its growth rate to the point that the zooplankton population begins to collapse. Therefore, accurately simulating blooms requires that the |${D}_{max}$| versus |$\rho$| trajectory must either spend a minimal amount of time in the domain of stability for the phytoplankton-only equilibrium or stay out of it completely. Ecologically this means zooplankton can only sustain a scarcity of phytoplankton (due to low sunlight) or excessive diffusion for a limited time.

Although the bloom mechanics described by the system of ordinary differential Equation (6) shed light on the dynamical causes of phytoplankton blooms in the North Atlantic according to the DRH, the timing of blooms simulated by the system (6) is early (Figs 8 and 11). A more accurate estimation of the seasonal changes in diffusion in the upper mixed layer of the ocean or of |${k}_W$| could be beneficial. Incorporating ocean temperature data into the phytoplankton growth term, as is done in Kuhn et al. (2015), could also affect the results. Non-dimensionalizing system (6) could be useful for simplifying the model and better understanding of the impact of parameter groups and parameter variations. The model described by the reaction–diffusion system (6) also neglects the sinking of phytoplankton, entrainment of nutrients from the deep ocean, and light attenuation dependence on phytoplankton concentration. Such processes are known to happen in the ocean and have been simulated by plankton ecology models.

CONCLUSIONS

Many theories about the underlying causes of seasonal, marine phytoplankton blooms have emerged in the last century. The CDH, CTH and disturbance-recoupling hypothesis (DRH) argue that physical processes, such as shoaling of mixed-layer depth and decreasing turbulent mixing (diffusion), cause blooms by boosting phytoplankton production. Behrenfeld and Boss (2014) argue that a disruption of the predator–prey balance is needed in the explanation of a phytoplankton bloom. In this paper, we follow up on this notion by identifying bloom mechanics in an ordinary differential equation system and in a reaction diffusion system that can cause phytoplankton blooms. Whatever the underlying causes of blooms may be, whether they can be explained by an existing theory or a combination of existing theories, we propose that phytoplankton blooms correspond to the transient nature of the attraction of a phytoplankton-only equilibrium state of the planktonic ecosystem. Many ordinary differential equation plankton ecology models that include grazing by zooplankton have a phytoplankton-only equilibrium. This steady state is often computed but not included in the discussion of how blooms arise. Ideally, a complete narrative of a bloom in such a model, with or without periodic forcing would include the role of the phytoplankton-only equilibrium in the dynamics of the bloom. This paper not only fills this gap but also explains the analogous role of this steady state in a depth-dependent NPZ model. We evaluated the predictive power of the dynamical properties in our model by analyzing the effect of seasonal data of solar radiation and mixed layer depth on the dynamical causes of phytoplankton blooms. Our model is a conceptual model which has limitations, many of which were discussed in Cowall et al. (2019). Additional limitations should be taken into consideration when interpreting the blooming process described in this paper. Because of our model’s simplicity, it fails to imitate phytoplankton blooms to the degree of accuracy of more complex models, such as that of Kuhn et al. (2015). Behrenfeld et al. (2017) demonstrate that phytoplankton division rates and loss rates are well-coupled. However, the bottom graphs of Figs 8 and 11 show a significant decoupling between phytoplankton biomass and zooplankton biomass. This is a consequence, in part, of the system (6) having a preference for a phytoplankton-only equilibrium under low sunlight and high diffusion, which diminishes the zooplankton concentration, which is congruent with the major tenets of DRH. Also, the zooplankton have a slower maximum growth rate and a higher mortality rate than do the phytoplankton, which causes a delay in the zooplankton’s recovery. If marine phytoplankton bloom cycles involve a significant reduction of zooplankton then the zooplankton would have to recover faster than they do in this model. The simplifying assumption that |$N+P+Z={N}_T$| implies that the system is closed; no new nutrients enter the system. Although this assumption allows for a thorough analysis of the model’s dynamical properties, the assumption is only ecologically reasonable for short periods of time. It is not a realistic assumption for the entire annual bloom cycle because nutrients enter the system through entrainment and leave the system by sinking throughout the year.

Despite its limitations, our model provides informative results because it allows for the analysis of its dynamics. These dynamical properties also appear to be present in more complex models. For example, Kuhn et al. (2015) used a parameter optimization scheme on an advection–diffusion–reaction nutrient-phytoplankton-zooplankton-detritus (NPZD) model to produce results that mimic phytoplankton blooms in the North Atlantic quite well. Their results, however, suggest that the dynamical processes described in this paper are involved. Before a bloom, as in this paper, their phytoplankton levels are low, and their zooplankton levels are low but reasonable. This suggests that low sunlight and high diffusive turbulence have pushed the system close to a point where a phytoplankton-only equilibrium is the attracting state. The following increase in phytoplankton biomass is likely due to the temporary influence of the existence of a saddle point phytoplankton-only equilibrium, which temporarily attracts then repels the system.

If it is true that phytoplankton blooms in the North Atlantic are primarily caused by a disruption in the predator–prey balance as suggested by Behrenfeld and Boss (2014), this paper provides a detailed explanation of the underlying mechanism for phytoplankton bloom formation with reference to attracting and repelling states of the plankton system. We postulate that the disruption in the predator–prey balance suggested by Behrenfeld and Boss (2014) is set up by winter conditions rendering low a concentration of phytoplankton and an even lower (but reasonable) concentration of zooplankton. The return of spring conditions (increased sunlight and decreased diffusion) causes the phytoplankton only equilibrium to temporarily attract and then repel the system toward a co-existing equilibrium in late spring. In autumn, the system is disrupted again by decreased light and increased diffusion, allowing for the possibility of the phytoplankton only equilibrium to attract the system in the spring, thus completing the cycle. Although our simple model (6) does not capture the secondary details of the bloom process, it can be useful for explaining the mechanics of more complex models and for understanding marine planktonic ecosystems.

FUNDING

This work was supported by the University of Delaware Department of Mathematical Sciences.