# Exploration of the critical depth hypothesis with a simple NPZ model

Marina Lévy

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marina Lévy. Exploration of the critical depth hypothesis with a simple NPZ model. ICES Journal of Marine Science, 2015, 72 (6), pp.1916-1925. 10.1093/icesjms/fsv016. hal-01233869

HAL Id: hal-01233869

https://hal.science/hal-01233869v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ICES Journal of Marine Science (2015), 72(6), 1916–1925. doi:10.1093/icesjms/fsv016

**Contribution to the Themed Section: 'Revisiting Sverdrup's Critical Depth Hypothesis'**
**Food for Thought**

**Exploration of the critical depth hypothesis
with a simple NPZ model**

Marina Lévy*

Sorbonne Université (UPMC, Paris 6)/CNRS/UPMC/IRD/MNHN, Laboratoire d'Océanographie et du Climat (LOCEAN), Institut Pierre, Simon Laplace (IPSL), Paris Cedex 05 75252, France

*Corresponding author: e-mail: marina.levy@upmc.fr

Lévy, M. Exploration of the critical depth hypothesis with a simple NPZ model. – ICES Journal of Marine Science, 72: 1916–1925.

Received 15 September 2014; revised 7 January 2015; accepted 12 January 2015; advance access publication 16 February 2015.

The critical depth hypothesis (CDH) is a predictive criteria for the onset of phytoplankton blooms that comes from the steady-state analytical solution of a simple mathematical model for phytoplankton growth presented by Sverdrup in 1953. Sverdrup's phytoplankton-only model is very elementary compared with state-of-the-art ecosystem models whose numerical solution in a time-varying environment do not systematically conform to the CDH. To highlight which model ingredients make the bloom onset deviate from the CDH, the complexity of Sverdrup's model is incrementally increased, and the impact that each new level of complexity introduced is analysed. Complexity is added both to the ecosystem model and to the parameterization of physical forcing. In the most complete experiment, the model is a one-dimensional Nutrient-Phytoplankton-Zooplankton model that includes seasonally varying mixed layer depth and surface irradiance, light and nutrient limitation, variable grazing, self-shading, export, and remineralization. When complexity is added to the ecosystem model, it is found that the model solution only marginally deviates from the CDH. But when the physical forcing is also changed, the model solution can conform to two competing theories for the onset of phytoplankton blooms—the critical turbulence hypothesis and the disturbance recovery hypothesis. The key roles of three physical ingredients on the bloom onset are highlighted: the intensity of vertical mixing at the end of winter, the seasonal evolution of the mixed-layer depth from the previous summer, and the seasonal evolution of surface irradiance.

**Keywords:** bloom onset, critical depth, critical turbulence, disturbance recovery, NPZ model.

**Introduction**

More than 60 years ago, Sverdrup (1953) presented a simple, water column-based, mathematical model that tied together the roles of vertical mixing, light attenuation with depth and seasonally increasing light into a predictive criteria for the onset of the spring bloom. The critical depth hypothesis (CDH) is this model's solution. The CDH predicts that blooms start when the seasonal mixed-layer shoals above a critical depth ($D_C$), and offers a simple analytical expression for this depth.

In more recent years, observational evidence that biomass could start accumulating before restratification (Townsend et al., 1994; Boss and Behrenfeld, 2010) led researchers to revisit the CDH. Two other hypotheses were subsequently proposed to explain winter, rather than spring, bloom onsets. The critical turbulence hypothesis (CTH) attributes winter blooms to physical drivers (Huisman et al., 1999; Chiswell, 2011; Taylor and Ferrari, 2011; Brody and Lozier,

2014; Ferrari et al., 2014) and states that, rather than a critical depth, there should exist a critical level of turbulence (parameterized by a turbulent diffusivity $\kappa$) below which blooms occur. The emphasis of the CTH is on the distinction between apparently mixed layers (i.e. where density is homogeneous on the vertical but the turbulence is not necessarily strong) and regions with active vertical mixing (active mixing layers, characterized by large values of turbulent diffusivity $\kappa$) (Brainerd and Gregg, 1995; Franks, 2015). Note that the CTH has also been described as the “stratification-onset model” (Chiswell, 2011). On the other hand, the disturbance recovery hypothesis (DRH) highlights the importance of grazing control on bloom dynamics (Evans and Parslow, 1985; Banse, 1994; Marra and Barber, 2005; Behrenfeld, 2010; Boss and Behrenfeld, 2010; Behrenfeld and Boss, 2014). In essence, the DRH posits that winter entrainment of plankton-free water dilutes plankton populations, leading to reduced grazer encounters, and consequently net growth in the integrated population.
Though Sverdrup did not pose it this way originally, his model can be understood in terms of a differential equation for the time evolution of phytoplankton concentration (P)

$$ \partial_t P = (\mu - l)P + \partial_z(\kappa \partial_z P). \quad (1) $$

The parameters $\mu$ and $l$ are the phytoplankton growth and loss rates, respectively. The last term accounts for vertical mixing, and $\kappa = \kappa(z)$ is the vertical mixing coefficient. Note that in this equation, transport processes other than vertical mixing are omitted and $P$ is assumed to vary only with depth and time. Sverdrup solves this equation at steady-state, resulting in a necessary condition for the depth-integrated phytoplankton growth to exceed its losses. Importantly, Sverdrup makes a number of assumptions before solving this equation, among which are a constant loss rate $l$, and of a value of $\kappa$ large enough to ensure homogeneous distribution of phytoplankton within the surface mixed layer. This last assumption allows him to circumvent the details of vertical mixing by focusing on integral properties in an homogeneous surface layer. Hence, by construction, Sverdrup's model solution does not account for variable grazing, or for low values of $\kappa$, and thus cannot reproduce the CTH or DRH.

Contemporary ecosystem models basically follow the foundation of Equation (1), but allow for more flexibility in the description of $\mu$ and $l$ than Sverdrup's model. This flexibility is allowed by the consideration of additional state variables that feedback on $P$ through their contribution to $\mu$ and $l$, such as resources ($N$) and predators ($Z$). Models vary in complexity in terms of the number of variables, as well as in their expressions for $\mu$ and $l$. Importantly, they are also embedded in a physical model, with explicit representation of vertical mixing.

Numerical solutions of these models do not always conform to the CDH. Using a Sverdrup-like one-variable (i.e. phytoplankton only) ecosystem model embedded in a high-resolution three-dimensional upper-ocean turbulence model, Taylor and Ferrari (2011) showed departure from the CDH when turbulent mixing was below a critical turbulence level, which supports the CTH. On the other hand, in a one-dimensional 15-variable ecosystem model for the western Subarctic Pacific, Yoshie et al. (2003) found an increase in the integrated phytoplankton biomass as the mixed layer deepened due to a decrease of the grazing pressure by dilution. A similar result was obtained for the Subarctic Atlantic with the more complex 24-variable Biogeochemical Element Cycling Community Climate System Model (Behrenfeld et al., 2013). Along the same line, the 24-variable ecosystem model PISCES, in an idealized one-dimensional setup for the Southern Ocean (Llort et al., 2015) also supported the DRH. Finally, other models, such as the four-variable ecosystem model embedded in a coarse-resolution ocean general circulation model used by Dutkiewicz et al. (2001) to simulate the subpolar North-Atlantic bloom onset, supported the CDH. The above studies, in aggregate, highlight that each of these seemingly competing spring bloom drivers may be achieved depending on the form of the model. However, the models are too different to identify how specific biological and physical factors influence bloom onset.

In this note, I examine the bloom onset problem with a one-dimensional Nutrient-Phytoplankton-Zooplankton (NPZ) model (Franks, 2002), forced by vertical mixing in a time-varying mixed layer. The NPZ model, with only three state variables, is one of the simplest ecosystem models. Yet this set-up, which allows for variable grazing and variable $\kappa$, is sufficiently complex to explore the CDH, CTH, and DRH in one framework. Using a unified framework

allows me to identify the model ingredients that make the bloom onset deviate from the CDH. I start from the simplest model of Sverdrup, and progressively add complexity. This enables me to explore the impact that each new level of complexity introduces in a series of incremental numerical experiments. Complexity is added both to the ecosystem model (number of state variables, formulation of $\mu$ and $l$) and to the parameterization of physical forcing (from a mixed layer at steady-state to a mixed layer following a full seasonal cycle, and from constant incident irradiance to seasonal irradiance). This illustrative set of experiments is far from exhaustive. The intention is to shed light on the important roles of both the physical forcing and model structure, as well as to bridge the gap between Sverdrup's model solution and the different solutions obtained with more sophisticated ecosystem models.

After recalling the details of Sverdrup's model and presenting the NPZ model, Sverdrup's model assumptions are relaxed one by one, thus progressively adding complexity in both the model and physical forcing. The most complete experiment is obtained with a time-variable 1D-NPZ model, which, unlike Sverdrup's model, accounts for time-variable self-shading, explicit grazing, nutrient and light limitation, variable vertical mixing, and time-variable mixed-layer depth. Model solutions are presented at each step of the process; the bloom onset is diagnosed as the time at which water column growth and losses come in balance.

## Method

### Sverdrup's model

Sverdrup's model is based on an earlier model proposed by Riley (1946). It considers that the phytoplankton growth rate $\mu$ is proportional to light $I$ that the light-extinction coefficient $k$ is constant and that the phytoplankton loss rate $l$ is equal to a constant mortality rate $m_p$. Under these assumptions, $\mu$ and $l$ in Equation (1) are

$$ \mu = \alpha I, \quad (2) $$

$$ l = m_p \quad (3) $$

with

$$ I = I_0 e^{-kz}, \quad (4) $$

where $I_0$ is the fraction of surface irradiance usable for photosynthesis. Assuming that vertical mixing $\kappa$ is strong enough to distribute the organisms evenly within a surface mixed layer allows one to integrate Equation (1) from the surface to the bottom ($D > 0$) of the mixed layer, resulting in the expression

$$ \partial_t \bar{P} = \frac{\alpha I_0}{kD} (1 - e^{-kD}) \bar{P} - m_p \bar{P} \quad (5) $$

with

$$ \bar{P} = \int_{0}^{D} P \, dz \quad (6) $$

the total phytoplankton biomass within the mixed layer. Setting $\partial_t \bar{P} = 0$ results in an equation for a critical mixed-layer depth $D_C$ for which depth-integrated growth equals depth-integrated losses, i.e.

$$ \frac{\alpha I_0}{k D_C} (1 - e^{-k D_C}) = m_p. \quad (7) $$

The value of $D_C$ depends on the three model parameters $\alpha$, $m_p$, and $k$, and on the surface irradiance $I_0$. In Sverdrup's paper, the number
<table><caption>**Table 1.** Sverdrup's model parameters.</caption><thead><tr><th>Parameter name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Phytoplankton growth rate per unit light</td><td>α</td><td>3.66 × 10<sup>-7</sup></td><td>s<sup>-1</sup> W<sup>-1</sup> m<sup>2</sup></td></tr><tr><td>Phytoplankton mortality rate</td><td>m<sub>p</sub></td><td>5.80 × 10<sup>-7</sup></td><td>s<sup>-1</sup></td></tr><tr><td>Light attenuation</td><td>k</td><td>0.05</td><td>m<sup>-1</sup></td></tr><tr><td>Surface irradiance</td><td>l<sub>0</sub></td><td>50</td><td>W m<sup>-2</sup></td></tr><tr><td>Compensation irradiance</td><td>l<sub>c</sub></td><td>1.58</td><td>W m<sup>-2</sup></td></tr><tr><td>Critical depth</td><td>D<sub>c</sub></td><td>270</td><td>m</td></tr></tbody></table>

of parameters is reduced by assuming that the compensation irradiance $I_c = m_p/\alpha$ is known, where $I_c$ is defined as the irradiance at the depth for which $\mu(z) = l(z)$ (see Siegel et al., 2002, for the evaluation of the compensation irradiance from ocean colour observations). We assign these parameters such that both $k$ and $I_c$ are close to the values proposed by Sverdrup in his paper (Table 1). Note that Sverdrup uses different values for $k$ and we choose a value in the lower range, which leads to a critical depth of $D_c = 270$ m from Equation (7). Results with a larger value of $k$, leading to a shallower critical depth, will also be presented.

## NPZ model

Complexity is added to Sverdrup's model by introducing two extra state variables, nutrients (N) and zooplankton (Z). In the context of the bloom onset, introducing Z enables the model to account for time-dependent phytoplankton losses, which is necessary to allow the emergence of a model solution that satisfies the DRH. The need to introduce N is not as straightforward since the bloom always appears during the winter-to-spring transition when nutrients are plentiful. However, as we will see, a complete description of the bloom requires one to integrate the model over the full seasonal cycle. Hence, the seasonal cycling of nutrients—essentially the entrainment between fall and winter followed by export and remineralization between spring and summer—must also be accounted for.

The phytoplankton and zooplankton equations follow the general form of Equation (1), and are given by

$$ \partial_t P = (\mu_p - l_p)P + \partial_z(\kappa \partial_z P) \quad (8) $$

$$ \partial_t Z = (\mu_z - l_z)Z + \partial_z(\kappa \partial_z Z) \quad (9) $$

with plankton growth and loss terms given by

$$ \mu_p = \alpha K_i (1 - e^{-(I/K_i)}) \frac{N}{N + K_n} \quad (10) $$

$$ l_p = m_p + g_z Z \quad (11) $$

$$ \mu_z = g_z P \quad (12) $$

$$ l_z = \tau_z + m_z Z + (1 - a_z) g_z P. \quad (13) $$

We note that the phytoplankton growth rate $\mu_p$ incorporates two limitation terms that were absent in Sverdrup's model and which vary between 0 and 1: light saturated growth ($1 - e^{-(I/K_i)}$) and nutrient limitation ($N/N + K$). Moreover, the phytoplankton loss rate $l_p$ includes a term due to grazing $g_z Z$ which depends on the abundance of zooplankton. The zooplankton loss rate $l_z$ includes excretion ($\tau_z$), mortality ($m_z Z$) and non-assimilated grazing

<table><caption>**Table 2.** NPZ model parameters.</caption><thead><tr><th>Parameter name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Phytoplankton growth rate per unit light</td><td>α</td><td>3.66 × 10<sup>-7</sup></td><td>s<sup>-1</sup> W<sup>-1</sup> m<sup>2</sup></td></tr><tr><td>Nutrient limitation half-saturation</td><td>K<sub>n</sub></td><td>0.7</td><td>mmol m<sup>-3</sup></td></tr><tr><td>Light limitation half-saturation</td><td>K<sub>i</sub></td><td>33.</td><td>W m<sup>-2</sup></td></tr><tr><td>Phytoplankton mortality rate</td><td>m<sub>p</sub></td><td>5.80 × 10<sup>-7</sup></td><td>s<sup>-1</sup></td></tr><tr><td>Background light attenuation</td><td>k<sub>b</sub></td><td>0.05/0.15</td><td>m<sup>-1</sup></td></tr><tr><td>light attenuation by phytoplankton</td><td>k<sub>p</sub></td><td>0.03</td><td>m<sup>-1</sup> mmol<sup>-1</sup> m<sup>3</sup></td></tr><tr><td>Zooplankton maximal grazing rate</td><td>g<sub>z</sub></td><td>9.26 × 10<sup>-6</sup></td><td>s<sup>-1</sup></td></tr><tr><td>Assimilated food fraction by zooplankton</td><td>a<sub>z</sub></td><td>0.7</td><td></td></tr><tr><td>Zooplankton excretion rate</td><td>τ<sub>z</sub></td><td>5.80 × 10<sup>-7</sup></td><td>s<sup>-1</sup></td></tr><tr><td>Zooplankton mortality rate</td><td>m<sub>z</sub></td><td>2.31 × 10<sup>-6</sup></td><td>s<sup>-1</sup> mmol<sup>-1</sup> m<sup>3</sup></td></tr><tr><td>Export ratio</td><td>f<sub>e</sub></td><td>0.2</td><td></td></tr><tr><td>Euphotic layer depth</td><td>z<sub>b</sub></td><td>150</td><td>m</td></tr><tr><td>Remineralization exponential decay</td><td>r</td><td>0.858</td><td></td></tr></tbody></table>

($1 - a_z) g_z P$). Light is attenuated with depth as in Equation (4), but unlike in Sverdrup's model, the light attenuation coefficient depends on $P$,

$$ k = k_b + k_p P. \quad (14) $$

The nutrient pool is depleted by photosynthesis ($-\mu_p P$) and resupplied by the combined losses of phytoplankton and zooplankton after they are remineralized ($R$). There are no external sources of $N$. To implicitly account for the export of organic material out the euphotic layer by the sinking of particles, only a fraction ($1 - f_e$)$R$ of the remineralized material is locally resupplied to the nutrient pool; the rest of the remineralized material is vertically integrated ($f_e \int_0^H R dz$, with $H$ the depth of the water column) and redistributed below the euphotic depth $z_b$, following the Martin curve $f(z)$ (Martin et al., 1987). The resulting equation for $N$ is

$$ \partial_t N = -\mu_p P + (1 - f_e) R + \partial_z f \left( f_e \int_0^H R dz \right) + \partial_z (\kappa \partial_z N) \quad (15) $$

with

$$ R = l_z Z + l_p P \quad (16) $$

and

$$ f(z) = \min \left( 1, \left( \frac{z}{z_b} \right)^{-r} \right), \quad f(H) = 0. \quad (17) $$

The NPZ model presented above is a common tool in oceanographic research, and has many variants (Franks, 2002). The functional forms of the growth and loss terms (Equations (10)–(13)) and model parameters (Table 2) used here come from of the more complex Locean Ocean Biogeochemical System for Ecosystem and Resources model, which has been calibrated against biogeochemical data of the Northeast Atlantic spring bloom (Karleskind et al., 2011).
This NPZ model can be reduced to Sverdrup’s model by setting to zero the grazing rate g, the phytoplankton light attenuation k_p, the nutrient half-saturation constant K_n, and by replacing the light limitation term (1 - e^(-I/(K_n))) with its first order expansion at low light, I/K_n. The consequences of these simplifications will be investigated incrementally in the set of experiments with label ‘B’, described below.

Equations (8), (9), and (15) are solved with an implicit-in-time scheme, using a time-step of 1 h. The vertical resolution is 5 m, and there are 70 equally spaced vertical layers. Initial conditions for N, P, and Z are set to 4, 0.1, and 0.01 mmol m⁻³, respectively, homogeneously over the whole column; these values correspond to winter conditions.

## Bloom onset diagnostics

The bloom onset is assessed by examining the seasonal evolution of the total phytoplankton growth (μ, green lines in Figure 1) and loss (l, blue lines in Figure 1) rates. The onset is defined as the time at

which phytoplankton net growth (μ̄ - l̄) switches from negative to positive value. As noted by Chiswell (2013), integrating over the depth of the mixed layer (as in Sverdrup, 1953; Behrenfeld, 2010; Behrenfeld et al., 2013) does not account for the phytoplankton biomass which is detrained, diffused and slowly decays below the mixed layer when the mixed-layer shallows. Thus here, the phytoplankton growth and loss rates are integrated over the whole depth H of the water column rather than over the mixed-layer depth D [as in Equation (6)], e.g.

$$ \bar{\mu} = \frac{1}{\bar{P}} \int_{0}^{H} \mu_{p} P \, dz \qquad (18) $$

$$ \bar{l} = \frac{1}{\bar{P}} \int_{0}^{H} l_{p} P \, dz \qquad (19) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** For each series of experiments A (left hand column), B (middle column), and C (right hand column), the first line (a-c) shows the seasonal evolution of the mixed-layer depth D (in black) and critical depth D_c (in red). The other lines show the seasonal evolution of the total water-column phytoplankton growth rate (μ, in green) and loss rate (l, in blue) for each experiment in the corresponding series. The actual bloom onset occurs when net growth is positive, i.e. when the green curve comes above the blue one. The theoretical bloom onset predicted by the CDH occurs when the black curve comes above the red one (indicated by the dotted vertical lines in the B and C series). In (a), the thick black line shows D for experiment A2 and the thin black line shows D in A1 and A3; D_c is the same for A1, A2, and A3. In (b), D is the same for the four experiments B1-B4; the plain red line is D_c for B1 and B2, the dotted line for B3 and the dashed line for B4. In (c), D is the same for the three experiments C1-C3; the plain red line is D_c in C1, the dotted line in C2, and the dashed line in C3. In (e), the thick green line is the growth rate in B1, and the thin green line is the growth rate in B2; the loss rate is identical in B1 and B2.</figcaption></figure>
such that

$$ \frac{1}{\bar{P}} \partial_t \bar{P} = \bar{\mu} - \bar{l} \qquad (20) $$

with

$$ \bar{P} = \int_{0}^{H} P \, dz. \qquad (21) $$

With these definitions, the total stock of biomass in the water column, $\bar{P}$, increases as soon as $\bar{\mu}$ exceeds $\bar{l}$.

## Relaxing Sverdrup's hypothesis

The analytical expression for a critical depth relies on a number of assumptions, clearly identified in Sverdrup's 1953 paper, and used to construct his model. Paraphrasing Sverdrup, these assumptions are:

1. Phytoplankton develops in a thoroughly mixed layer at the surface of the ocean where turbulence is large enough to distribute the organisms evenly on the vertical.

2. Inorganic nutrients are fully available in the mixed layer and do not limit productivity.

3. Productivity is proportional to the light energy available for photosynthesis.

4. The light-extinction coefficient is constant in the mixed layer.

5. Grazers are not yet present in pre-bloom conditions.

It should be noted that these assumptions greatly facilitate the analytical derivation of the critical depth but are not necessary to simulate a bloom numerically. Here, using a set of numerical experiments (Table 3), the CDH is challenged by relaxing the above assumptions, one by one.

## Physical forcing

The physical drivers of the model are the vertical mixing coefficient $\kappa$ and surface irradiance $I_0$. In all experiments, $\kappa$ is set to $\kappa_{max}$ inside an upper mixed layer of depth $D$ and to a minimum value of $10^{-5} \text{ m}^2 \text{ s}^{-1}$ below.

Three series of experiments (labelled A, B, and C) are conducted in which three different time evolutions of the mixed layer are considered (Figure 1a-c, black curves). These experiments are characteristic of the situation in the North Atlantic. In the A series, the mixed-layer depth is kept constant in time. In the B series, the winter-to-spring transition of the mixed layer is included. Finally

in the C series, the seasonal cycle of the mixed layer from the previous summer is included. Surface irradiance is kept constant in time, except in the last two experiments of the C series (C2 and C3) where seasonal variations are also included (see Table 3).

## Results

### Series A: Constant mixed-layer depth

In the A series of experiments (Figure 1, left hand column), the mixed-layer depth is kept constant in time and Sverdrup's model is used. This series corresponds to the analytical solution of Sverdrup and serves to test that the numerical model solution is faithful to the analytical solution. The turbulent diffusivity $\kappa_{max}$ is set to $1 \text{ m}^2 \text{ s}^{-1}$, such that assumption #1 is met. This value is typical of a convective situation and corresponds to a time scale of the order of 0.1 d for a mixed-layer depth of 100 m. To achieve the test, two experiments are compared, one (A1) where $D$ is slightly deeper than $D_C$, the other (A2) where $D$ is slightly shallower. Results for A1 and A2 confirm that our numerical model reproduces the situation explored by Sverdrup: in A2, $\bar{\mu}$ is larger than $\bar{l}$ (Figure 1g): the condition for a bloom onset is satisfied; this is not the case in A1 (Figure 1d). Note that in A1 and A2, the net growth is constant and positive, so phytoplankton biomass would grow indefinitely with time.

To relax assumption #1 and test how the strength of turbulence in the mixed layer affects the bloom, we conduct experiment A3, similar to A1 except that $\kappa_{max}$ switches to a smaller value in spring ($10^{-2} \text{ m}^2 \text{ s}^{-1}$ after 1 March). This value now corresponds to a time scale of 10 d for a mixed-layer depth of 100 m. Experiment A3 mimics a situation where the mixed layer remains deep but turbulence decreases. Such a decrease has been attributed to net cooling subsidence at the end of winter (Taylor and Ferrari, 2011) and lowered windstress (Chiswell, 2011). We found that in A3, $\bar{\mu}$ becomes larger than $\bar{l}$ as soon as vertical mixing is reduced (Figure 1j). In this situation, the vertical phytoplankton distribution is no longer homogeneous in the mixed layer, but exhibits a marked surface maxima (not shown). A situation similar to this, where density was fairly homogeneous within the top 300 m and phytoplankton developed in a thin surface layer was observed by Chiswell (2011) in the bloom East of New Zealand. Experiment A3 thus breaks the CDH but satisfies the CTH: net growth becomes positive when vertical turbulence falls below a critical value. Repeating this experiment with different values of $\kappa_{max}$ would allow one to determine the critical turbulent diffusivity numerically. Note that the value of $\kappa_{max}$ that I have used here is quite high, and lower critical turbulence levels would be obtained

<table><caption>Table 3. Set of experiments.</caption><thead><tr><th></th><th><i>D</i><sub>max</sub> (m)</th><th><i>D</i><sub>min</sub> (m)</th><th><i>κ</i><sub>max</sub> (m<sup>2</sup> s<sup>-1</sup>)</th><th><i>l</i><sub>0</sub><sup>min</sup> (W m<sup>-2</sup>)</th><th><i>l</i><sub>0</sub><sup>max</sup> (W m<sup>-2</sup>)</th><th>self-shading</th><th>Grazing</th><th>I and N limitations</th><th><i>k</i><sub>b</sub> (m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>A1</td><td>300</td><td>300</td><td>1.</td><td>50.</td><td>50.</td><td>N</td><td>N</td><td>N</td><td>0.05</td></tr><tr><td>A2</td><td>220</td><td>220</td><td>1.</td><td>50.</td><td>50.</td><td>N</td><td>N</td><td>N</td><td>0.05</td></tr><tr><td>A3</td><td>300</td><td>300</td><td>1e<sup>-2</sup></td><td>50.</td><td>50.</td><td>N</td><td>N</td><td>N</td><td>0.05</td></tr><tr><td>B1</td><td>300</td><td>25</td><td>1.</td><td>50.</td><td>50.</td><td>N</td><td>N</td><td>N</td><td>0.05</td></tr><tr><td>B2</td><td>300</td><td>25</td><td>1.</td><td>50.</td><td>50.</td><td>N</td><td>N</td><td>Y</td><td>0.05</td></tr><tr><td>B3</td><td>300</td><td>25</td><td>1.</td><td>50.</td><td>50.</td><td>Y</td><td>N</td><td>Y</td><td>0.05</td></tr><tr><td>B4</td><td>300</td><td>25</td><td>1.</td><td>50.</td><td>50.</td><td>Y</td><td>Y</td><td>Y</td><td>0.05</td></tr><tr><td>C1</td><td>300</td><td>25</td><td>1.</td><td>50.</td><td>50.</td><td>Y</td><td>Y</td><td>Y</td><td>0.05</td></tr><tr><td>C2</td><td>300</td><td>25</td><td>1.</td><td>25.</td><td>200.</td><td>Y</td><td>Y</td><td>Y</td><td>0.05</td></tr><tr><td>C3</td><td>300</td><td>25</td><td>1.</td><td>25.</td><td>200.</td><td>Y</td><td>Y</td><td>Y</td><td>0.15</td></tr></tbody></table>

The first six columns are the parameters of the physical forcing (see also Figure 1), the last three are the parameters of the biogeochemical model. The A series uses Sverdrup's model (no self-shading, no grazing, no nutrient-N and light-I limitation terms). The C series uses the full NPZ model (i.e. with self-shading, grazing, and nutrients). In the B series, the model progressively transits from Sverdrup's to NPZ.
with a larger light attenuation coefficient (Huisman et al., 1999, their Fig. 4).

## Series B: Winter-to-spring transition of the mixed layer

The idea behind the CDH is that the mixed layer crosses the critical depth between winter and spring, which initiates the bloom. Implicitly, this suggests time variations in the mixed-layer depth, but in the steady-state analytical derivation of the critical depth, this effect is not explicitly accounted for. The B series of experiments focuses on a transient winter-to-spring transition: starting from unfavourable winter conditions, the mixed layer crosses the critical depth and shoals from $D_{max}$ to $D_{min}$ over time $t_s$, with $D_{min} < D_C < D_{max}$. The four experiments B1–B4 of this series share this transient physical forcing (Figure 1, central column), whereas the ecosystem model is progressively complexified from Sverdrup's model (B1) to the NPZ model (B4).

In B1, the ecosystem model is the same as in the A series (i.e. Sverdrup's model), which enables the transient and steady-state solution of the same equation to be compared. As expected, the shoaling causes $\bar{\mu}$ to increase above $\bar{l}$, in agreement with the CDH. To further assess the CDH, we compare the exact timing of the theoretical and actual bloom onsets in B1. The theoretical onset (marked by the thin vertical lines) occurs when the mixed-layer depth (Figure 1, in black) and critical depth (in red) intersect. The actual onset occurs when the growth rate (green curve) and loss rate (blue curve) intersect. Our B1 experiment reveals a time shift of ~1 week between the theoretical and actual bloom onset dates, and thus highlights a small deviation from the CDH (Figure 1e). Interestingly, this deviation does not come from the model assumptions, but from the way the original model of Sverdrup was solved, i.e. at steady state. It is worth noting that the magnitude of this deviation (1 week) is much smaller than the precision in the bloom onset detection date from satellite data (1 month, Cole et al., 2012; Brody et al., 2013). We also note that, as in the A series, the post-bloom situation in B1 is unrealistic since net phytoplankton growth is not arrested, due to the lack of nutrient limitation.

Next, in B2, B3, and B4, Sverdrup's assumptions #2–#5 regarding the form of the ecosystem model are relaxed sequentially. In B2, light saturated growth and nutrient limitation are accounted for, this is equivalent to an NPZ model. The time shift between the theoretical and actual onsets is extended by another week compared with B1 (Figure 1e). The larger shift is due to a reduced winter growth rate in response to both the nutrient limitation term (which is slightly <1) and more complex light dependence of growth (by a fairly similar amount, not shown). We also note that the post-bloom situation in B2 is more realistic than in B1, since $\bar{\mu}$ goes back to low values after the bloom, due to nutrient limitation. Assumptions #2 and #3 do, therefore, cause moderate deviation from the CDH.

Experiment B3 allows for active self-shading ($k_p$ is no longer 0), leading to time- variations in the critical depth (since $k$ increases with increasing levels of phytoplankton). We find that $k$ deviates from its background value mostly after the bloom onset, once biomass has had time to significantly accumulate. This change in $k$ is reflected in the post-bloom shoaling of the critical depth, from 270 to 80 m (Figure 1b). The shape of the growth rate is also marginally affected compared with B2 during post-bloom conditions, in response to changes in light penetration (Figure 1h). Overall, however, the bloom onset in B3 is fairly similar to that in B2, justifying Sverdrup's constant-light-extinction hypothesis (#4).

Finally, accounting for grazing (experiment B4) significantly increases the loss rate after the bloom, but not before and therefore

has little consequence on the actual onset (Figure 1k). We also note out-of-phase oscillations of the growth and loss rates in post-bloom conditions, which are typical of predator–prey interactions, with the signature of an oscillating critical depth.

A preliminary conclusion from the B set of experiments is that, at least in the chosen set of model parameters, Sverdrup's assumptions about the very simple formulations of $\mu$ and $l$ do not seriously challenge the CDH. At most, there is a time shift of 1–2 weeks between the theoretical and actual onset, and this time shift mainly ensues from bloom dynamics. Most importantly, accounting for explicit grazing (assumption #5) is not sufficient to significantly deviate from the CDH or to conform to the DRH. As shown next, this is not the case when the full seasonal is included.

## Series C: Full seasonal cycle of the mixed layer

Because of the absence of disturbance in winter, spring shoaling of the mixed layer is the only physical factor that can trigger the bloom in the A and B experiments. To address this issue, the NPZ ecosystem model is forced with a repeating seasonal cycle of the mixed layer in the C set of experiments. Starting from summer, when the mixed layer is at its shallowest ($D_{min}$), the mixed layer progressively deepens through fall and winter until it reaches $D_{max}$, then shallows more rapidly during spring (in time $t_s$) until it comes back to its initial summer value $D_{min}$ (Figure 1c). The model is integrated for 2 years so that the summer initial conditions of the second year are consistent with the physical forcing. Results from the second year are shown in the right hand column of Figure 1. Both the model and the physical forcing are identical in B4 and C1 after 1 March (i.e. after the deepest mixed layer is reached); the two experiments differ in their physical forcing during the summer-to-winter transition.

In experiment C1, phytoplankton growth and losses are in close balance during summer, when the mixed layer is stratified and the physical forcing is steady. Phytoplankton and zooplankton are confined within the mixed layer, and there is a strong vertical gradient in nutrients (Figure 2). The deepening of the mixed layer on model date November 1 creates a disturbance that drives the system out of this summer equilibrium. Nutrients are entrained within the mixed layer, and the phytoplankton and zooplankton populations are diluted. A remarkable consequence of this perturbation is the increase in phytoplankton biomass seen in November and December (solid curve, Figure 2). This increase is not seen in the surface phytoplankton concentration (dashed curve, Figure 2), which, on the contrary, decreases in response to dilution. This situation is reminiscent of the winter situation observed in the Sub-Arctic Atlantic with profiling floats (Boss and Behrenfeld, 2010). Clearly, in C1, the winter onset of net growth is associated with the dilution and decrease of the zooplankton biomass (Figure 2). Both $\bar{\mu}$ and $\bar{l}$ decrease, with the decrease in $\bar{l}$ slightly larger than that of $\bar{\mu}$ (Figure 1f). This slight imbalance between $\bar{\mu}$ and $\bar{l}$ enables the increase of the phytoplankton biomass: dilution decreases the grazing rate at a rate faster than it decreases the growth rate. Thus the winter situation in C1 satisfies the DRH: biomass starts accumulating in November, at a time when the mixed-layer depth is still above the critical depth, and thus disagrees with the CDH.

However, unlike the situation described by Behrenfeld (2010), biomass accumulation does not extend throughout winter; after January, grazing becomes negligible and the phytoplankton loss rate saturates (at the value of the mortality rate) while $\bar{\mu}$ continues to decrease, causing net growth to eventually become negative. Chiswell (2011) suggested that the CDH could perhaps be used to explain such negative net growth in winter. For experiment C1 to
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Seasonal evolution of the vertical distribution of nutrients, phytoplankton, and zooplankton in experiment C1 (colour contours). The overlaid thick black line shows the mixed-layer depth *D*. For phytoplankton and zooplankton, also shown are the seasonal evolution of the total depth-integrated biomass (solid line, top panel, in mmol m⁻²) and of the surface biomass (dashed line, top panel, in mmol m⁻³).</figcaption></figure>

be consistent with Chiswell's hypothesis, we would expect this negative net growth to be driven by the deepening of the mixed layer below the critical depth; however, the mixed layer crosses the critical depth in February, 1 month after negative net growth has initiated.

From March onward, when the mixed layer stratifies, the evolution becomes very similar to that in experiment B4 with the occurrence of a second growth period in spring, in closer agreement with the CDH. It should be noted that, as expected, the shoaling of the mixed layer induces the detrainment of phytoplankton and zooplankton biomass below the mixed layer, such that, during the onset and climax of the bloom, phytoplankton biomass is not confined to the mixed layer, as hypothesized by Sverdrup, but shows a slowly decaying profile below the mixed layer (Figure 2). We also

note that a peculiarity of C1 is that nutrients are never exhausted from the surface layer (Figure 1). This is due to the weak surface illumination $I_0$ which is held constant at a low winter value. Thus the fall and winter conditions in C1 are close to those described by Sverdrup, i.e. not nutrient limited.

In C2, we add another degree of realism by allowing $I_0$ to vary seasonally in a sinusoidal manner between winter ($I_0^{min}$, 21 December) and summer ($I_0^{max}$, 21 June) extrema. The generally larger values of $I_0$ lead to notably deeper critical depths, such that in C2 the mixed layer is almost always above the critical depth (Figure 1c). The larger surface illumination in spring and summer (compared with all previous experiments) now allows more efficient exhaustion of mixed-layer nutrients by summer (not shown). In C2, growth and loss rates are in equilibrium in summer, and both decrease, due to
the decreasing amount of light (Figure 1i). Most notably in C2, the disturbance driven by the deepening of the mixed layer on November 1 causes a rapid increase of the growth rate, in contrast to C1 where the growth rate decreased (Figure 1f). This increase ensues from the resupply of nutrients by entrainment, which alleviates nutrient limitation (not shown), and leads to a positive net growth. However, as in C1, this increase is short lived. As the mixed layer continues to deepen, light limitation increases until the date of the winter solstice (21 December) causing a rapid decay of the growth rate. In consequence, this first growth period is limited to fall (to the month of November), while in C1 it covered fall and winter (until January). Moreover, unlike in C1, this fall bloom satisfies neither the DRH nor the CDH: it is closer to a so-called entrainment bloom (i.e. Lévy et al., 2005; Martinez et al., 2011; Chiswell and Grieve, 2013).

The seasonal increase in surface light in experiment C2 drives a slow increase in the growth rate that starts at the winter solstice (21 December), i.e. much earlier than in C1 where the increase of $\bar{\mu}$ followed the mixed-layer stratification (1 March). Incidentally, this leads to a theoretical spring bloom onset in closer agreement with the CDH than in any other case: the actual bloom onset, which was late by 2 weeks in most of the previous experiments, now occurs less than a week before the theoretical bloom onset date. Accounting for variable surface irradiance thus adds a degree of freedom which counterbalances the deviations from the CDH identified before. Note that an experiment similar to C2 but with a more rapid restratification time scale $t_s$ (2 d instead of 30 d), intended to mimic the rapid decrease of turbulence in the mixed layer following the spring sign change in surface buoyancy flux (Taylor and Ferrari, 2011), has the main consequence of a more abrupt bloom onset, and even closer agreement between the actual and theoretical bloom onsets (not shown).

Finally in C3, the conditions are identical with those of C2 except that a stronger background light absorption coefficient is used ($0.15 m^{-1}$ instead of $0.05 m^{-1}$), leading to a critical depth shallower by ~100 m. The overall stronger light absorption leads to a stronger light limitation in fall which prevents the appearance of the fall bloom; in November, the growth rate decreases abruptly when the mixed layer started deepening. In this last experiment, there was only one growth period, in spring, that satisfied the CDH, although the bloom started ~2 weeks after the mixed layer became shallower than the critical depth.

Experiments C1, C2, and C3 thus illustrate three different fall/winter situations: one (C1) characteristic of the DRH, with a decrease of the loss rate that exceeds the decrease of the growth rate; one (C2) characteristic of a fall bloom, with an increase of the growth rate due to the entrainment of nutrients; and one (C3) characteristic of the CDH, with a sharp decrease of the growth rate as soon as the mixed layer became deeper than the critical depth. The spring situations, however, are more similar, with a bloom that starts shortly after the mixed layer shallowed in all three cases.

## Conclusions

The analytical expression for a critical maximum mixed-layer depth for bloom onset was formulated by Sverdrup (1953) from a phytoplankton-only model, solved at steady state. Since then, both *in situ* and satellite observations have revealed many instances where the bloom occurred before the onset of seasonal stratification, thus challenging the CDH (Townsend et al., 1994; Behrenfeld, 2010; Boss and Behrenfeld, 2010). These observations have led to a wide debate regarding Sverdrup's model assumptions (a recent review

of the debate is provided by Lindemann and St John, 2014) and to the formulation of two alternative theories that challenged some aspects of Sverdrup's hypothesis.

More precisely, the CTH recognizes that apparently mixed layers (i.e. where temperature and salinity are fairly homogeneous) are not always regions of intense vertical mixing. This situation occurs frequently at the end of winter, when the mixed layer is well established, and when cessation of surface cooling and reduction of wind stress reduces significantly vertical mixing within this layer compared with winter values. On the other end, the DRH challenges Sverdrup's assumption of a constant loss (mortality) rate. In particular, Behrenfeld (2010) suggests that a decrease in the loss rate is as plausible as an increase in the growth rate in initiating the bloom. In winter, such a decrease could be initiated by the dilution of predators and prey when the mixed layer entrains plankton-free waters.

Conceptual models of the bloom now incorporate all three hypotheses (Chiswell, 2011; Lindemann and St John, 2014). In parallel, numerical models have shown their ability to simulate them all (e.g. Dutkiewicz et al., 2001; Behrenfeld et al., 2013). It is uncertain which is more likely occurring in nature, but understanding is critical if one wants to project how phytoplankton dynamics will evolve in a changing ocean (Lindemann and St John, 2014).

Beyond all debate, Sverdrup's approach was indisputably pedagogical: it illustrated for the first time how a basic ecosystem model could be used to address a key ecological problem, namely the onset of phytoplankton blooms. In this regard, building on Sverdrup's approach, I presented a unifying framework in which the different concepts could be tested and rationalized. The framework is a one-dimensional NPZ model, forced by vertical mixing in a time-varying mixed layer. It builds on the initial phytoplankton model of Sverdrup's and adds the ingredients necessary for testing the CTH and DRH, i.e. explicit vertical mixing with variable strength $\kappa$ and variable loss through incorporation of grazing. In addition, this framework also tested other assumptions of Sverdrup, i.e. steady state, constant-light absorption, the absence of nutrient limitation, and a productivity proportional to light.

The set of experiments presented here demonstrates that all three mechanisms for stimulating phytoplankton blooms can be modelled within this 1D-NPZ framework. It showed that adding complexity to the ecosystem model is not sufficient to make the model solution deviate from the CDH: the understanding of bloom dynamics requires a comprehensive representation of the physical drivers, and this representation was oversimplified in Sverdrup's model. The following lessons have been learned, and it is to be hoped will be useful for future bloom studies.

First, it is crucial to make a clear distinction between the mixed layer and the mixing layer, particularly when analysing observations; a number of prior studies have pointed that the seasonal thermocline (defined on a density criteria) is not always a good proxy of the mixing layer (Franks, 2015). This distinction is probably more natural to modellers, who can easily access to the model vertical mixing coefficient. Observational estimates of vertical mixing are sparse, though it may be possible to use the large array of Argo floats, from which the mixed-layer depth can be extracted, in conjunction with satellite observations of the bloom (Sallée et al., 2015). An interesting alternative is to compare the bloom onset with the surface net heat flux rather than to a critical depth (Ferrari et al., 2014).

Second, it is crucial to account for the evolution of the physical parameters from the previous summer—and thus of the full seasonal cycle—to understand bloom dynamics. Many of the recent studies on
the bloom onset focus on the winter-to-spring transition (Brody and Lozier, 2014; Ferrari et al., 2014). However, as nicely illustrated by the conceptual model of Chiswell (2011), the summer-to-winter transition can also be associated with the accumulation of phytoplankton biomass. These early onsets may arise from a diverse set of processes, such as entrainment of nutrient (Chiswell, 2011), dilution of grazers (Behrenfeld, 2010), or increasing levels of solar radiation (Llort et al., 2015). In the experiments presented here, the fall/winter production period is clearly separated from the spring bloom by a period of negative net growth in winter. However, the two periods can also merge in a single, long period of net growth and this has been reported both from observations (Lévy et al., 2005; Behrenfeld, 2010; Sallée et al., 2015) and from different sets of model experiments (Karleskind et al., 2011; Behrenfeld et al., 2013; Llort et al., 2015).

Third, the seasonal evolution of surface irradiance is a key, and often neglected ingredient of the bloom timing. It is puzzling that in our most realistic experiments, when seasonal variations of solar radiation are accounted for, the actual onset is in closer agreement with Sverdrup's theoretical onset date, whereas these variations were not accounted for in the derivation of the critical depth. We find that this good agreement is due to counter-balancing effects, where solar radiation shifts the bloom onset toward the winter solstice, i.e. before the theoretical CDH date, while unsteady bloom dynamics shift the onset past this theoretical date.

Finally, variable loss—through incorporation of grazing—is necessary but not sufficient to conform to the DRH. In fact in all of our experiments, the loss rate is constant from winter to spring, in agreement with Sverdrup's assumption that grazers play a minor role on the onset of blooms in spring. Winter bloom conforming to the DRH only emerges in the specific situation where nutrients are still available at the end of summer and light levels are low, a situation that can be found at subpolar latitudes. Within the parameter range explored here, this winter onset disappears as soon as the seasonal evolution of surface irradiance is accounted for.

To finish, I remind the reader that the conclusions drawn from this illustrative set of experiments are limited by the fact that neither the physical nor the biological parameter spaces have been exhaustively explored, and that the role of other potentially important mechanisms such as acclimatation of the cells to external conditions or the physiology of sinking (Lindemann and St John, 2014) has not been addressed.

## Acknowledgements

I thank the Center for Atmosphere Ocean Science, Courant Institute of Mathematical Sciences, New York University, and in particular Shafer Smith, for inviting me to give lectures in their department, which provided the inspiration for this note. Many thanks also to J. Llort, J.-B. Sallee, L. Bopp, J.-M. Andre, and A. Tagliabue for their insightful comments, to R. Ferrari for lively discussions on the bloom onset and to three anonymous reviewers for their helpful comments.

## References

* Banse, K. 1994. Grazing and zooplankton production as key controls of phytoplankton production in the open ocean. Oceanography, 7: 13–20.

* Behrenfeld, M. 2010. Abandoning Sverdrup's Critical Depth Hypothesis on phytoplankton blooms. Ecology, 91: 977–989.

* Behrenfeld, M. J., and Boss, E. S. 2014. Resurrecting the ecological underpinnings of ocean plankton blooms. Annual Review of Marine Science, 6: 167–194.

* Behrenfeld, M. J., Doney, S. C., Lima, I., Boss, E. S., and Siegel, D. A. 2013. Annual cycles of ecological disturbance and recovery underlying the subarctic Atlantic spring plankton bloom. Global Biogeochemical Cycles, 27: 526–540.

* Boss, E., and Behrenfeld, M. 2010. In situ evaluation of the initiation of the North Atlantic phytoplankton bloom. Geophysical Research Letters, 37: L18603.

* Brainerd, K. E., and Gregg, M. C. 1995. Surface mixed and mixing layer depths. Deep Sea Research Part I, 42: 1521–1543.

* Brody, S. R., and Lozier, M. S. 2014. Changes in dominant mixing length scales as a driver of subpolar phytoplankton bloom initiation in the North Atlantic. Geophysical Research Letters, 41, doi:10.1002/2014GL059707.

* Brody, S. R., Lozier, M. S., and Dunne, J. P. 2013. A comparison of methods to determine phytoplankton bloom initiation. Journal of Geophysical Research-Oceans, 118: 2345–2357.

* Chiswell, S. 2011. Annual cycles and spring blooms in phytoplankton: don't abandon Sverdrup completely. Marine Ecology-Progress Series, 443: 39–50.

* Chiswell, S. M. 2013. Comment on "Annual cycles of ecological disturbance and recovery underlying the subarctic Atlantic spring plankton bloom". Global Biogeochemical Cycles, doi: 10.1002/2013GB004681.

* Chiswell, S. M., and Grieve, J. B. 2013. Climatology of surface chlorophyll a, autumn-winter and spring blooms in the south-west Pacific Ocean. Journal of Geophysical Research: Oceans, 118: 1003–1018.

* Cole, H., Henson, S., Martin, A., and Yool, A. 2012. Mind the gap: the impact of missing data on the calculation of phytoplankton physiology metrics. Journal of Geophysical Research, 117: C08030.

* Dutkiewicz, S., Follows, M., Marshall, J., and Gregg, W. 2001. Interannual variability of phytoplankton abundances in the North Atlantic. Journal of Geophysical Research, 48: 2323–2344.

* Evans, T. G., and Parslow, J. S. 1985. A model of annual plankton cycles. Biological Oceanography, 3: 327–347.

* Ferrari, R., Merrifield, S., and Taylor, J. 2014. Shutdown of convection triggers increase of surface chlorophyll. Journal of Marine Systems, S0924-7963(14)00038-4, doi: 10.1016/j.jmarsys.2014.02.009.

* Franks, P. J. 2002. NPZ models of plankton dynamics: their construction, coupling to physics, and application. Journal of Oceanography, 58: 379–387.

* Franks, P. J. 2015. Has Sverdrup's critical depth hypothesis been tested? Mixed layers vs. turbulent layer. ICES Journal of Marine Science Journal, 72: 1897–1907.

* Huisman, J., van Oostveen, P., and Weissing, F. J. 1999. Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms. Limnology and Oceanography, 44: 1781–1787.

* Karleskind, P., Lévy, M., and Memery, L. 2011. Subduction of carbon, nitrogen, and oxygen in the northeast Atlantic. Journal of Geophysical Research, 116: C02025.

* Lévy, M., Lehahn, Y., André, J., Mémery, L., Loisel, H., and Heifetz, E. 2005. Production regimes in the Northeast Atlantic: a study based on Sea-viewing Wide Field-of-view Sensor (SeaWiFS) chlorophyll and ocean general circulation model mixed layer depth. Journal of Geophysical Research, 110: 1998–2006.

* Lindemann, C., and St John, M. A. 2014. A seasonal diary of phytoplankton in the North Atlantic. Frontiers in Marine Science, 1: 1–6.

* Llort, J., Lévy, M., Sallée, J-B., and Tagliabue, A. 2015. Onset, intensification, and decline of phytoplankton blooms in the Southern Ocean. ICES Journal of Marine Science, 72: 1971–1984.

* Marra, J., and Barber, R. T. 2005. Primary productivity in the Arabian Sea: a synthesis of JGOFS data. Progress in Oceanography, 65: 159–175.
Martin, J. H., Knauer, G. A., Karl, D. M., and Broenkow, W. W. 1987.
 VERTEX: carbon cycling in the northeast Pacific. Deep Sea Res.,
34: 267-285.

Martinez, E., Antoine, D., D'Ortenzio, F., and de Boyer Montégut, C.
2011. Phytoplankton spring and fall blooms in the North Atlantic
in the 1980s and 2000s. Journal of Geophysical Research, 116: C11029.

Riley, G. A. 1946. Factors controlling phytoplankton populations on Georges Bank. Journal of Marine Research, 6: 54–73.

Sallée, J. B., Llort, J., Tagliabue, A., and Lévy, M. 2015. Characterisation
of distinct bloom phenology regimes in the Southern Ocean. ICES
Journal of Marine Science, 72: 1985–1998.

Siegel, D., Doney, S., and Yoder, J. 2002. The North Atlantic spring phytoplankton bloom and Sverdrup's critical depth hypothesis. Science, 296: 730.

Sverdrup, H. 1953. On conditions for the vernal blooming of phyto-
plankton. Journal du Conseil / Conseil Permanent International
pour l'Exploration de la Mer, 18: 287-295.

Taylor, J. R., and Ferrari, R. 2011. Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms. Limnology and Oceanography, 56: 2293–2307.

Townsend, D., Cammen, L., Holligan, P., Campbell, D., and Pettigrew, N. 1994. Causes and consequences of variability in the timing of spring phytoplankton blooms. Deep Sea Research Part I, 41: 747.

Yoshie, N., Yamanaka, Y., Kishi, M. J., and Saito, H. 2003. One dimensional ecosystem model simulation of the effects of vertical dilution by the winter mixing on the spring diatom bloom. Journal of Oceanography, 59: 563–571.

Handling editor: Rubao Ji