Onset, intensiﬁcation, and decline of phytoplankton blooms in the Southern Ocean

TheseasonalcycleofphytoplanktonbiomassintheSouthernOcean(SO)ischaracterizedbyaperiodofrapidaccumulation,knownasbloom,thatistypicalofhigh-latituderegions.RecentstudieshaveillustratedhowspatialandtemporaldynamicsofbloomsintheSOaremorecomplexthaninotheroceans.Thiscomplexityislikelyrelatedtodifferencesinverticalmixingandtheironavailability.Inthiswork,weexaminethesensitivityofbloomdynamicstochangesinverticalmixingandironavailabilityusingabiogeochemicalmodel.Underidealizedphysicalforcing,weproduceseasonalcyclesofphytoplanktonforanensembleofSOscenariosandwedescribethebloomdynamicsintermsofthenetbiomassaccumulationrate.Basedonthismetric,wedeﬁnethreecrucialbloomphases:theonset,theclimax,andtheapex.Fortheensembleofmodelledblooms,onsets always occur in winter and can be either bottom-up (increase in productivity) or top-down (decrease in grazing) controlled. Climaxes are mostly found in spring and their magnitudes are bottom-up controlled. Apexes are always found in late spring and strongly top-down controlled. Our results show that while a “strict” onset deﬁnition is consistent with a winter onset, the surface spring bloom is associated with the climax of the integrated bloom. Furthermore, we demonstrate that onset phase can be distinguished from climax phase using appropriate bloom detection methodsbasedonsurfacesatellite-basedproducts.TheensembleoftheseresultssuggeststhatSverdrup’sbloomingconditionsarenotindicativeofthebloomonsetbutoftheclimax.Weconcludethattherecentbloomonsetdebatemaypartlybeduetoaconfusionbetweenwhatisdeﬁnedhereasthebloomonsetandtheclimax,andthattheSOobservedcomplexityisduetothefactorsthatcontroltheclimax.


Introduction
The Southern Ocean (SO) is the largest high-nutrient lowchlorophyll region in the world's ocean. Its relatively low productivity has been attributed to a combination of iron scarcity (Martin et al., 1990), elevated grazing, and light limitation (Boyd and Ellwood, 2010). Despite these unfavourable biological growth conditions, large accumulations of phytoplankton biomass, or blooms, are observed in surface waters each spring over wide areas of the SO (Thomalla et al., 2011). The distribution of these blooms is very patchy in space and time. Hot-spots of phytoplankton accumulation are mainly seen where sources of iron are significant, i.e. in the lee of Islands (Moore and Abbott, 2000;Arrigo et al., 2008). Additionally, the bloom onset dates are rather spread in time, from October to January (Thomalla et al., 2011), and do not show a clear latitudinal pattern. This is unlike the pattern for North Atlantic spring bloom which is zonal and propagates from South to North (Siegel, 2002;Lévy et al., 2005). The variability of bloom dynamics in the SO in terms of their amplitude, timing and location has been mainly documented from ocean-colour observations (Moore and Abbott, 2000;Arrigo et al., 2008;Thomalla et al., 2011). However, the drivers of the observed variability remains unclear. We hypothesize that patchy environmental conditions, involving zonally asymmetric mixed-layer distributions (Sallée et al., 2010) combined with equally complex dissolved iron distributions (Tagliabue et al., , 2014 are mostly responsible for the complex SO patterns. Our current understanding of phytoplankton bloom dynamics mainly comes from works based on the North Atlantic. In this region, the mixed-layer, nutrient, and atmospheric environment are largely different from in the southern hemisphere. Historically, the emergence of blooms in the North Atlantic has been related to thinning of the ocean surface layer where turbulence is active (hereafter referred to as "turbulent or mixing-layer", or XLD, to be distinguished from the usual "mixed-layer depth", or MLD, which represents the upper-layer where hydrographical properties are well mixed; Franks, 2015). This thinning implied an increase of averaged exposure of phytoplankton cells to light (Gran and Braarud, 1935;Riley, 1942). Along with this bottom-up view, Sverdrup (1953) proposed that bloom would start when surface mixing-layer crosses a critical depth above which integrated phytoplankton growth would overcome phytoplankton losses (Siegel, 2002). It must be emphasized to note that Sverdrup (1953)'s hypothesis (also known as the Critical Depth hypothesis) is founded on the assumption that "within the top layer the turbulence is strong enough to distribute the plankton organisms evenly through the layer " (Sverdrup, 1953). This assumption is crucial to understand the Critical Depth hypothesis as points out that what matters for phytoplankton is the vertical profile of turbulent mixing, rather than the hydrographical properties of the water column. Indeed, the relevant parameter for the Sverdrup (1953)'s hypothesis is the XLD rather than the MLD (Franks, 2015).
As an alternative to the bottom-up understanding of ocean blooms, a top-down view has also appeared. This view proposes that, in essence, the causes of phytoplankton concentrations cannot be fully understood without considering the role of their main predator, zooplankton (Banse, 1994). The top-down hypothesis has gained interest recently via a series of papers that challenged the prevailing "bottom-up" paradigm . Using various tools (satellite data: Behrenfeld, 2010; float data: Boss and Behrenfeld, 2010; and model estimates: , it has been suggested that the North Atlantic spring bloom does not initiate in spring alongside thinning mixing-layers, but rather in winter when mixing-layer is deepening. This winter initiation is consistent with the hypothesis that dilution enables phytoplankton to better escape their predators and accumulate biomass (Evans and Parslow, 1985;Yoshie et al., 2003;Marra and Barber, 2004).
In this context, our primary objective is to examine the drivers of phytoplankton blooms over the full range of SO environmental conditions. In particular, we examine how different environmental conditions (mixing-layer depth, ferricline, and solar radiation) result in bottom-up or top-down control. To that end, we extend a framework in which the rate of net biomass accumulation (r) results from a competition between growth (m) and loss (l ) of phytoplankton: r = m − l = (1/P)(dP/dt), with P the total biomass of phytoplankton present in the water column (Riley, 1942;Sverdrup, 1953;Behrenfeld, 2010). At equilibrium, phytoplankton growth and loss are in balance and phytoplankton population remains stable. This balance can be disturbed by a sudden change in iron supply, light conditions, or stratification in isolation or in combination, which would then likely affect m and l in different ways . Our overarching question is how such perturbations modify m and l at seasonal scale, and which term is the most sensitive to a given perturbation and the most effective at driving variations in r (i.e. phytoplankton population fluctuations). One of the key aspects of our approach is that we examine three important phases in the annual cycle of blooms: the bloom onset, climax, and apex. Using time evolution of r, we define these bloom phases as follows: 1. The bloom onset is when total biomass starts to accumulate.
2. The bloom climax occurs when the rate of biomass accumulation is maximal. After this instant, accumulation of phytoplankton continues but at a slower rate because ecosystem has yet started its way to readjustment (i.e. to recoupling).
3. Finally, the bloom apex marks the peak in total biomass, i.e. the time after which l . m (recoupling is actually achieved) and accumulation starts decreasing. Figure 1 illustrates these three phases associated with the time of minimum integrated biomass (onset), of maximum slope of integrated biomass (climax), and of maximum integrated biomass (apex). This distinction complements previous studies on bloom dynamics that focused either exclusively on onset (Sverdrup, 1953;Behrenfeld, 2010) or on climax (Lozier et al., 2011;Ferrari et al., 2014). We address the question of bloom drivers in the SO in the framework of a numerical model. This model uses a state-of-the-art biogeochemical network (Aumont and Bopp, 2006) within a vertically discretized 1D water column configuration where vertical mixing is the main physical process. Aiming to examine the full range of SO conditions, we perform statistical analyses using an ensemble of 1200 model simulations with distinct seasonal cycles of mixing-layer depth, ferricline, and solar radiation. Bloom onset, climax, and apex dates are diagnosed for each run in the ensemble. Simultaneously, the distributions of variables such as iron supply, mixing-layer depth, light, phytoplankton growth, and loss rates for each of the three bloom phases are attributed to drivers. We investigate when bottom-up or top-down controls prevailed. This model allows us to test existing theories on bloom onset in an idealized and comprehensive framework, and to discuss them in the context of theSO. Furthermore, the completeness of the model data allows us to compare the onset and climax dates with the dates at which two different satellite-based bloom detection methods identify bloom initiation. These bloom detection methods are also compared with the date of the bloom onset predicted by Sverdrup (1953)'s hypothesis with the aim to evaluate at which point the validation (or refusal) of this hypothesis is influenced by the bloom detection method.

Biogeochemical model
The model was set up to represent the Permanent Open Ocean Zone (POOZ) of the SO, away from ice formation and melting, where nitrate and silica do not limit productivity. Our goal in this study is to untangle how the different phases of a bloom (onset, climax, and apex) are controlled by their physical and biogeochemical environment. As such, we deliberately chose to reduce the complexity of the problem by considering a 1D physical framework (e.g. lateral advection is neglected). Varying vertical diffusion reproduces seasonal cycle of the mixing-layer depth. Along with this idealized physical configuration, we model the associated biogeochemical activity with the model PISCES (Aumont and Bopp, 2006). PISCES contains 24 biogeochemical tracers with five nutrients able to limit phytoplankton growth: nitrate, phosphate, ammonium, iron, and silicate. The iron pool is explicitly modelled and controlled by a range of processes such as phytoplankton uptake, bacterial uptake, zooplankton, and bacterial recycling, remineralization and scavenging. In addition, four living pools are represented: two phytoplankton size classes (small and large) and two zooplankton size classes (microzooplankton and mesozooplankton). Large phytoplankton differs from the small phytoplankton by higher requirements in iron and a greater iron half-saturation constant. Grazing pressure on each phytoplankton is also differentiated by size: microzooplankton (Z) preferentially grazes small phytoplankton while mesozooplankton (M) preferentially grazes the larger phytoplankton but also microzooplankton.
Prognostic equation for each phytoplankton group (i ¼ 1,2) is: where P i is the phytoplankton biomass, m i is the growth rate, g i represents the grazing rate, and m i is the mortality rate. The last right hand side term is the effect of vertical diffusion over biomass due to vertical mixing of intensity k z . The growth rate (m i ) is computed as follows: where f i (T ) is the dependence of the growth rate with temperature (Eppley, 1972), h(z) is a penalization term for deep mixing, PAR i is function of the shortwave radiation at the surface, a i is the initial slope of the photosynthesis-irradiance curve, Q Chl i the Chl:C quota for each phytoplankton, and L i is the nutrient limitation. Nutrients limitation in our set-up can only be due to Fe, so Iron limitation is formulated following a Quota approach (McCarthy, 1980;Droop, 1983) with the term Q Fe i,opt allowing luxury uptake (as in Buitenhuis, 2010). Grazing rate by zooplankton is computed as a Michaelis -Menton parametrization with a phytoplankon feeding threshold (Aumont and Bopp, 2006). It can be formulated as follows: where Z stands for both micro-and mesozooplankton, g 0,Z m is the maximum grazing rate at 08C, f Z (T) is a temperature-dependent function, K Z graz is the corresponding half-saturation constant, g Z Pi is the preference of each zooplankton for each phytoplankton, and F Z,Pi thresh is a function computing the corresponding feeding thresholds Onset, intensification, and decline for each couple phyto/zoo. Encounter probability between grazers and phytoplankton is not explicitly computed but results from the volumetric concentration of the four biological compartments. The model equations were computed on a regular vertical grid of 74 + 1 vertical levels (constant spacing of 7 m for the first 74 levels and a last one of 500 m depth) and a time step of 20 min. The vertical mixing coefficient (k z ), temperature and surface photosynthetic available radiation (PAR) were analytically prescribed every 6 h. Vertical mixing coefficient were assumed to be constant and equal to 1 m 2 s 21 within a surface mixing layer of depth XLD, and equal to 10 25 m 2 s 21 below. Hence, in this study, XLD was not calculated but imposed through the k z vertical profile. The extremely large value of k z within the XLD guaranteed that turbulence was strong enough to homogeneously distribute phytoplankton (i.e. Sverdrup, 1953s assumption). Following Lévy (2015), who highlighted the need to represent the full seasonal cycle of the XLD to study the bloom, we imposed an idealized seasonal cycle of the XLD divided into three phases (red curve in Figure 1): † A fall/winter phase of convection and progressive XLD deepening. † A spring phase where cessation of convection led to the thinning of the XLD. † A summer phase with a relatively shallow and constant XLD.
To ensure that the timing of these phases and the magnitude of the XLD were relatively realistic, we used data estimates derived from Argo data (Sallée et al., 2010). These data provided us with an estimate of the depth of the seasonal thermocline in the SO, which we assumed to reflect the mixing depth. Sub-seasonal variability in the XLD was not accounted for. For temperature and surface PAR, we used a smoothed climatological seasonal cycle constructed from observations (DFS3-ERA40; Broderau et al., 2008) averaged over the 40 -608S latitudinal band. The summer initial condition for dissolved Fe profile was constructed by assuming low concentrations (0.03 nMolFe/l) above a prescribed ferricline depth (Z Fe ), and larger concentrations (0.5 nMolFe/l) below. The depth of the ferricline for the initial condition iron profile is one of the parameters we varied in our set of simulations. It is generally understood that the summer Fe profile is set by a combination of remineralization, scavenging and physical supplies by lateral sources. While remineralization and scavenging are parametrized in PISCES, there remains a large degree of uncertainty in the parameterization. In addition, lateral supplies were not explicitly accounted for. To overcome these issues and to allow the model to reach a repeating and stationary seasonal cycle, the dissolved iron profile was restored towards its initial value at the end of each summer.
Initial vertical profiles for macronutrient (i.e. nitrates, phosphates, and silicates) were constructed based on the winter mean profiles collected during the KERFIX project (Jeandel et al., 1998). This project aimed to monitor ocean-atmosphere CO 2 and O 2 exchanges and related processes with a time-series station (called KERFIX station) located at 50840 ′ S-68825 ′ E, 60 miles southwest of the Kerguelen Islands (SO). From January 1990 to March 1995, regular monthly measurements of physical and biogeochemical water properties were carried out at KERFIX station. As for iron, macronutrients were restored to the initial profiles at the end of each summer. Initial conditions for the four living compartments were set to low values for the first year. The simulations were integrated for 3 years, starting in austral summer (15 February), with outputs saved at daily frequency. A repeating seasonal cycle was reached after 2 years and results are based on the third year of simulation. As an example, Figure 1 shows a complete seasonal cycle of integrated biomass (black curve), surface Chl (green), XLD (red), and summer Z Fe , for one of these runs.

Ensemble runs
An ensemble of runs was performed by varying the XLD and summer Z Fe in the range of values found in the POOZ. Specifically, these XLD and Z Fe were modified based on the following variables: (i) The winter maximal mixing depth (XLD max ) (ii) The summer minimal mixing depth (XLD min ) (iii) The date at which XLD max was reached (tXLD max ) (iv) The date at which XLD min was reached (tXLD min ) The variables i -v were set based on a discrete set of observed values, with equal weight given to each discrete value (see ranges in Figure 1, blue arrows). The ranges for Z Fe were established based on a recent compilation of dissolved iron measurements , the ranges of values used to set XLD were based on .500 000 density profiles sampled in the SO by Argo floats data (Sallée et al., 2010). Our choice of variables led to an ensemble of almost 1200 different scenarios that combine different values of the above parameters i -v covering a wide range of XLD and Z Fe observed in the SO.
In our model, the amount of Fe injected at the surface each year was not prescribed: Fe was entrained in the mixing layer during the deepening phase. Thus, the Fe supply depended on the relative depths of the ferricline and XLD max ( Figure 2). The relationship between both variables was however non-linear, due to the effects of consumption/remineralization by the biological community and the rate of stratification/destratification. A peculiarity of the SO is that the ferricline depth is often found below the maximum  winter mixed-layer depth (Tagliabue et al., 2014), implying some regions permanently Fe-limited despite strong and deep winter mixing. In our scenarios, this situation occurred when Z Fe was greater than XLD max ( Figure 2). Note that this is different from other high-latitude productive regions such as the North Atlantic, where nitrate is the main limiting nutrient. In the North Atlantic, the depth of winter mixing and the convective nitrate supplies are more tightly correlated than its SO counterpart (depth of winter mixing and convective Fe supplies).

Bloom onset, climax, and apex
The bloom phenology was decomposed into three main events defined by the rate of net biomass accumulation (r), which writes as: where H is the depth of the water column. For simplicity, hereinafter we will refer to P int as P. From the seasonal evolution of r, we defined the bloom onset (total biomass starts to increase: P min , r = 0, and r ′ . 0; hereafter all temporal derivatives are marked by a prime, i.e. r ′ = ∂r/∂t); the bloom climax (the rate of accumulation is maximum: r = r max and r ′ = 0) and the bloom apex (bloom peaks in total biomass: P = P max , r = 0 and r ′ , 0). See Figure 1 where the three steps are reported. We note here that r was computed in our study from a total water column integral, which slightly differs from what is done in Sverdrup (1953) or Behrenfeld (2010) where P was integrated down to the base of the mixing layer. As Chiswell (2013) pointed out, mixing layer integration of P might be misleading when the mixing layer restratifies as plankton is not conserved in the XLD. Integrating over the whole water column overcomes the discontinuity issue pointed out by Chiswell (2013).
Integrating Equation (1) and dividing all terms by P, r can be written as the integrated balance between phytoplankton source (i.e. growth rate, m) and sinks (i.e. grazing and mortality rates, g and m, respectively): Hence, the evolution of modelled ecosystem in the water column can be synthesized as: where m is the mean growth rate of the total depth integrated phytoplankton community, and l is the sum of grazing and mortality.

Bloom timing
In our set of experiments, the time of the deepest convection (tXLD max ) can vary by up to 2 months between experiments, and the time at which summer stratification is reached (tXLD min ) by 1 month (Figure 1). To account for this variability, the timing of the different bloom phases was not only measured relatively to the day of the year but also relatively to the phase of the physical forcing (i.e. relative to the time of tXLD min and tXLD max ). In this sense, bloom phases before tXLD max occurred in "winter", when the mixing layer is still deepening. Similarly, bloom phases occurring between tXLD max and tXLD min are occurred during the spring thinning of the mixing layer.

Bottom-up vs. top-down control
In this paper, we aim at investigating whether the bloom seasonal cycle is controlled by bottom-up or top-down processes. Here, we detail how the relative intensity of m ′ and l ′ can be used to link onset, climax, and apex to their bottom-up or top-down controls.
Onset (r ¼ 0, r ′ . 0) occurs when gains first overcome losses. At onset integrated phytoplankton biomass is minimal and losses are always decreasing. Under these circumstances, two possible mechanisms can cause the bloom onset ( Figure 3): 1. Growth regime: The growth rate has already started to increase while the loss rate is still decreasing or stable. The growth will then become larger than loss at some point leading to the initiation of net biomass accumulation in the water column (i.e. r . 0). The onset is controlled by the growth (i.e. light and nutrients) and therefore is bottom-up driven. Analytically, this regime can be expressed as: 2. Dilution regime: Growth and loss rate are decreasing due to nutrients depletion, low light conditions, and XLD deepening. The latter causes the dilution of plankton (both, phyto-and grazers) when increasing the volume of water in the mixing layer. This process strongly decreases the prey -predator encounter probability causing a faster decrease on loss rate than on growth rate. The grazing pressure relaxation allows phytoplankton ecosystem to start increasing. This regime is top-down controlled and it corresponds to bloom onset scenario described in Marra and Barber (2004), Behrenfeld (2010), and Boss and Behrenfeld (2010). It can be analytically expressed by Climax (r = r max , r ′ = 0) marks the instant of the fastest population increase. Consequently, climax is the inflection point in the seasonal evolution of biomass: From an ecosystem point of view, this means that the trend in loss rate overcomes the trend in growth rate (l ′ ≥ m ′ ). Therefore, the Onset, intensification, and decline bloom climax marks the beginning of the recoupling process that leads to the re-equilibrium of the system. Climax can be achieved in two different ways: either m ′ becoming negative due to nutrient limitation (bottom-up control) or l ′ becoming greater than m ′ (top-down control). Finally, apex (r ¼ 0, r ′ , 0) marks the actual time of recoupling, when losses first equals gains. At apex, losses are always increasing l ′ . 0. Apex can be reached when growth is still increasing, and biomass accumulation is stopped by grazers (i.e. l ′ . m ′ with m ′ . 0); we refer to this case as the top-down controlled. In the bottom-up case, the apex is mainly due to change on the growth rate trend (i.e. m ′ , 0). This situation is often caused by the nutrients depletion in the mixing layer. We note that both top-down and bottom-up controls can mutually act together. Our analysis only points out the dominant process at play.

Bloom onset detection and Sverdrup's hypothesis
With the aim to evaluate how onset and climax phases can be detected using ocean-colour data, we applied two different bloom detection methods to model outputs. Several bloom detection methods exist in the literature and it has been shown that bloom detection dates can strongly differ depending on which method is applied (Ji et al., 2010;Brody et al., 2013). Here we use two methods, a surface biomass-(sP) and a surface chlorophyll-(sChl) based methods, that have already been implemented in literature with ocean-colour data (Behrenfeld, 2010;Brody et al., 2013). Despite these methods are designed to be applied with ocean-colour data, in our case bloom detection dates are obtained from modelled sP and sChl and compared to the actual onset and climax computed from the full vertical profile. The two bloom detection methods are defined as follows: 1. P*-method: date at which P * ′ . 0 Behrenfeld (2010), where: 2. sChl-method: Date of maximal sChl" (Sallée et al., 2015) The P*-method is based on a depth integrated view of the bloom: it estimates the amount of biomass within the water column assuming that phytoplankton is homogeneously mixed in the ocean upper layer and that the amount of biomass below is negligible. Onset is then detected when integrated biomass starts to increase. On the other hand, the sChl-method is only based on the surface imprint of the bloom and onset detection is based on the rate of change of sChl. Interestingly, when used in the literature, the P*-method resulted on bloom onset detected in winter (Behrenfeld, 2010) while sChl-method detects bloom onsets mostly in spring (White et al., 2009;Sallée et al., 2015). A number of studies have addressed high-latitude blooms using ocean-colour data with the aim to validate (or to reject) Sverdrup (1953)'s hypothesis (Siegel, 2002;Behrenfeld, 2010;Chiswell, 2011). Here, we aim to quantify at which point the validation (or refusal) of the Sverdrup (1953)'s hypothesis depends on the bloom detection method implemented. We took advantage of the model data completeness to compute the critical depth (Z c ) based on the two main Sverdrup (1953)'s assumptions: strongly turbulent surface layer and constant mortality (assumed to be the main loss term in winter-early spring). The formal expression of Sverdrup (1953) critical depth (as presented in Lévy, 2015) aI 0 kZ c (1 − e −kZc ) = m (12) can be approximated to where m 0 is the phytoplankton growth rate at surface (i.e. m 0 = aI 0 ), k is the light attenuation coefficient (in m 21 ), and m is the mortality rate. The variable Z c is computed at each time step using the model outputs for m 0 , m averaged from August to October and a light attenuation coefficient of k ¼ 0.05 m 21 , constant throughout the year (i.e. no phytoplankton shelf-shading; Lévy, 2015). Determining the date at which XLD = Z c for each modelled bloom, we are able to investigate whether Sverdrup's bloom conditions are satisfied or not for the two surface bloom detection methods presented above (P*-method and sChl-method).

Abrupt and smooth blooms
Two types of bloom phenology emerge from our 1200 runs ensemble: abrupt blooms, characterized by a sudden and very strong intensification of biomass accumulation, and smooth blooms, which display a smoother biomass accumulation. In fact, there is a continuous range of possible phenologies between these two bloom types and, hence, no objective method that distinguishes them. Nevertheless, as abrupt blooms reach, by definition, a higher value of r at climax, the 20% of blooms with the largest r max were identified as abrupt, and the remaining 80% as smooth in the following analysis. Importantly, these two types of SO bloom phenologies are also identifiable from ocean-colour observations (Sallée et al., 2015). For illustrative purpose, we will describe an example of an abrupt and smooth bloom taken from our results (Figure 4). In the abrupt case example (Figure 4a, c, and e), the XLD reaches 400 m in winter and the summer ferricline is located at 150 m. This deep winter mixing causes strong light limitation over 6 months of the year (from May to October; yellow surface in Figure 4e). Simultaneously, Fe limitation (red surface in Figure 4e) declines as XLD becomes deeper than Z Fe and entrains Fe to the surface. The bloom onset occurs around 1 July and is followed by 2.5 months of a low and stable positive r during which the XLD continues to progressively deepen. We will refer to this period as the plateau. Then, when tXLD max is reached and the XLD starts shallowing, r rapidly increases until climax (r max ) is reached on 20 October. In this scenario, the climax is an abrupt and strong peak that occurs during the period of XLD restratification. This date also marks the start of a large and rapid increase in both integrated biomass (P) and surface chlorophyll (sChl) (black and green lines in Figure 4a). Apex is reached 10 days after the climax (1 November) associated to a rapid decline in r, which is driven by decreasing m as Fe-limitation becomes important, as well as increasing losses (l ′ ≫ 0). Following apex, growth/loss equilibrium (r ≈ 0) is re-established over the summer (i.e. grazers-prey recoupling). This type of bloom is characteristic of high-latitude regions like the North Atlantic and iron rich waters of the SO (Waniek, 2003).
In the smooth bloom case example (Figure 4b, d, and f), XLD reaches 200 m in winter, with a long period of restratification (.2 months) and a relatively deep mixing during summer (65 m). The variable Z Fe is deeper than XLD max , so mixing does not reach the largest Fe stocks, which maintains a significant Fe limitation all yearround. As for abrupt bloom case, onset occurs in mid-June and is followed by a plateau phase which lasts during XLD deepening. In contrast to the abrupt bloom (Figure 4a, c, and e), the plateau has a much smoother shape and, instead of switching to a high r phase, is followed by a phase where r declines slowly. Climax is in late-July, in the midst of the plateau and just before tXLD max ; the accumulation intensity at climax is almost four times lower than for the abrupt bloom case (r smooth max = 0.02 d −1 compared with r abrupt max = 0.075 d −1 ). Restratification is not associated with an increase in accumulation, which is likely due to prevailing Fe limitation. The climax lead to the beginning of recoupling and apex is Onset, intensification, and decline reached after 2 months (around mid-September), associated to a strong grazing pressure that overcomes the steady increase in growth rate (at r = 0 : l ′ . m ′ . 0, Figure 4b-d). Overall, the temporal changes in r are much less pronounced than for the abrupt case (Figure 4a-c) and the apex seems to arise from the long process of re-equilibrium that begins just after the onset, due to maintained top-down control. The variables P and sChl display a weak seasonal cycle with maximal values between September and October. The maximal value of sChl is reached just before P max , which in contrast to the abrupt bloom case, has a broader and lower peak. This second type of bloom are often observed in subtropical regions and low iron concentration areas of the SO.
Bloom seasonal cycle: onset, climax, and apex Onset A remarkably consistent result over our 1200 simulations is that in 100% of the situations that we have explored, the bloom onset always occurs in winter, when the mixing layer is deepening (Figure 5a). The median value of the onset date is 2 months before tXLD max and 4 months before tXLD min . Abrupt blooms tend to initiate earlier than smooth blooms, and the spread of the time of initiation is wider for smooth blooms than for abrupt blooms.
As introduced above, these winter onsets can either be growth or dilution driven. Here we assess each of the blooms in terms of the rate of change in losses and growth (l ′ and m ′ , respectively) at the time of onset ( Figure 6). In this way, we can discriminate between the bottom-up (increase in growth rate or Growth regime) and top-down (decrease in grazing or Dilution regime) control. Overall, across the 1200 scenarios of our study we find an equal distribution between Growth regimes (m ′ . 0, red quadrant in Figure 6) and Dilution regimes (m ′ , 0, blue semi-quadrant in Figure 6). No significant differences in the onset regime are found for abrupt blooms (represented by triangles in Figure 6).
In terms of absolute timing, all onsets occurring after the winter solstice are associated to Growth regimes and those occurring before the winter solstice, to the Dilution regime (see colour bar in Figure 6). Thus, in our simulations (and unlike the results of  in a model of the North Atlantic) dilution is not always efficient enough to initiate the bloom. We found that the efficiency of dilution at initiating the bloom is related to the speed of destratification of the mixing layer. Dilution is efficient  when the destratification is strong enough to dilute predators, but weak enough to retain favourable light conditions for phytoplankton. In our simulations, we find a threshold between dilution and Growth regimes of at ≈2 m d 21 (not shown). When the destratification rate is greater than ≈2 m d 21 , the decrease in losses due to dilution is not strong enough to overcome the decrease in growth and to cause the bloom to onset. In this case, the onset is delayed until light conditions become more favourable (i.e. after winter solstice), switching to a Growth regime. Alternatively, when the destratification rate is less than ≈2 m d 21 the opposite occurs, favouring Dilution regime.

Climax
After the bloom onset, biomass accumulation increases until it reaches a maximum rate that we define as the climax (time of r = r max ). We note here that this date is different from the date of maximal biomass stock: accumulation continues until r switches back to a negative value. Instead, climax refers to the maximum increase rate of integrated biomass. In contrast to the bloom onset, for the large majority of blooms (≈80%), climax is reached during the phase of XLD retreat, i.e. after XLD max (Figure 5c). The remaining 20% of blooms are characterized by a climax before or when the mixing-layer reaches its maximum depth. However, all these blooms with a climax before XLD max are smooth blooms. Abrupt blooms, associated with, by definition, a large (i.e. an "intense" climax) have their climax occurring after XLD max .
Interestingly, for all seasonal cycles analysed, we found the accumulation reaches its maximum before the surface layer re-stratifies to its minimal value (Figure 5d). In summary, we find that in 80% of our simulation climax occurs during the spring stratification (i.e. before tXLD min and after tXLD max ). In addition, climax associated with all the abrupt blooms occurs at the time where the mixing layer reaches its minimal value (i.e. at tXLD min and after tXLD max ; dark green in Figure 6c and d). These results suggest a relationship between climax date and intensity, and the surface layer re-stratification period, which in turns points out the possible importance of light on biomass accumulation rate. Therefore, we find that the faster is the restratification of the surface layer, the larger is the maximum accumulation rate (Figure 7), with the total Fe input playing a secondary role. We note that for a given restratification speed (except for the very low restratification, slower than 5 m d 21 ), r max is tightly linked to the total Fe input in the surface layer. We interpret the tight relationship between r max and the restratification as an indication of bottom-up control: a rapid improvement in light conditions leads to a parallel increase in growth rate (m), which quickly translates into a rapid elevation of r, as grazers are not able to respond at the same rate. The amount of Fe available in winter thus works as a catalyst, allowing phytoplankton to take an optimal benefit of the increase in light conditions.

Apex
The apex date is reached when loss rate first overcomes growth rate (l = m) causing biomass accumulation to arrest (r ¼ 0). For 90% of the scenarios we analysed, apex occurs after the mixing layer reaches its maximum (median value 1.5 months later; Figure 5e). In addition, apex is reached, in .75% of bloom scenarios, before the date of minimum mixing depth (median value 1 month before; Figure 5f). However, apex for abrupt blooms occurs after tXLD min (hence also after tXLD max ; Figure 5e and f). Such blooms can be viewed as examples of "bloom and bust" dynamics occurring during a rapid restratification that causes a sudden drop in r following the climax (compare Figure 5d and f; dark colours, or example of Figure 4b).
By assessing the state of the ecosystem at the apex date, we can better understand the processes leading to this stage. At apex, l ′ is always positive (Figure 8), which means grazing pressure is increasing. In contrast, m ′ can be either positive (but necessarily lower than l ′ ) or negative (Figure 8). Over the ensemble of scenarios, 68% have an increase in growth (m ′ . 0) at the time of apex, and 32% have a decrease in growth (m ′ , 0) at the time of apex. However, l ′ is always greater than m ′ (up to a factor 10 in some cases) indicating strong top-down control for our entire suite of scenarios. The earliest apex dates are always associated with m ′ . 0 (Figure 8), while after around day 260 ( mid-September), m ′ can either be positive or negative at recoupling. The highest loss rates (l ′ . 100 d −1 ) are always associated to abrupt blooms (see triangles in Figure 8).  Onset, intensification, and decline

Relating bloom phases and bloom detection methods based on surface Chl
Many of the seasonal cycles generated by our model are characterized by a plateau phase, with positive, but weak and relatively constant accumulation between onset and climax (see Section Abrupt and smooth blooms and Figure 4). While in terms of biomass accumulation, the bloom has definitely started (i.e. r), there is little accumulation of chlorophyll at the ocean surface (sChl; see examples in Figure 4, green dashed curve). Clearly, detecting bloom onset from surface chlorophyll observations would be challenging in such cases. In this section, we aim to evaluate two bloom detection methods that have been previously applied to surface observations of Chl accumulation (e.g. Behrenfeld, 2010;Brody et al., 2013), see what phase of the bloom they detect (onset or climax) and if they comply to the critical depth hypothesis or to the dilution/recoupling hypothesis.
The biomass-based bloom detection method (P*-method; see Section Methods for details) detects dates that coincide, for 85% of blooms, to the onset computed from the full vertical profile (Figure 9a). This result is not surprising given the model set-up: in our model, the mixing layer is very strongly mixed, so we expect P to be relatively constant over XLD and very weak below the mixing layer (so P * ≈ P). The accuracy of the P*-method is therefore strongly tied to the choice of the mixing depth over which P* is computed: it must be a strongly mixed surface layer which is not always well described by typical mixed-layer depth criterion (e.g. Taylor and Ferrari, 2011). Arguably, P*-method will work better during the convective phase of the surface layer, when mixed-layer actively mixes. The P*-method detects onset dates 2 months before the climax Figure 8. Growth and losses rate trends at the date of the apex for the ensemble of modelled blooms. Circles stand for smooth blooms and triangles for abrupt blooms. The absolute apex date (day of the year) is represented by coloured symbols. Figure 9. The histogram of the ensemble of modelled blooms representing: (a) bloom detection date using P*-criterion normalized to the actual onset date, (b) bloom detection date using sChl-criterion normalized to the actual onset date, (c) bloom detection date using P*-criterion normalized to the actual climax date, (d) bloom detection date using sChl-criterion normalized to the actual climax date. Median value is represented by a vertical dashed line. Abrupt blooms distribution are in dark colour and smooth blooms distribution in light colour.
( Figure 9b), in agreement with the typical time difference found between onset and climax (Figure 5a and c). On the other hand, the sChl-method detects dates that are only very weakly related to the actual onsets in our model scenarios (Figure 9c). Onset dates derived from this method are 1-4 months later than the actual onset, with a median value of 2.5 months later. Interestingly, the dates detected by sChl-method are more closely related to the bloom climax (Figure 9d), with this bias even clearer for abrupt blooms (dark green on Figure 9d).
From these results, we conclude that, usually, the P*-method is reliable on detecting the bloom onset while the sChl-method mainly detects the bloom climax.

Evaluating Sverdrup's bloom conditions from space
A significant part of the recent works that contributed to the bloom onset debate is based on ocean-colour data (Siegel, 2002;Behrenfeld, 2010;Chiswell, 2011). However, some authors have shown how bloom detection from ocean-colour data may be strongly influenced by the time-series data gaps (Cole et al., 2012) and specially by the detection method applied (Ji et al., 2010;Brody et al., 2013). With the aim to evaluate at which point the choice on the detection method influence the validation or refusal of the Sverdrup (1953)'s hypothesis, here used model data to compare the dates of bloom detection with the bloom onset date predicted by the critical depth hypothesis (Figure 10). Such a comparison was done in a similar way as for bloom phases: normalizing the dates of detection by the date at which XLD reaches the critical depth (i.e. XLD = Z c ).
This comparison shown that the P*-method detected dates 2 -4.5 months before Sverdrup (1953)'s conditions were satisfied (median value 3.6 months before; Figure 10a). Similarly, close to 95% of the dates associated with sChl-method were before XLD became shallower than Z c . However, for abrupt blooms, the dates detected by the sChl-method were well distributed around the XLD = Z c date (Figure 10b, black bars).

Discussion
During the last 20 years, many studies based on SO bloom dynamics have been conducted. Most of them rely on satellite ocean-colour observations (Moore et al., 1999;Venables et al., 2007;Fauchereau et al., 2011) except specific locations where mooring observations have been sampled (Jeandel et al., 1998;Abbott et al., 2000;Weeding and Trull, 2013), occasional oceanographic surveys (Boyd et al., 2000;Pollard et al., 2007;Blain et al., 2007), and recent datasets obtained by elephant seals equipped with CTD and fluorescence sensors (Blain et al., 2013). While in situ observations usually offer water column measurements, they are limited to specific regions and last for only a few weeks/months. In contrast, satellite-based chlorophyll data provide substantial spatial and temporal coverage. However, they are limited to interpret bloom dynamics solely based on its surface imprint.
Satellite-based analysis of high-latitude bloom onset often relate the increase of surface chlorophyll to the stratification of the mixed layer in spring. From this temporal correlation, authors conclude that alleviation of light limitation in the surface ocean layer is the main bloom trigger (Nelson and Smith, 1991;Siegel, 2002). Such results are based on the seminal concepts of Gran and Braarud (1935) and Riley (1942) and theoretically supported by Sverdrup's hypothesis (Sverdrup, 1953). More recently, combined analyses of satellite and model data identified onset based on its "strict" definition: the date at which integrated gains overcome losses (Behrenfeld, 2010). In this case, onset is systematically found in winter presumably caused by a fast decrease on grazing pressure during MLD deepening. This apparent inconsistency between the two results is subject to much debate (Chiswell, 2011;Ferrari et al., 2014).
Our results shed light on the current debate by describing the bloom as a sequence of three distinct phases: an onset, a climax, and an apex. While a "strict" onset definition is consistent with a winter onset (in agreement with Behrenfeld, 2010; the surface spring bloom is associated with the climax of the integrated bloom, which is the rapid accumulation occurring after the winter onset. One advantage of using a model approach is that it allows us to investigate the mechanisms that drive each of the bloom phases. Interestingly, two possible winter bloom triggers have been identified: grazer-prey dilution and winter net growth ( Figure 6). In addition, we find that dilution is only efficient when the destratification of the mixing layer is not too fast. When destratification is rapid, grazers are diluted, but the phytoplankton growth is reduced even Figure 10. The histogram of the ensemble of modelled blooms representing: (a) bloom detection date using P*-criterion normalized by the date at which Sverdrup's conditions are satisfied, (b) bloom detection date using sChl-criterion normalized by the date at which Sverdrup's conditions are satisfied. Median value is represented by a vertical dashed line. Abrupt blooms distribution are in dark colour and smooth blooms distribution in light colour.
Onset, intensification, and decline more strongly due to light limitation. In such cases, the winter onset is delayed to later in the season, when light conditions improve following the winter solstice. Furthermore, the climax phase is clearly bottom-up controlled and is the only phase of the bloom for which we identified a significant role of iron (enhancing the speed at which phytoplantkon accumulates). Finally, the date at which net accumulation stops (i.e. the apex) is strongly top-down controlled (Figure 8). The complete recoupling and the thereafter readjustment is influenced by many complex biogeochemical processes involving remineralization, aggregation of particles or virus infection (Boyd et al., 2012). We want to stress here that any given model, even if containing a significant level of complexity, would be suspect in the representation of these processes, which are still poorly understood. This is particularly true in the SO, where the complex cycle of iron is involved.
To complement the idealized study on the bloom dynamics and with the aim to evaluate our conclusions in satellite-based studies, we showed that it is possible to differentiate onset and climax phases using two different bloom detection methods (Figure 9a and b). In Sallée et al. (2015), these two detection methods have been used to estimate onset and climaxes dates from ocean-colour data in the SO. However, it must be noted here that (as discussed in Chiswell, 2013;Sallée et al., 2015) the method proposed to detect the onset (P * = MLD × sP, Behrenfeld, 2010) can only be successfully applied during winter destratification if three conditions are satisfied: ocean-colour satellite data are available in winter, the timing and magnitude of MLD can be accurately estimated, and the MLD is actively mixed (i.e. MLD ≈ XLD).
In the final part of this paper, we investigated at which point bloom detection methods agree with critical depth hypothesis. To do so, we computed the critical depth (Z c ) using model outputs and we compared the date at which XLD crosses Z c to the dates of bloom detection. Our results showed that Sverdrup (1953)'s blooming conditions coincided well with the dates detected by sChl-method for the case of abrupt blooms (Figure 10b), and hence with the climax phase (which was proven to be mainly top-down driven; Figure 7). Altogether, our results suggest that the top-down mechanisms identified by Gran and Braarud (1935), Riley (1942), and Sverdrup (1953) are not indicative of the bloom onset, but they are still crucial to bloom dynamics as they presumably control the climax phase. We therefore suggest that much of the debate regarding winter vs. spring onset mostly results from confusions on the definition of the word "onset". It must be emphasized to note that what originally made blooms such an attractive phenomenon was "the sudden appearance of an enormous numbers of diatoms in early spring" (Bigelow, 1926). Therefore, in our opinion, the key phase of the bloom is arguably the climax, not the onset. Indeed, it is the bloom climax (and its associated surface signal) what actually defines the observed spatial heterogeneity of SO blooms (Thomalla et al., 2011). We conclude then that the observed differences on spatial distribution of surface spring blooms between the North Atlantic and the SO regions are indicative of differences on the factors that control the climax phase; i.e. the XLD dynamics in spring and the nutriment limitation (iron for the SO). This conclusion is coherent with the fact that these two factors, and specially the coupling between them, present unique characteristics in the SO (Tagliabue et al., 2014).
Our results and conclusions are based on an idealized model where strong assumptions were applied to minimize the degrees of freedom and ease results interpretation. These simplifications and assumptions must be taken into account when interpreting the results. First of all, the seasonal cycle is modelled in a 1D water column where lateral advection is not considered. Even if this may have important consequences on nutrients/iron transport, our approach is supported by recent works on iron supply in the SO (Tagliabue et al., , 2014. Second, in our model vertical mixing is assumed to be very strong and homogeneous from the surface to a depth level (XLD) and very low below. This highly turbulent mixing layer is a reasonable assumption for the SO where winds are generally strong and sustain efficient turbulent vertical mixing. However, we did not address the sources of vertical mixing and the possible sub-seasonal variations in mixing. We note here that in the present study we analysed bloom in relation to the mixing layer depth (XLD) which is not necessarily the same as the mixed-layer depth. While Sverdrup (1953) referred to the seasonal thermocline (classically associated to the mixed layer), recent studies have focused the interest on the upper-layer mixing (based on the critical turbulence hypothesis of Huisman et al., 1999) and the mechanisms able to reduce it: positive heat fluxes (Taylor and Ferrari, 2011;Ferrari et al., 2014), wind reduction (Chiswell, 2011), or sub-mesoscale eddies (Lévy et al., 2001;Mahadevan et al., 2012). In our study, we avoided such controversy by imposing a very strong mixed upper-layer with no sub-seasonal variability. In the real ocean (and especially in the SO), the phytoplankton activity (and therefore, the bloom dynamics) is highly affected by atmospheric and oceanic physical events ranging from synoptic (Waniek, 2003) to sub-mesoscale  scales and from day to week. Such events are arguably an important source of variability when addressing the phytoplankton seasonal cycle with ocean-colour satellite data.
The idealized seasonal cycle of XLD used in this study is based on Argo observations and present three main phases during the seasonal cycle: a deepening phase (autumn -winter), a quicker shallowing phase (spring) and a transition phase (in summer) between the shallowing and the deepening. This transition phase present stable or slightly decreasing MLD (Sallée et al., 2010;Figure 1c). However, the evolution of XLD throughout the seasonal cycle is actually more complex in the SO (e.g. Swart et al., 2014), with short (i.e. sub-seasonal) and rapid (days to weeks) deepening/shallowing events. Such events are likely to influence the integrated accumulation of phytoplankton and the dates at which the onset, climax, and apex are reached. Among the three phases, the climax is by far the most sensitive to rapid changes on stratification (Figure 7). On the other hand, our results suggest that sub-seasonal XLD variability may weakly affect the onset and apex dates. The reason is that onset and apex controls (grazers dilution or low net growth for the former, and grazing pressure for the latter) are mainly related to the phytozooplankton coupling which is much less sensitive to rapid ( day) changes in mixing depth.
Finally, it must be emphasized to point out that in our set of experiments the limiting nutrient was always dissolved iron while is known that in Fe-rich water of the SO, diatoms can also be limited by silicic acid (Boyd and Ellwood, 2010). We therefore note that our results should not be extrapolated to any other location or SO regions where main iron supply is not winter mixing (in the lee of Island and shallow plateaus), where silicic acid is a limiting factor, nor to higher latitudes (.708S) where the role of light and ice seasonal cycle can be critical on the phytoplankton bloom phenology.

Conclusions
Implementing an ensemble of .1200 idealized physical scenarios for an isolated water column, a complex biogeochemical model has been forced with the aim to reproduce the plankton seasonal cycle for a collection of open waters/ice-free SO spots. Daily frequency model outputs covering a large spectra of the variables involved in the phytoplanktonic bloom allowed us to target the question of bloom formation mechanisms from different focus.
Three crucial stages of bloom seasonal evolution have been defined: onset, climax, and apex date. All onsets occurred in winter and a large majority ( 80%) of the climax, in spring. For the onset, upper-layer mixing (or XLD) appeared as a key component on tilting the system to be bottom-up or top-down controlled. For climax, the amount of Fe (and thus the relative depth between mixing layer and ferricline) seemed to play a secondary but significant role on the intensity of accumulation. Concerning apex, permanent top-down control was identified.
Two bloom detection criteria were tested using model surface chlorophyll and mixed-layer integrated biomass estimated from surface values. The biomass-based method appeared as a good proxy for detecting bloom onset while the method based on surface chlorophyll was reliable on detecting the climax. Finally, we compared the date at which XLD crosses the critical depth with the dates of bloom detection by the biomass-and the sChlbased methods. Sverdrup (1953)'s blooming conditions fairly coincided with the dates of detection using the sChl-method, and therefore with the bloom climax.
Our results suggest the existence of bottom-up as well as topdown drivers of the different phases of the blooms. It also enlightens the apparent controversy between onset/surface bloom detection and it shows that how different criteria can be used to answer different questions.