Synoptic-scale analysis of mechanisms driving surface chlorophyll dynamics in the North Atlantic

. Several hypotheses have been proposed for the onset of the spring phytoplankton bloom in the North Atlantic. Our main objective is to examine which bottom-up processes can best predict the annual increase in surface phytoplankton concentration in the North Atlantic by applying novel phenology algorithms to ocean colour data. We construct indicator ﬁelds and time series which, in various combinations, provide models consistent with the principle dynamics previously proposed. Using a multimodel inference approach, we investigate the evidence supporting these models and how it varies in space. We show that, in terms of bottom-up processes alone, there is a dominant physical mechanism, namely mixed-layer shoaling, that best predicts the interannual variation in the initial increase in surface chlorophyll across large sectors of the North Atlantic. We further show that different regions are governed by different physical phenomena and that wind-driven mixing is a common component, with either


Introduction
About half of global primary production is performed by marine phytoplankton. Phytoplankton production fuels marine ecosystems and the harvesting of marine living resources, as well as playing an important role in global carbon cycling (Field et al., 1998). In many parts of the world's oceans, marine primary production undergoes a distinct seasonal cycle, with the major part of production occurring in the spring bloom (Longhurst, 1995;Martinez et al., 2011;Platt et al., 2010). This seasonal cycle is particularly apparent in the North Atlantic (Yoder et al., 1993), where it imprints seasonal variations in species abundance and annual routines (e.g. spawning, migration) throughout the marine food web from zooplankton (Gaard, 2000;Gislason and Silva, 2012;Heath et al., 2000) to fish (Badcock and Merrett, 1976;Trenkel et al., 2014) and marine mammals (Pauly et al., 1998). In the North Atlantic, the progression of primary production throughout the year, and its variation between years, is commonly used as a proxy for ecosystem state (Frajka-Williams and Rhines, 2010;Lévy et al., 2005;Townsend et al., 1994). The North Atlantic spring bloom is an important biological event and has attracted considerable scientific attention during the last decades (Behrenfeld, 2010;Chiswell et al., 2013;Platt et al., 2003).
Phenology is the term used to describe the study of the timing of annually recurring biological events, such as the observed "greening" of the surface ocean, an indicator of bloom initiation. Phenology provides a standard way of understanding the cascading fluctuations throughout the food web. To achieve this, a good phenology metric should be accurate, precise, and sensitive to the underlying environmental processes, both physical or biological (Ferreira et al., 2014). Much of the recent interest in spring bloom dynamics (Behrenfeld, 2010;Chiswell et al., 2013) concerns the mechanisms that influence different characteristics of the annual cycle.
Chlorophyll concentration is, arguably, the most important ecological variable setting the pace of life in temperate and high-latitude seas. In this study, we use surface chlorophyll concentrations as derived from satellite ocean colour to detect spring bloom initiation (Cole et al., 2012;Sasaoka et al., 2011;Brody et al., 2013). We thus assume that the chlorophyll concentration at the surface represents that of the surface mixed layer (Evans and Parslow, 1985). While we note that some aspects of bloom dynamics are more properly described by integrating phytoplankton biomass over the mixed layer (Behrenfeld, 2010), it is the surface chlorophyll that is the most readily accessible via the highly resolved (both spatially and temporally) ocean colour products.
There are essentially three environmental processes that can change the surface chlorophyll concentration: phytoplankton growth through light and nutrients; loss terms, such as respiration, grazing, coagulation, and sinking; and dilution through mixed-layer deepening. These processes are particularly important during two key phases of the seasonal cycle: (1) events that lead to an increase in phytoplankton biomass -bloom initiation -and (2) conditions that halt the net increase in biomass -the peak of the bloom. Phytoplankton biomass will increase whenever the growth rate exceeds the loss rate (Sverdrup, 1953). This picture, as regards the distinction between biomass and surface chlorophyll concentration, is somewhat complicated by dilution; a deepening mixed layer dilutes the concentration but has no effect on the biomass, a process that has repercussions on the feeding success and thus population dynamics of grazers. However, a shoaling mixed layer has no direct influence on the concentration but removes biomass to some extent. These processes and their implications for phytoplankton, the resources they rely on, and their grazers have been carefully considered in recent reanalyses of spring bloom dynamics (Behrenfeld et al., 2013a;Lindemann and St John, 2014).
It is also fair to say that the annual trajectory of phytoplankton biomass and surface phytoplankton concentration follow different dynamics (Chiswell et al., 2013). While we recognise that phytoplankton biomass variation is an important aspect of spring bloom dynamics, in this paper, we examine which fundamental physical processes may best predict the timing of the increase in surface phytoplankton concentrations. Furthermore, we do so since ocean surface colour is a readily available synoptic-scale observable spanning many years of measurements. The interannual variability in bloom timing is evaluated in terms of how much the increase in surface layer chlorophyll is advanced or delayed compared to the day of the climatological maximum rate of increase.

Mixed-layer shoaling
Over the years, several theories have been put forward which, in one way or the other, try to model the growth and loss rates in terms of fundamental processes (Table 1 and Fig. 1). The classic application of the growth-loss view of bloom initiation relates to when photosynthetic production of organic matter surpasses respiration (Sverdrup, 1953), where respiration refers to all losses and is constant. This hypothesis is commonly referred to as the "critical-depth hypothesis", which states that a bloom begins when the surface mixed layer shoals to a depth above the critical depth (where integrated production equals losses). The shoaling of the mixed layer means that individual phytoplankton cells remain in the euphotic zone for longer (Chiswell, 2011;Siegel et al., 2002;Sverdrup, 1953). By extension, this suggests that the light intensity integrated over the mixed layer is the most relevant factor driving phytoplankton blooms in the North Atlantic. Here, we term this hypothesis the "critical-depth model" (Table 1).

Active mixing
Mixed-layer shoaling, however, is not the only process which can increase the residence time of primary producers in the well-lit surface ocean. Similar effects can be driven by periods of low surface mixing (Townsend et al., 1992). This has led to a series of alternative interpretations, which highlight active mixing (specifically the lack thereof) as a key factor (Huisman et al., 1999;Taylor and Ferrari, 2011a;Townsend et al., 1994).
One of the first quantitative studies (Townsend et al., 1994) examined the combined effects of wind-driven mixing and light, the hypothesis being that blooms can occur during periods when light is low but increasing and turbulent mixing weakens. These conditions can be met well before the surface mixed layer begins to shoal. We call this the "critical-lightexposure model" (Table 1).
This type of reasoning can also lead to only the competing effects of stratification by solar heating and of destratification by wind-driven mixing being considered. This view encapsulates the key elements of the "critical-turbulence model" (Huisman et al., 1999, where brief interludes in mix-  Huisman et al. (1999;  Critical light exposure L: light intensity. M: wind-driven mixing R ∼ α 3a L + α 3b M + β 3 Townsend et al. (1994) Critical heat flux Q: heat flux. M: wind-driven mixing R ∼ α 3a Q 0 + α 3b M + β 3 Taylor and Ferrari (2011a, b) ing and heating produce a stable layer in which phytoplankton cells are retained within the euphotic layer. Thus, a balance between heat flux and wind-driven mixing may explain North Atlantic phytoplankton seasonality (Table 1). More recently, Taylor and Ferrari (2011b) have shown that blooms may be detected much earlier than the shoaling of the mixed-layer depth, and it has been proposed that blooms can be initiated as soon as deep convection ceases (Taylor and Ferrari, 2011a), that is, as soon as the ocean experiences a net inward heat flux. In this context, the timing of the transition from net cooling to net warming is a key element linked to the variability in phytoplankton seasonality. We term this the "critical-heat-flux model" (Table 1).

Other processes not considered
There have been theories also focusing on specific regional effects. For instance, Mahadevan et al. (2012) were able to link bloom onset to eddy-driven stratification, prior to net warming. Fronts were also found to trigger high-latitude blooms by reduced mixing, which explains high chlorophyll levels in light-limited regions (Taylor and Ferrari, 2011b). Other studies (Frajka-Williams and Rhines, 2010;Garçon et al., 2001;McGillicuddy et al., 2007) have also linked spring bloom initiation to offshore advection, eddy-induced upwelling, or river runoff. Finally, oceanic convection has been found responsible for significant vertical transport, thus maintaining a winter stock of phytoplankton in the deep mixed layer that can potentially reseed the spring bloom (Backhaus et al., 1999(Backhaus et al., , 2003. Behrenfeld (2010) adopted a different approach by examining the influence of dynamic top-down controls, suggesting the "dilution-recoupling hypothesis". This is a concept that is implicit in the model of Evans and Parslow (1985). The hypothesis of Behrenfeld (2010) proposes that a vertically integrated biomass increases in midwinter with the increase of day length, even when the mixed-layer depth is at its deepest, and reaches its maximum with the recoupling of grazers due to stratification. Unfortunately, as also noted by Behrenfeld (2010), data on top-down controls remain elusive at the spatial and temporal resolutions necessary to test this hypothesis against the complex structure of North Atlantic phytoplankton seasonality.

When and why does a surface bloom start?
As noted by Cole et al. (2015), assessing the drivers of bloom initiation variability may lead to an understanding of what starts the bloom in the first place. Despite all of the abovementioned hypotheses, there is still no clear consensus regarding a single main driver of North Atlantic spring blooms. Additionally, the spatial application of these theories may not hold true in smaller regions, where local forcing plays a more important role. Nonetheless, the key process, common to all hypotheses of surface bloom initiation, is based on the spring stabilisation of the water column, when both light and nutrients are at sufficient levels, whether by mixed-layer shoaling (Sverdrup, 1953) or by weakening turbulent mixing (Huisman et al., 1999Taylor and Ferrari, 2011a, b;Townsend et al., 1994). Their main differences reside in the physical proxy for bloom initiation: what physical indicator best predicts bloom timing?
While there are a number of metrics that can be used to delineate bloom initiation (Yoder and Kennelly, 2003;Rolinski et al., 2007;Siegel et al., 2002), our goal is to seek a metric that can be credibly related to the processes proposed above, i.e. those that relate to the preconditioning of the water column prior to surface bloom initiation. In this regard, any metric that uses the bloom peak (such as the popular 5 % above annual median; Siegel et al., 2002), or seasonally integrated chlorophyll, will be handicapped because it inherently takes into account not only what starts the bloom but also what terminates it some weeks or months later. We seek instead a phenology metric that is not confounded by the bloom peak, does not require winter values, and is a straightforward indicator of the greening of the surface ocean as observed from space. Our metric is based on how advanced or delayed the development of surface chlorophyll concentration is in a particular year compared to the climatological date and rate of maximum concentration increase.
We construct four models based on the literature using a range of physical observations, primarily from satellite but also from model data, and describe key processes observed www.biogeosciences.net/12/3641/2015/ Biogeosciences, 12, 3641-3653, 2015 in the North Atlantic (Table 1 and Fig. 1). In each case, we make the models as simple as possible -capturing the essential process dynamics in terms of, at most, two observable or estimated fields only. We use the information theoretic (IT) approach to investigate which model for surface blooms has the most support within the North Atlantic. The IT approach is a very useful tool when comparing different models. In particular, it provides a rigorous framework for evaluating the evidence in support of competing models. It does so by defining a priori a set of "multiple working hypotheses" rather than a single alternative to the null hypothesis. The IT approach is then followed by expressing each hypothesis in quantitative terms that represent a hypothesis's strength of evidence to be used further in the model selection (Burnham et al., 2011).
We conduct our study focusing on bottom-up controls that may trigger a North Atlantic phytoplankton surface bloom, and we thus neglect the effect of top-down controls (grazing, Behrenfeld, 2010;Evans and Parslow, 1985;Irigoien et al., 2005). Information on top-down controls is not available at the spatial and temporal coverage needed to assess mesoscale physical forcing. In addition, as Chiswell (2011) shows, the seasonal cycle of surface chlorophyll differs from the vertically integrated chlorophyll. The dilution-recoupling hypothesis of Behrenfeld (2010) applied to vertically integrated chlorophyll blooms, while the other hypotheses (Huisman et al., 1999Platt et al., 1991;Siegel et al., 2002;Sverdrup, 1953;Taylor and Ferrari, 2011a, b;Townsend et al., 1994) can be applied to surface chlorophyll. Our aim is to compare the latter, in which it is assumed that surface blooms only take off when the surface waters stabilise.

Information theoretic (IT) approach
The main aspects of the IT framework (Akaike, 1973;Burnham et al., 2011;Burnham and Anderson, 2002) in the context of our study include (1) identifying plausible mechanistic hypotheses and (2) a strong reliance on the quantitative evidence of factor(s) affecting a response variable, rather than a formal assessment of the statistical significance of such a factor or factor(s). In our study, (1) is expressed through mathematical descriptions of the different hypotheses to be tested (see Table 1 and Sect. 2.2), while (2) is covered by ranking the spatial evidence of the models using the concept of model selection and multimodel inference (see Burnham et al. (2011) and Sect. 2.5).

Physical mechanisms
We are particularly interested in knowing how much information from raw data is correlated to surface chlorophyll. Raw data refers to the original data in their simplest form, without preprocessing. Therefore, we quantitatively translate the fundamental physical processes that can be used to predict a phytoplankton surface bloom in the North Atlantic into simple and straightforward models (Table 1 and Fig. 1). Critical depth -A bloom initiates if the mixed-layer depth (MLD, H ) shoals below the critical depth, so light (photosynthetically active radiation, PAR, L) becomes available to phytoplankton cells (Fig. 1a) (Sverdrup, 1953;Siegel et al., 2002;Platt et al., 1991). Therefore, L integrated over H provides an estimate of the light available at the euphotic depth for phytoplankton to grow.
Critical turbulence -A bloom initiates if there is a balance between buoyancy (heat flux, Q) and wind-driven mixing (M; Fig. 1b) (Huisman et al., 1999. Critical light exposure -A bloom initiates if winddriven mixing (M) is at a level low enough to allow cells to experience surface light conditions (L; Fig. 1c) (Townsend et al., 1994).
Critical heat flux -Bloom initiation is associated with the date when net warming starts (Q ≥ 0), and low wind-driven mixing (M) increases the residence time of phytoplankton in the euphotic layer ( Fig. 1d) (Taylor and Ferrari, 2011a, b).

Data sets
In order to gather the information necessary to formulate the models for the North Atlantic domain, we used satellite observations (chlorophyll concentration, attenuation coefficient and photosynthetically active radiation), model estimations for the variables where satellite data were not available (mixed-layer depth), and model and observational merged data (wind stress and heat flux).
We used products derived from the European Node for Global Ocean Colour (GlobColour Project, http://www. globcolour.info/). The GlobColour Project blends observational data from the Sea-viewing Wide Field-of-view Sensor (SeaWiFS), the Moderate Resolution Imaging Spectroradiometer (MODIS-AQUA), and the Medium Resolution Imaging Spectrometer (MERIS) instruments by using the Garver-Siegel-Maritorena (GSM) algorithm (Maritorena et al., 2002) to generate a merged, global ocean colour product. Combining the three sensors increases the data coverage in both time and space, thus providing significantly elevated spatio-temporal coverage (Maritorena et al., 2010) and making it a common choice for phenology studies (Cole et al., 2012;Kahru et al., 2011). For this study, we chose to use daily, 0.25 • -resolution level-3 mean chlorophyll concentration (C) and attenuation coefficient (K d ) products (based on the analysis performed by Ferreira et al., 2014), from 1998 to 2010 inclusive, thus providing a total of 13 years of data. The surface photosynthetically active radiation (PAR, L) was obtained from the SeaWiFs data centre (http:// oceancolor.gsfc.nasa.gov/). We used a daily, 9 km resolution product from 1998 to 2010. These data were further gridded onto 0.25 • grid using linear interpolation to match the spatial resolution of the other data sets.

Biogeosciences
The mixed-layer PAR (L H ) was defined as L integrated from the surface to the depth of the mixed layer H : using the relevant K d reported by Irwin et al. (2012) and Cole et al. (2015). Mixed-layer depth (MLD, H ) data were obtained from TOPAZ 4 reanalysis (Towards an Operational Prediction system for the North Atlantic European coastal Zones; (Sakov et al., 2012). The TOPAZ system is a coupled ocean-seaice data simulation system for the North Atlantic and Arctic Ocean with a resolution of 12-16 km, and it is the main forecasting system for the Arctic Ocean in Copernicus (http://www.myocean.eu) and the Norwegian contribution to the Global Ocean Data Assimilation Experiment (GO-DAE). It uses the Hybrid Coordinate Ocean Model (HY-COM, http://hycom.org/hycom/) (Bleck, 2002). HYCOM is coupled to an elastic-viscous-plastic (EVP) sea ice model (Hunke and Dukowicz, 1997) and a thermodynamic module (Drange et al., 1996). The model assimilates sea surface temperature, altimetry, ice concentration, ice drift, and available in situ measurements with the ensemble Kalman filter (Evensen, 2003). The model daily output is binned onto a 0.25 • regular grid. The MLD is calculated using a density criterion with a threshold of 0.01 kg m −3 (Petrenko et al., 2013) from 1998 to 2010.
Wind stress (τ wind ) is used as a measure for wind-driven mixing (M) (Simpson et al., 1981;Taboada and Anadón, 2014) and was estimated by using M ∝ |τ wind | 3/2 , which is proportional to the power exerted by the wind on the surface ocean and the turbulent kinetic energy used in the calculations of Brody and Lozier (2014) calculations of the mixing length scale.
Both τ wind and heat flux (Q) data were gathered at a spatial resolution of 1.875 • × 1.905 • from the National Centers for Environmental Prediction (NCAR) and the National Center for Atmospheric Research (NCEP) (Kalnay et al., 1996). These data sets were further gridded onto 0.25 • grid using linear interpolation to match the spatial resolution of the other data sets.
All data sets started on 1 October 1997. We only focused on latitudes north of 40 • N due to the fact that lower latitudes have a less well-defined seasonal cycle (Follows and Dutkiewicz, 2011;Brody and Lozier, 2014).

Metrics
One of the fundamental aspects of spring bloom is the rapid increase in surface chlorophyll concentration; a phenomenon that can be interpreted as bloom initiation. In this work, we choose a bloom initiation metric that relates to how advanced or delayed the surface chlorophyll concentration is in a particular year, compared to the climatological date of maximum surface concentration increase. We term this the rate of change phenology anomaly (RPA, R). This metric has the advantage of not depending on the maximum chlorophyll concentration (an indicator of the peak of the bloom). Neither does it depend on winter values, which are usually missing from remote-sensing products (Ferreira et al., 2014), or on vertical integration (Behrenfeld, 2010), both of which introduce extraneous factors into the mechanistic reasoning as to the onset of a bloom. These are all limitations that occur in many other metrics used in the literature (Brody and Lozier, 2014;Sharples et al., 2006;Siegel et al., 2002). We decided to use an anomaly of surface chlorophyll because it is a more relevant measure with regard to higher trophic levels and is, we believe, closer to bloom preconditioning. Additionally, in order to use an integrated chlorophyll field, we would need to use modelled mixed-layer depth, which is incompatible with testing one of our key models.
At each location x, y (each 0.25 • pixel), we estimate the climatological pattern of surface chlorophyll concentration C(x, y, t) by applying a generalized additive model (GAM) to the observations from 1998 to 2010 (Fig. 2). We then calculate the day of the year when the climatological mean exhibits the maximum rate of increase, g(x, y) = max{dC/dt}. We define the climatological date of maximum increase as T 0 = t : dC/dt = g, and we define the climatological chlorophyll concentration on that day asC 0 =C(x, y, T 0 ). For each year i and location x, y, we fit a GAM with a smooth spline to the period T 0 ± 15 days for observed surface chlorophyll to produce C i (x, y, t). Lastly we define the rate of change phenology anomaly as R i (x, y) = C i (x,y,T 0 )−C(x,y,T 0 ) g(x,y) . Thus, the RPA metric R i (x, y) is a value in days and relates to how advanced or delayed the seasonal development of chlorophyll concentration is in each year i compared to the climatology of the bloom. We set a threshold that at least three observations must exist within the 30-day window for the RPA method to be valid. We apply spatial kriging with a maximum radius of 250 km to fill in pixels where the method cannot be used, e.g. due to missing data around T 0 in some years or low seasonality.
We investigate the spatially dependent ranking of the models (Table 1 and Fig. 1) using the IT approach. Thus, we construct indicator fields and time series which, in various combinations, provide models consistent with the principle physical dynamics observed in the North Atlantic. At each location, we apply a centred moving average of 30 days to the physical driver observations, and these averages will be www.biogeosciences.net/12/3641/2015/ Biogeosciences, 12, 3641-3653, 2015 Calculation of the rate of change phenology anomaly for each location x, y, i.e. each 0.25 • pixel (R i (x, y)). (a) Each seasonal cycle (solid black and blue lines) is used to estimate the climatology (C(x, y, t), dark red line). (b) The maximum increase g inC and the day on which it occurs (T 0 ) are used as a reference to estimate how delayed or advanced surface bloom is each year.
(c) A 30-day window around the T 0 is isolated for each year's seasonal cycle. R i (x, y) is estimated from the difference between annual C i (T 0 ), climatologyC(T 0 ), and g. The R i (x, y) is thus a value in days.
referred to as L , L H , M , and Q . We also use Q 0 for the date when Q becomes positive (start of net warming) and remains positive for 7 consecutive days. We further applied an inverse distance-weighted interpolation (using the weighted average of the values at the known pixels) to all thresholds to fill in the pixels where the thresholds could not be estimated. All pixels in waters shallower than 200 m were removed as coastal regions have higher associated biases (Maritorena et al., 2010) due to high turbidity and consequent different optical properties (McCain et al., 2006;Longhurst et al., 1995;Sathyendranath et al., 2001;Antoine et al., 1996).

Analysis
There are several model selection tools that can be used for comparing and ranking models. In our IT approach, we used the Akaike information criterion (AIC) (Burnham et al., 2011), which is based on the residual sum of squares (RSS) from each model. By comparing and ranking the evidence from different models, their relative importance can be quantified. Since we only aimed at assessing 13 years of data (from 1998 to 2010), we used the AICc. The AICc is the AIC corrected for small samples. Theoretically, as sample size increases, AICc converges to AIC. Another model selection unit is the Akaike weight, which can be either based on the AIC or the AICc. The Akaike weight is a value between 0 and 1 representing the weighted mean probability of  Table 1. Based on the weight of each model, we were also able to select the most strongly supported model for each 0.25 • pixel.

Results
From the four hypotheses considered (critical depth, critical turbulence, critical light exposure, and critical heat flux) within each 0.25 • pixel, the one with the highest Akaike weight is selected as the winning hypothesis (Fig. 3), where we see that the critical depth seems to be the most frequent winning hypothesis.
The spatial distribution of winning hypotheses shows no systematic pattern with regard to basin, depth, or latitude (Fig. 3). We also ran this analysis with two other bloom timing metrics: 5 % above annual median (Henson et al., 2010;Racault et al., 2012;Siegel et al., 2002) and maximum increase in chlorophyll concentration (Rolinski et al., 2007;Sharples et al., 2006;Wiltshire et al., 2008). We found similar results: no systematic pattern (results not shown).
In spite of the general dominance of the critical-depth hypothesis, there are, however, regions that show some coherency: the critical turbulence appears to be well supported mainly off Newfoundland; the critical heat flux has local sup-  port north of Iceland and in the Labrador Sea; the critical light exposure appears to have a wider distribution with very low frequencies. The spatial distribution of Akaike weights (Fig. A1) indicates the strength of support for the "winning" hypothesis. There are regions where the weights are close to 1, indicating that the corresponding models are clear winners. Some of these regions are the same as the ones observed in Fig. 3, for instance, offshore of Newfoundland, suggesting strong support for the critical-turbulence hypothesis in this region. A pixel-wise multimodel inference approach also allows the quantification of the number of occurrences of each of the four alternative hypotheses as the winning one (Fig. 3). There are no clear differences in the ranking units of the three less frequent hypothesis (0.15, 0.11, and 0.07), whilst the critical depth showed a higher-ranking unit (0.67).

Biogeosciences
To better understand the effect of each physical component (L H , L , M , Q , Q 0 ) within the four hypotheses (Fig. 1), we built single-variable models (linear regressions) using each component as the variable for each location (Fig. 4). The most frequent winning physical driver based on the Akaike weights is heat flux Q . Its spatial distribution dominates off Newfoundland, in the subpolar gyre and intermediate gyre regions, and in the Bay of Biscay. Its dominance is, however, only slightly greater than the other physical components.

Discussion
The phenology of spring bloom characteristics (e.g. initiation, peak) is thought to be controlled by a number of mechanisms including bottom-up and top-down processes. Here we specifically set out to test various bottom-up processes that can be used as indicators of phytoplankton surface blooms, testing several simplified hypotheses across a broad extent of the North Atlantic. In this regard, spring surface bloom initiation is problematic in that defining it has as much to do with what limits the bloom amplitude as what starts it in the first place. Moreover, limiting factor(s) can be the ultimate switching mechanism needed for a bloom to start. Instead, we seek to explain what bottom-up processes determine the interannual variability in bloom development around the time when, climatologically, one would expect the maximum rate of increase in surface chlorophyll concentration. By quantifying each physical mechanism independently, we observe that, even though there is no clear losing mechanism in the North Atlantic domain, the classical theory (critical depth) of Sverdrup (1953) still dominates, i.e. it provides superior evidence supporting the interannual variability in timing across the greatest range of space in the North Atlantic (Fig. 3).
All of the four alternative hypotheses are expressed as simple interpretations of what potentially drives the surface blooms in the North Atlantic on the mesoscale (Fig. 1). The models are constructed so as to be as simple as possible, using at most two out of four physical observables (light intensity, light intensity integrated over the mixed-layer depth, wind-driven mixing, and heat flux) in various combinations. Each model is based on one of the two classes of mechanisms discussed in the introduction: mixed-layer shoaling (critical depth) or active mixing (critical turbulence, critical light exposure, and critical heat flux). Our study shows the strength of the critical-depth model and indicates a dominance of the mixed-layer shoaling over the active mixing mechanism, but not everywhere.
There is an apparent inconsistency between our results and some recently reported results, notably those of Cole et al. (2015) and Brody and Lozier (2014). In the former, the strongest relationship with bloom initiation was found for the date of zero heat flux (Q 0 ), while in the latter it was for the shoaling of mixing length (essentially, heat flux tempered by wind stress and stratification). There are, however, several reasons why the results may differ. Firstly, Brody and Lozier (2014) tested the climatological bloom initiation date against the various drivers in a spatial context rather than testing the interannual variations in a temporal context as we do here. In contrast, Cole et al. (2015), while maintaining the temporal aspect, reduced each seasonal cycle of potential drivers to a single annual metric, e.g. the date when the mixed-layer depth shoals most rapidly. Precisely how these different aggregation processes influence the outcome of statistical treatments remains unresolved. More importantly, the bloom initiation metrics chosen by each of these studies are also different. Cole et al. (2015) chose the 5 % above annual median as their metric , a metric that may be less than reliable with regard to bloom initiation. Brody and Lozier (2014) used the date of the first increase of surface chlorophyll concentration (F 0 ), given by F 0 = t : dC/dt = 0 rather than our date of maximum increase T 0 = t : dC/dt = g. While it may be debated which of these have greater significance (and for which ecosystem process), it also underscores an important issue: that different milestones in the seasonal development of the spring bloom may well come under the influence of different dynamics. In our study, even though the critical-depth hypothesis is the winner (most spatially frequent), the spatial distribution of the winning model shows regions where the mixed-layer shoaling mechanism seems not to be supported. For instance, there is a dominance of the critical-turbulence and criticallight-exposure models in the Bay of Biscay. This may be due to the high degree of upwelling in this region; hence the failure of the critical-depth hypothesis to predict surface bloom dynamics. Another example occurs east of Newfoundland, where the critical-turbulence and critical-heat-flux hypotheses dominate. Both of these hypotheses have wind-driven mixing as a common parameter. In addition, heat flux and light intensity are also key individual drivers in this region (as confirmed by the map in Fig. 4). These findings suggest that spring bloom seasonality in these regions may be driven by periods of reduced active turbulent mixing, increasing exposure to light (Huisman et al., 1999;Taboada and Anadón, 2014). The region off Newfoundland is also very energetic (high physical forcing), highly influenced by the subpolar gyre, and serves as a path for the northward movement of the Gulf Stream waters. The failure of critical depth to explain the bloom dynamics in this region may be due to the subduction of cold waters from the subpolar gyre and the warm waters from the North Atlantic drift. This may explain why the critical turbulence and the critical heat flux were dominant in the region east of Newfoundland and into the central North Atlantic. These 3-D processes should be tested in the future to help understand the dynamics of the North Atlantic system.
The explanatory power of the hypotheses that assume the mechanism of active mixing (critical turbulence, critical light exposure, and critical heat flux) is fairly evenly distributed (Fig. 3). These three hypotheses seem to operate with a switch-on mechanism, i.e. a number of conditions have to be met for bloom growth, and any one may be the critical condition that triggers the growth spurt. This interpretation is supported by comparing Figs. 3 and 4, where the criticaldepth model is a clear winner in the model intercomparison but only scores averagely when tested against individual parameters. In this case, the limiting conditions appear to be either light intensity or heat flux (since all three have winddriven mixing as a common parameter). Our results show that there is no clear winning hypothesis among these three ac-tive mixing models, but there is a bias towards mechanisms involving heat flux (Fig. 4). This finding is supported by Taylor and Ferrari (2011a), where a bloom develops due to the start of net warming, weakening turbulent mixing, and a subsequent increase in the residence time of phytoplankton cells within the euphotic layer. In order for this to happen, a standing stock of phytoplankton cells needs to exist a priori. The "seed stock" consists of the remnants from the previous year that survived all winter at depth due to convection. As suggested by Backhaus et al. (2003Backhaus et al. ( , 1999 and Chiswell (2011), deep convection spreads out the overwintering remnants, but, as soon as stratification begins, those at the surface start to bloom. From our results (Fig. 3), we confirm that heat flux is a strong physical driver. Thus, in regions where the critical depth is not the winning model, the active mixing mechanism (either triggered by light intensity or heat flux) seems to play an important role.
The second most common physical property was winddriven mixing (Fig. 4), and this is the common parameter in the models concerning the active mixing mechanism. In the past, the importance of wind-driven mixing has been shown by Huisman et al. (1999), ,  and confirmed by Taylor and Ferrari (2011a, b). The first group of authors stresses a balance between wind-driven mixing and sinking rates, so that an intermediate mixing allows both enough surface nutrient replenishment and sufficient average light exposure. Recently, Taboada and Anadón (2014) suggested that wind forcing (wind stress as a proxy for wind surface mixing) played a key role in bloom timing and magnitude (see their Fig. 5a and c). The results shown by these authors are based on single-parameter hypotheses (not including heat flux) and confirm that spring blooms are triggered by different physical properties in different mesoscale regions. Our results are thus in agreement where wind stress is found as a common parameter within the North Atlantic domain.
Winds have essentially two effects: turbulent mixing (Backhaus et al., 2003;Townsend et al., 1994) which is only shallow (around 50 m in midlatitudes) and surface cooling which promotes deep convection (Backhaus et al., 2003;Brody and Lozier, 2014). Together with the cessation of convective overturn, wind stress decreases during the spring. Deep mixing is therefore no longer active, and there is a shift from a deep mixed regime to a shallow light-driven regime. However, it is important to note that the depth of the mixed layer is not the same as the depth of the vertical mixing of plankton (Chiswell, 2011). These two depths only match when vertical mixing is at its limit (Taylor and Ferrari, 2011a). In the presence of low vertical mixing, a surface bloom can initiate even if critical-depth conditions (Sverdrup, 1953) are not met, i.e. even if the thermocline is deeper than the critical depth. This mechanism is presented by Chiswell (2011) as the "stratification-onset model", regarding which the author contends that the critical-depth hypothesis is valid during autumn and winter, when the deep-Biogeosciences, 12, 3641-3653, 2015 www.biogeosciences.net/12/3641/2015/ ening thermocline may suppress production due to the downward mixing of plankton, but not in spring, since the upper layers are not well mixed in plankton. The model is consistent with the findings by Taylor and Ferrari (2011a), in which surface stratification results from the cessation of convective overturn and low wind stress. In our study, we show that the critical-depth hypothesis is still valuable for predicting phytoplankton spring surface blooms in the North Atlantic. Our findings, however, are underlain by assumptions that are worth considering. Firstly, we based the critical-depth hypothesis on Sverdrup's classical theory, thus only accounting for L H . This makes the model inherently simpler. The other three hypotheses use two parameters separately and are therefore somewhat handicapped (higher penalty due to higher number of parameters) when compared to the critical depth. We believe that this type of study would improve if similar combinations were found for the remaining hypotheses: critical turbulence, critical light exposure, and critical heat flux. For this reason, we tried to use a two-parameter approach (considering H and L separately) for the criticaldepth hypothesis, so that the four models would have the same number of parameters and the AICc weights would thus be comparable. The critical depth explained by L H alone showed itself to be inherently superior (with a much stronger signal) than the combined H and L model; thus, we chose to keep our interpretation of the critical-depth hypothesis using L H . This underscores the point that physical reasoning can significantly improve in improving model predictions.
Secondly, we recognise that our study assumes that the same mechanism predicts surface bloom timing at a given location for the entire time frame (from 1998 to 2010). However, it is conceivable that different mechanisms may be best predictors in different years. Considering the high variability in the spatial distribution of the models (Fig. 3), it is reasonable to expect a similar high temporal variability. In the same way we observe that different mechanisms dominate in different regions; intuitively, one can assume that different mechanisms will also dominate in different years. Indeed, given the scatter in winning models, it is entirely conceivable that bloom timing is governed by a limiting factor and that multiple conditions have to be met, any one of which may be the trigger in any given year or location.
Thirdly, we also recognise that our study fails to assess top-down mechanisms. A key hypothesis that has been attempted by Brody and Lozier (2014) is the dilutionrecoupling hypothesis (Behrenfeld, 2010). Brody and Lozier (2014) found very little correspondence between seasonal thermocline increases and integrated chlorophyll increases. However, as they noted, in order to successfully study this hypothesis, one would require temporally and spatially distributed data on grazing pressure and encounter rates between grazers and phytoplankton. Since such highly resolved data sets are not available, top-down mechanisms cannot be properly assessed at this time.

Conclusions
The complexity of spring bloom dynamics in the North Atlantic has been discussed since Sverdrup (1953) published the "critical-depth hypothesis". The discussion took a different direction when Behrenfeld (2010) suggested a topdown control of the phytoplankton seasonal cycle with the "dilution-recoupling hypothesis". Various studies followed the same line of thought (Behrenfeld et al., 2013a, b, c;Irigoien et al., 2005). However, bottom-up factors are still the most studied (Huisman et al., 1999;Siegel et al., 2002;Taylor and Ferrari, 2011a;Townsend et al., 1994), especially because data are more readily available than for top-down factors. All these theories (Fig. 1) are not necessarily contradictory. Instead, each one adds a missing element necessary to fully understand spring bloom dynamics (Lindemann and St John, 2014). Even though satellite observations have provided great insight over the last decades, the picture is still one of complexity. Our study thus confirms that a single hypothesis for what drives a North Atlantic spring bloom may be too simplistic.
A consensus is yet to be reached regarding the onset of spring phytoplankton blooms in the North Atlantic. Every theory published in the literature claims to best predict the timing of the spring bloom. However, one cannot adopt a single hypothesis simply because all of the theories seem to apply, either on shorter temporal or spatial scales. By revisiting four of the main hypotheses on the subject, we are able to confirm that phytoplankton surface bloom dynamics in the highly variable North Atlantic are far too complex to be driven by the same mechanism in all places and in all years. We show that, in terms of bottom-up processes alone, there is a dominant physical mechanism (mixed-layer shoaling) that best predicts the growing phase of North Atlantic phytoplankton blooms on the mesoscale. However, some regions show coherent patterns, supporting the idea that there are distinct physical phenomena driving spring surface blooms rather than a single one. We believe these findings to be relevant to the ongoing discussion on North Atlantic bloom onset.