Basin-wide mechanisms for spring bloom initiation: how typical is the North Atlantic?

Harriet S. Cole1,2, Stephanie Henson3, Adrian P. Martin3, and Andrew Yool3 Ocean and Earth Science, National Oceanography Centre Southampton, University of Southampton Waterfront Campus, European Way, Southampton SO14 3ZH, UK Now at Marine Laboratory, Marine Scotland Science, 375 Victoria Road, Aberdeen AB11 9DB, UK National Oceanography Centre, University of Southampton Waterfront Campus, European Way, Southampton SO14 3ZH, UK *Corresponding author: tel: +1 224 295508; fax: +1 224 295511; e-mail: h.cole@marlab.ac.uk


Introduction
Phenology is the study of the timing of periodic biological events, such as the annual phytoplankton bloom, and has led to a number of ecological and biogeochemical insights. For example, the timing of greatest food availability is important for grazers in addition to food abundance as summarized in the match-mismatch hypothesis (Cushing, 1990). This hypothesis states that interannual variability in the timing of the bloom results in years where the bloom coincides with larval hatching (a match) and years when their timing is not synchronous (mismatch). This affects larval survival and recruitment rates of several commercially important species (Cushing, 1990; increase as the longer period of stratification extends the growing season with blooms starting earlier and ending later (Doney, 2006). In contrast, productivity in the subtropics is expected to decrease as nutrient limitation will be more severe and longer lasting. Thus, bloom timing metrics can be used as additional monitoring indicators to detect changes in the pelagic ecosystem . However, before trying to predict these changes we must first understand the drivers of contemporary interannual variability.
Several recent studies have related the observed variability in bloom initiation to annual, seasonal, or shorter term means in a physical parameter. In the North Atlantic, bloom initiation is reported to vary by 2 -3 weeks relative to the mean and has been related to the winter mean net heat flux (NHF) and windspeed (Henson et al., 2006) as well as the windspeed during the bloom (Ueyama and Monger, 2005). Variability in initiation timing has also been linked to large-scale climate indices such as the North Atlantic Oscillation (Henson et al., 2009;Zhai et al., 2013). In the Arctic, bloom initiation has been related to the timing of sea ice melt, changes in which have driven an advance in bloom timing of 50 d since 1997 (Kahru et al., 2011).
In the North Pacific, several drivers of bloom timing variability have been identified. Sasaoka et al. (2011) found bloom initiation varied by 5 weeks in the coastal subpolar North Pacific and was related to the winter Southern Oscillation Index. In El Niño years (negative SOI), sea surface temperature (SST) in these regions is warmer than average consistent with earlier stratification and alleviated light limitation resulting in an earlier bloom. Conversely, in open ocean regions of the North Pacific blooms occurred earlier in La Niña years when SST was cooler, mixing stronger, and iron and nutrient entrainment greater. Additionally, bloom timing variability has been linked to large-scale climate indices such as the Pacific Decadal Oscillation and changes in SST and mixed layer depth (MLD; Sasaoka et al., 2011;Chiba et al., 2012).
Bloom initiation variability in the Southern Ocean has only recently been examined, and seasonal and intraseasonal variability in MLD, heat fluxes, and iron supply have been identified as potential drivers (Thomalla et al., 2011). Other regions have been shown to be more variable, especially boundary regions between subpolar and subtropical gyres. In these areas, bloom initiation dates vary by 20 weeks in the North Atlantic (Henson et al., 2009) and Southern Ocean (Thomalla et al., 2011).
Studying the drivers of variability in bloom initiation may also provide insights into the mechanisms that lead to the onset of the spring bloom in the first place. Additionally, using interannual variability in the timing of the bloom is likely to be a stronger test of theories of bloom initiation than using a single year, which may have experienced anomalous conditions. Currently, there are several theories concerning the conditions that cause the bloom to start in subpolar regions. The critical depth hypothesis states that the bloom starts when the MLD shoals to a point where the average irradiance in the mixed layer is high enough for net growth to occur in the phytoplankton population (i.e. photosynthesis . respiration; Sverdrup, 1953). An extension of this, the critical turbulence hypothesis states that blooms may start when mixed layers are still very deep, if the rate of near surface vertical mixing has reduced sufficiently. This means that phytoplankton are no longer rapidly mixed out of the euphotic zone and are able to accumulate and grow in the near surface (Huisman et al., 1999;Taylor and Ferrari, 2011). Turbulent mixing may be reduced through an increase in the heat flux into the ocean or a reduction in windspeed. A further theory is the dilution-recoupling hypothesis, which states that the bloom starts when the mixed layer is at its maximum depth, when phytoplankton and grazers are sufficiently diluted so encounter rates are minimal and growth can overcome grazing losses (Behrenfeld, 2010;Boss and Behrenfeld, 2010). On smaller spatial scales, eddies have been observed to play a role in starting the spring bloom by restratifying the water column in spring, earlier and more quickly than by solar heating alone (Mahadevan et al., 2012).
Significantly, all these hypotheses were developed using observations from the North Atlantic. However, each of the subpolar basins present their own character and are hydrodynamically different from each other; the North Atlantic has very deep winter mixed layers, the North Pacific has a permanent halocline and is iron limited, as is the Southern Ocean. Thus, although the theories discussed above are mechanisms that may be present across all the basins, there may not be a single globally consistent mechanism for bloom initiation. This begs the question of whether the dominant mechanism for bloom initiation is the same everywhere, or if our view of spring bloom initiation is skewed by the dominance of studies focused on the North Atlantic.
This study aims to use satellite-derived chlorophyll data and output from a data-assimilating model to quantify variability in the date of bloom initiation in the subpolar North Atlantic, North Pacific, and Southern Oceans. The bloom initiation dates will be used to identify dominant, basin-wide drivers of variability in the onset of the bloom and to assess if bloom timing responds to interannual variability in forcing in a similar manner. In this way, we will address the question, "How typical is the North Atlantic?"

Datasets
Bloom initiation timing was calculated from a satellite-derived chlorophyll dataset, GlobColour, and from the output of the NASA Ocean Biogeochemistry Model (NOBM). Both bloom metric datasets were regridded to 1 × 18 grid to match the physical datasets. The GlobColour dataset was available at 1 × 18 resolution and as 8 d composites and was downloaded from http://globcolour. info. It combines data from three ocean colour sensors: SeaWiFS, MODIS, and MERIS. Data covering the period 2002-2009 are used here as all three sensors were operational over this period and to match the availability of the other datasets being used. GlobColour is reported to perform better than the individual sensors when compared with the SeaWiFS Bio-optical Archive and Storage System (SeaBASS) database, an in situ dataset used for comparison with remote sensing products (Durand, 2007). For GlobColour chlorophyll, the match up statistics with in situ data are similar to SeaWiFS with an root mean square log error of 30% (SeaWiFS: 31%) and a coefficient of determination of 0.73 (SeaWiFS: 0.76; Gregg and Casey, 2004;Durand, 2007). However, GlobColour improves coverage in both space and time, since it includes three datasets (Durand, 2007).
NOBM is a dynamic plankton ecosystem model that assimilates SeaWiFS data as surface chlorophyll. The ability of NOBM to reproduce seasonal blooms in chlorophyll was assessed by Cole et al. (2012). NOBM shows high fidelity to seasonal characteristics, although absolute concentrations of chlorophyll can be slightly different (Cole et al., 2012). This means that while the magnitude of the bloom in NOBM may be different from that recorded by SeaWiFS, the timing of the bloom remains the same in both datasets. Since the bloom timing is based on the relative change in chlorophyll concentration, accuracy in the timing rather than the magnitude is of greater consequence for this study.
The analysis was performed on both the GlobColour and NOBM chlorophyll datasets. The GlobColour dataset extends further towards the poles giving wider global coverage and spans more years ( ) than NOBM (2002( -2007. However, NOBM is used here to more accurately determine interannual variability in bloom timing (as it is a gap-free dataset vs. the satellite-derived GlobColour which has missing data), as well as to provide a second estimate and thus improve the robustness of this study. The chlorophyll output from NOBM has previously been used to quantify the impact of missing data in the satellite record on the accuracy of bloom initiation and peak dates (Cole et al., 2012). Cole et al. (2012) found that missing data resulted in errors of 30 d in bloom initiation timing further justifying the use of NOBM output in this analysis. All the physical datasets used were observational and were used in conjunction with bloom timing dates from both GlobColour and NOBM.
The physical parameters considered were NHF, MLD, and mean mixed layer photosynthetically active radiation (PAR ML ). The NHF data were downloaded from the WHOI OAflux project website (http://oaflux.whoi.edu) as a 1 × 18 dataset and with daily temporal resolution, which was averaged to 8 d to match the chlorophyll data. The MLD was calculated from a global dataset of regularly gridded temperature and salinity profiles, obtained from the Coriolis project (http://www.coriolis.eu.org). These profiles were collected from ARGO floats, expendable bathythermographs, conductivity temperature and depth, profiles and moorings. These were available on a 1/10 by 1/108 grid but were regridded to 1 × 18 to match the other datasets. Density was calculated from the regridded profiles and the MLD was defined as the depth at which the density had increased by 0.03 kg m 23 relative to a reference depth of 10 m to avoid diurnal effects. PAR data were downloaded from the MODIS ocean colour webpage (http://oceancolor.gsfc.nasa.gov) as 8 d composites and were regridded to a spatial resolution of 1 × 18. The attenuation coefficient at 490 nm was obtained from the GlobColour webpage and is derived from the chlorophyll concentrations (Barrot, 2010). This was converted to Kd PAR (attenuation coefficient for all PAR wavelengths) using the conversion detailed in Morel et al. (2007). The mean irradiance in the mixed layer was calculated using Equation (1).
where PAR 0 is the surface photosynthetically available irradiance, k the attenuation coefficient (Kd PAR ), and D the depth of the mixed layer as described above. Subpolar regions were identified using correlations between the chlorophyll and MLD climatological annual cycles at each pixel (as in Henson et al., 2009). Pixels poleward of 408N/S with a negative correlation (shallower MLD associated with higher chlorophyll) were defined as being subpolar. Pixels without a strong seasonal cycle, defined as having a coefficient of variation (CV 2 variance/ mean) value of 0.35 or lower, were excluded from the analysis (Cole et al., 2012). Defining the subpolar regions was necessary as some subtropical regions have a strong seasonal cycle but the mechanisms underlying it are different from subpolar blooms.

Timing of bloom initiation
The initiation date was defined as the date on which the chlorophyll concentration exceeded a threshold value more than 5% of the annual median. This criterion has been frequently used previously to define the bloom initiation date (Siegel et al., 2002;Henson et al., 2009;Racault et al., 2012;Brody et al., 2013). As January is not always an appropriate date from which to start a "bloom year" (e.g. in the southern hemisphere), the date of the peak (defined as the date of the maximum chlorophyll value in each year) was used as a reference point and the date of bloom initiation is found in the preceding 6 months. The bloom initiation was found by searching backwards in time from the peak identifying the last datum before the threshold value was exceeded. The bloom initiation date was calculated for each year at each pixel in the North Atlantic, North Pacific, and Southern Ocean.

Physical timing metrics
Many of the previous studies discussed in the introduction that address bloom initiation processes focus on seasonal or annual means in parameters such as MLD, NHF, or windspeed. Here, the focus will be on the relationships between interannual variability in bloom initiation and the timing of a change in the physical environment (e.g. timing of MLD shoaling). Strong correlations between the interannual variability in bloom timing and physical forcing may lead to alternative insights into the mechanisms that control bloom initiation.
A suite of physical timing metrics (see Table 1 for definitions) were devised based on extant theories. The timing metrics were calculated for each year at each pixel. In the Southern Ocean, timeseries were adjusted to run from July to June before calculating the metric dates to match the seasonal chlorophyll cycle. The date the gradient in PAR ML became positive. This must occur between the beginning of the annual cycle and the date of the peak value of PAR ML Fastest rise in PAR ML The date of largest increase in PAR ML between the beginning of the annual cycle and the maximum PAR ML value MLD becomes shallower than euphotic depth (MLD , Zeu) Date on which the MLD becomes shallower than the euphotic depth (define as 1% of surface irradiance) for at least 16 d (two time-steps). This served as a proxy for the MLD crossing the critical depth and light levels ceasing to be limiting Sverdrup (1953) Basin-wide mechanisms for spring bloom initiation From the annual cycle of MLD the dates of maximum MLD and MLD shoaling were calculated. From the NHF annual cycle, the date the NHF became positive was calculated. This latter metric is a proxy for a change from winter convection to less turbulent conditions (Taylor and Ferrari, 2011). The dates of interest from the annual cycle of PAR ML include the date PAR ML starts to increase and the date of fastest rise in PAR ML . These are proxies for when light levels become favourable for growth and when light is most rapidly becoming less limiting. Additionally, the depth of the euphotic zone, defined as the depth at which 1% of the surface PAR still remained, was calculated. This was used to calculate the date on which the MLD became shallower than the euphotic zone depth as another proxy for the alleviation of light limitation.
Latitudinal means of the physical and bloom timing metrics were calculated for each degree of latitude across the three basins. In addition, each ocean basin was split into a number of 108 longitude by 108 latitude boxes, for which interannual anomalies in the physical and bloom metric dates were calculated (Figure 1). The size of these boxes was chosen to filter out mesoscale features and to focus on large-scale variability. Averaging the metric dates over these large boxes also removes the effect of temporal autocorrelation on the correlations. Similarly, the effect of spatial autocorrelation is reduced for the correlations of the latitudinal means due to the coarse resolution averaging used. For all basins, the anomalies in bloom and physical metric dates were approximately normally distributed.

Results
Generally, the date of bloom initiation in all three basins becomes progressively later towards the poles (Figure 1). In the North Atlantic, blooms initiate during February-March at 40 -458N progressing to initiation dates in May, June, and July at higher latitudes. In the North Pacific, although on average blooms occur later at higher latitudes, blooms are seen to initiate earlier (April/May) near the coast progressing to later blooms offshore (June/July). In the Southern Ocean, bloom timing gets progressively later towards the pole, with blooms starting in August/September at 408S and January/February south of 658S.
As an example of the degree to which the timing of the bloom and physical metrics coincide, average time-series from the North Atlantic (45-808N 260 to 08E), North Pacific (50-708N 2120 to 1208E), and Southern Ocean (280 to2508N 2180 to 1808E) are shown in Figure 2. Broadly speaking, bloom initiation occurred just after the NHF turns positive, the MLD was shoaling and PAR ML was increasing in all three basins except for the Southern Ocean. There are clear differences in the variability of the physical parameters between the three basins though the NHF annual cycle was fairly similar between the three basins ranging from 2200 to 200 W m 22 . The largest interbasin differences were seen in the maximum depth reached by the mixed layer. In the North Atlantic, the mean maximum MLD varied between 150 and 200 m was 50 m at its deepest in the North Pacific and 110 m in the Southern Ocean. The maximum The critical depth and critical turbulence hypothesis both indicate that if a physical timing event is a driver of spatial variability in bloom timing then the two metrics are likely to occur at similar times in a given year, though potentially with a lag, and have similar latitudinal gradients ( Figure 3; Table 2, see Supplementary Figure S2 for plots of further physical metrics). In the North Atlantic, the date of bloom initiation closely follows the progression of the date the NHF turns positive and progresses polewards at a similar rate (initiation: 5.72 km d 21 , NHF turns positive: 6.72 km d 21 ) and they are strongly correlated (r ¼ 0.95). The bloom initiation also has very similar timing to the date the ML irradiance starts to increase and the date of MLD shoaling. Both are strongly correlated (r ¼ 0.87 and r ¼ 0.86, respectively). The other physical timing metrics are also highly correlated with bloom initiation (most have r ≥ 0.67 and are statistically significant; Table 2) in the North Atlantic. However, in the North Pacific, the metrics do not match well and are generally weakly correlated (r ≤ 0.32), except the date at which the MLD becomes shallower than the euphotic depth (r ¼ 0.52). In the Southern Ocean, the date of bloom initiation has a similar progression to that of the date the NHF turns positive, is highly correlated (r ¼ 0.98), and progresses towards the pole at similar speeds (initiation: 22.54 km d 21 , NHF turns positive: 22.53 km d 21 ). Furthermore, the date the PAR starts to rise is also highly correlated with bloom initiation and progresses at a similar speed ( On an interannual basis, mostly weak correlations are seen between the interannual anomalies in bloom initiation and in the physical timing metrics. However, a reasonable correlation (r ¼ 0.61) is seen between interannual variability in the timing of bloom initiation and the date at which the NHF turns positive in the North Atlantic ( Figure 4, Table 3, see Supplementary Figure  S3 for scatterplots of further physical metrics). Weak but significant (p , 0.05) correlations with bloom initiation timing are also seen with the date of maximum MLD (r ¼ 0.21) and MLD shoaling (r ¼ 0.25), MLD . Zeu (r ¼ 0.29), and PAR starts to increase (r ¼ 0.29) in the North Atlantic. No strong positive correlations were seen in the North Pacific. The only significant, positive correlation seen in the Southern Ocean was with the date of NHF turning positive but was very weak (,0.12).

Drivers of variability in bloom initiation
In the North Atlantic, the timing of bloom initiation coincides with the period of increasing stratification and decreasing mixing as the NHF becomes positive, the MLD shoals, and PAR ML increases ( Figure 2). This suggests that they may potentially be prominent drivers for interannual variability in bloom initiation. This is supported by the synchronous timing and high correlation in the latitudinal progression of bloom initiation, NHF turning positive, MLD shoaling, and the date the PAR ML starts to increase.
In comparison, bloom initiation in the North Pacific had neither high correlation, nor synchronous timing with, any of the physical timing metrics. However, it was moderately correlated and had synchronous timing, with the date the MLD becomes shallower than the euphotic depth. This suggests that variability in the date of bloom initiation is weakly driven by the alleviation of light limitation as the MLD shoals and the surface irradiance increases as spring progresses.
Similar to the North Atlantic, the date the NHF turns positive and the PAR ML starts to increase are indicated as drivers of bloom timing variability in the Southern Ocean. Though many of the other metrics in the North Atlantic and Southern Ocean were moderately to highly correlated, not all had synchronous timing across the whole basin.
It should be noted that at higher latitudes in the Southern Ocean the uncertainty in bloom initiation date due to missing data was Basin-wide mechanisms for spring bloom initiation PAR ML begins to increase for each of the three subpolar basins. The light grey shaded area associated with the mean bloom initiation date represents the uncertainty in bloom initiation date arising from missing data in the chlorophyll time-series from which bloom initiation is calculated (Cole et al., 2012). The dark grey (blue in online colour version) shaded regions around the mean physical timing dates represent the variability (1 s.d.) in the physical timing metric dates at that latitude. Pixels with shallow bathymetry (,200 m) or without a significant seasonal cycle (coefficient of variation ,0.35) were removed before longitudinal averaging. The vertical axis in all panels is the average metric date (dd-mmm). A colour version of this figure is available online. high. However, uncertainty in the date may also arise from the definition used to calculate bloom initiation. The definition used here for bloom initiation was found to be similar to other definitions of bloom timing (Supplementary Figure S1). Not all the metrics identified as strong drivers of spatial variability were strong drivers of interannual variability. The majority of the physical timing metrics were only weakly correlated with interannual variability in bloom initiation, if at all. The strongest relationships were seen in the North Atlantic where the date that the NHF turns positive was seen to be a strong predictor for the bloom initiation date (Figure 4). The date that the NHF turns positive explains 37% of the variability in bloom initiation date. In comparison, the relationship between NHF turning positive and bloom initiation was very weak in the North Pacific and Southern Ocean. The other metric identified as a consistent driver of spatial variation in bloom initiation was the date PAR ML starts to rise in both the North Atlantic and Southern Ocean. However, only a weak (if significant) relationship with interannual variability in bloom initiation was found in the North Atlantic, whereas there was no relationship between these metrics in the other two ocean basins. Other weak correlations were seen with the date of MLD shoaling, maximum MLD, and MLD , Zeu, though only in the North Atlantic.

Spatial vs. interannual variability in bloom onset
Interestingly, many of the dominant drivers of latitudinal variability are not found to be dominant drivers of interannual variability  Figure 1; each box has 6 datapoints. Fitted lines are shown for statistically significant correlations and are generated using a type II regression which accounts for variability in both x and y parameters. and even then not consistently across the basins. Similarly, in the subpolar North Atlantic, Follows and Dutkiewicz (2002) found that the relationships between proxies for mixing (e.g. bloom period NHF and windspeed) were strongly correlated with bloom period chlorophyll concentration spatially, but not on an interannual basis. Follows and Dutkiewicz (2002) proposed that this lack of correlation was due to a number of factors such as changes in insolation, mesoscale variability, and grazing pressure having a stronger influence on interannual variability in bloom magnitude. Thus, the lack of an obviously dominant driver of interannual variability in bloom initiation may be because a combination of variables is responsible but that each, individually, only has a small influence on the initiation date or that noise in the datasets used here obscure the interannual relationship. There are many variables not included in this set of timing metrics, especially those representing ecological/biogeochemical processes such as grazing or micronutrient availability which may also have an influence on the interannual variability in timing of bloom onset. Alternatively, a lack of strong relationships in both spatial and interannual variability may indicate that these physical timing metrics are not drivers of bloom timing.

Mechanisms for subpolar bloom initiation
Two of the main theories for bloom initiation are the critical depth hypothesis (Sverdrup, 1953), for which MLD shoaling is a proxy, and the critical turbulence hypothesis (Taylor and Ferrari, 2011), for which NHF turning positive is a proxy. Both the date of NHF turning positive and MLD shoaling were strongly related to bloom initiation, in a latitudinal sense, in the North Atlantic. This was also seen in the Southern Ocean though the date of MLD shoaling showed a weaker relationship than in the North Atlantic. However, only the North Atlantic shows a basin-wide interannual relationship between the NHF turning positive and bloom initiation. This nevertheless suggests that the critical turbulence hypothesis is a better description of the processes that lead to a bloom than the critical depth hypothesis, at least for the North Atlantic.
To investigate the role of MLD and NHF in bloom initiation further, the conditions in the North Atlantic at the time of bloom initiation were examined in more detail ( Figure 5 and Table 4). The median conditions in the North Atlantic when blooms started were a low and positive NHF ( 31 W m 22 ), a relatively shallow MLD  ( 80 m), and low PAR ML ( 3 W m 22 ). In accordance with the critical turbulence hypothesis, a number of blooms (19%) were seen to start in deep mixed layers (deeper than median value) but when the NHF was low and positive ( 25 W m 22 ) ( Figure 5). Conversely, 19% of the blooms were also seen to start when MLDs were shallow (shallower than median value) with NHF mostly ranging from 0 to 100 W m 22 . The greatest number of blooms (51%) occurs when the NHF ranges from 50 to 150 W m 22 and the MLD is between 40 and 240 m deep (the red pixels in Figure 5a). Additionally, PAR ML at the start of the bloom was seen to range from very low values, indicating deep MLD, to very high values when the MLD was very shallow. Based on these results, it would seem that there is evidence to support both the critical depth and critical turbulence hypotheses. However, though MLDs are shallow when the bloom starts the NHF is also typically close to zero. This suggests that the critical turbulence hypothesis could also be valid where both a reduction in turbulence and MLD shoaling come together to initiate a bloom. The distribution in PAR ML values allows further examination of the critical depth hypothesis. The compensation irradiance is the irradiance at which photosynthesis is equal to respiration. Assuming that at bloom initiation, the MLD is equal to the critical depth, the mean mixed layer irradiance (PAR ML ) will be equal to the community compensation irradiance, as this is the light level needed for phytoplankton growth to overcome community losses and rapidly increase phytoplankton standing stocks, according to the critical depth hypothesis (Sverdrup, 1953).
The large range (0.17 -37.52 W m 22 ) in PAR ML at bloom initiation across the North Atlantic suggests that the critical depth hypothesis is not a complete description of the conditions necessary for initiating a bloom. This is because bloom initiation does not seem to be dependent on the mean mixed layer irradiance reaching a certain threshold level that is sufficient for net population growth as bloom initiation occurs under a large range of PAR ML values. Blooms starting in low PAR ML suggest that phytoplankton are not necessarily mixed throughout the whole mixed layer, due to low turbulence mixing, so that their exposure to sufficient irradiance levels for net growth is increased relative to if they were evenly distributed throughout the ML. Conversely, some blooms initiated when PAR ML was higher than average. This suggests that another limiting factor delayed the start of the bloom since the mean mixed layer irradiance would not be expected to be limiting at these levels.
In the North Pacific, there is weak evidence supporting the critical turbulence hypothesis. Interannual relationships between NHF turning positive and bloom initiation are very weak and bloom initiation occurs in a small number (8%) of pixels when NHF is low and positive ( 25 W m 22 ) and the MLD still relatively deep (Figure 5b). Furthermore, the median value of NHF at bloom initiation is much higher than in the Atlantic ( 122 W m 22 ), MLD much shallower ( 21 m), and PAR ML much higher ( 11 W m 22 ). This suggests that for the majority of the basin the initiation of the bloom is delayed compared with the North Atlantic, occurring later in the season when ocean warming is more advanced, the MLD shallower and as a result the PAR ML higher. This is supported by the higher NHF values (.100 Wm 22 ) and shallower MLD (,33 m) under which the majority of blooms in the North Pacific occur (red pixels in Figure 5b). Furthermore, this is reflected in the map of mean bloom initiation dates in Figure 1.
In Basin-wide mechanisms for spring bloom initiation 2037 neither the MLD shoaling nor NHF turning positive were strongly correlated with interannual variability in the bloom initiation date.
Overall, the NHF turning positive is a stronger predictor for the date of bloom initiation than the shoaling of the MLD but this relationship is seen to dominate the whole basin only in the North Atlantic. The strength of the relationship between NHF turning positive and bloom initiation in the North Atlantic is again obvious in its spatial pattern (mean correlation ¼ 0.60) (Figure 6) where strong and significant positive correlations occur across the basin. Conversely, there is no coherent pattern in the correlations between the onset of positive NHF and bloom initiation in the North Pacific. In the Southern Ocean, strong and significant positive correlations occur in larger patches near the Antarctic landmass. In all three basins, the correlations with the onset of positive NHF are more likely to be stronger and positive than the correlations with the shoaling of the MLD.

How typical is the North Atlantic?
Our results raise the question: why do the other basins behave differently from the North Atlantic? Several theories exist for the differences in bloom dynamics between the basins. The bloom in the North Pacific occurs later in the year compared with the North Atlantic. This has been previously reported and several hypotheses have been proposed to explain this observation (Heinrich, 1962;Miller, 1993;Fasham, 1995;Banse and English, 1999;Boyd and Harrison, 1999). Some have focused on the role of grazing on bloom development. The life history of the dominant copepod species in the Pacific means that adults emerge from winter hibernation and immediately lay their eggs, ensuring that their larvae are ready to feed once the bloom starts, whereas in the North Atlantic the adult copepods need to feed on the bloom before being able to lay eggs (Parsons and Lalli, 1988). The differences in the dominant copepod between the basins results from a combination of environmental factors and the differences in phytoplankton production and species composition (Parsons and Lalli, 1988). Alternatively, the elevated light levels resulting from the relatively shallow winter MLD in the North Pacific may sustain primary production over winter. This allows microzooplankton populations to be maintained at levels high enough to graze down the bloom as soon as it starts (Evans and Parslow, 1985;Fasham, 1995;Boyd and Harrison, 1999).
Another theory has suggested that iron limitation may play a role in limiting phytoplankton growth in spring in both the North Pacific and Southern Ocean. Both of these oceans are high nutrient-low chlorophyll regions, and so processes that deliver iron to these regions may dictate the start date of the bloom by increasing growth rates to a point where phytoplankton can escape grazing pressure (Landry et al., 1993;Banse and English, 1999).
The role of wind mixing has not been assessed here but is another parameter that modulates turbulence and may also be a significant driver of bloom timing. A reduction in windspeed may lead to a drop in turbulent mixing and a shoaling MLD and has previously been reported as a predictor for bloom onset around New Zealand and in the Irminger Basin (Henson et al., 2006;Chiswell et al., 2013).
Finally, large differences in the physical environments in the three basins have been observed here. For example, there is a deeper winter MLD in the North Atlantic, although the NHF annual cycle is broadly the same across the basins. The North Atlantic is saltier than the North Pacific due to the northward advection of warm and salty water from the subtropics by the Atlantic Meridional Overturning Circulation (Schmittner et al., 2005). This makes it easier for cooling to overturn the water column, giving a deeper MLD in the Atlantic.
The North Atlantic may in fact be a poor model system for theoretical understanding of bloom initiation. This is because the North Atlantic results presented here identify the onset of positive NHF as a basin wide driver of variability in bloom initiation though a similar result is not seen in the North Pacific or Southern Ocean. Furthermore, the timing of MLD shoaling was uncorrelated in the North Pacific and Southern Ocean but showed a weak relationship with bloom initiation timing in the North Atlantic. Thus, it would seem that observations from the majority of the world's subpolar regions do not support either the critical turbulence or critical depth hypothesis and that the North Atlantic response may be unique. Further evaluation of the current theories on bloom initiation in the North Pacific and Southern Ocean is needed to properly explain the observed differences from the North Atlantic. This could include the differences in the character of atmospheric forcing, underlying oceanographic conditions (e.g. saltier vs. fresher), physical heterogeneity across the basins, seasonal differences in micronutrient availability and in grazing pressure. This may lead to a greater understanding of the mechanisms that lead to the onset of the spring bloom.
To definitively identify the mechanisms that initiate the spring bloom, in situ measurements of additional parameters (e.g. turbulence, nutrients, grazing rates) would be needed. Alternatively, future work could examine the role of mechanisms in different regions using appropriately validated, high resolution, free-running models. Additionally, three-dimensional spatial scale processes at the mesoscale may also influence bloom timing. For example, the results of Mahadevan et al. (2012) demonstrate the role of frontal slumping as a mechanism for restratifying the water column and initiating the bloom. Furthermore, the exact mechanisms for the survival of phytoplankton over winter are still not known in great detail which, if examined further, may help to identify mechanisms for bloom initiation. One such theory is that although phytoplankton are convectively mixed down to great depths they are also returned to the surface euphotic zone with sufficient frequency that production can be sustained over winter (Backhaus et al., 1999;Backhaus et al., 2003). Alternatively, the mixed layer eddies observed by Mahadevan et al. (2012) may sustain production in winter by intermittently stratifying some of the water column.

Limitations
This study has used global datasets of NHF and MLD but it is important to note that the datasets are not "equal" in their relationship with the bloom timing metrics. The NHF dataset is measured synoptically from satellite and reanalysis products with each pixel being an areal average. Conversely, the MLD was calculated from sporadic, non-uniform point samples of temperature and salinity. Thus, they cannot be easily ranked against each other because the way the datasets are derived may result in a limit on the variability that one dataset can show. Furthermore, since many of the physical metrics are interdependent, it is difficult to separate out the effects of just one on bloom timing. As the NHF starts to increase from its winter minimum, the MLD begins to shoal as well because the reduction in cooling and mixing encourages restratification of the water column (Huisman et al., 1999). The results here show that the MLD shoaling occurs, on average, within +20 d of the NHF turning positive. Thus, it is difficult to definitively say which mechanism, MLD shoaling or a reduction in mixing, starts the bloom. In a comparison study of three different bloom timing metric 2038 H. S. Cole et al. definitions, Brody et al. (2013) demonstrated that all three gave bloom initiation dates that were approximately synchronous with the NHF turning positive. This suggests that our result is robust (to some extent) to the choice of bloom metric definition, although a slightly different period and different datasets are used here (note that the range in bloom initiation date with latitude calculated in this study from different metric definitions can be seen in Supplementary Figure S1).

Conclusion
In conclusion, evidence has been presented that the critical turbulence hypothesis is the most likely mechanism for bloom initiation in the North Atlantic. The influence of this mechanism appears much weaker in the North Pacific and Southern Ocean though smaller areas of the Southern Ocean did show strong correlation with NHF becoming positive. This lack of consistency in bloom timing response across these three subpolar regions indicates that the North Atlantic is not a universal model system for developing a general theoretical understanding of the mechanisms that lead to the onset of the spring bloom. Further investigation of the environmental and ecological differences between the subpolar basins may lead to a greater understanding of the environmental controls on bloom timing.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.