Phytoplankton spring bloom initiation: The impact of atmospheric forcing and light in the temperate North Atlantic Ocean

The spring bloom dominates the annual cycle of phytoplankton abundance in large regions of the world oceans. The mechanisms that trigger blooms have been studied for decades, but are still keenly debated, due in part to a lack of data on phytoplankton stocks in winter and early spring. Now however autonomous underwater gliders can provide high-resolution sampling of the upper ocean over inter-seasonal timescales and advance our understanding of spring blooms. In this study, we analyze bio-optical and physical observations collected by gliders at the Porcupine Abyssal Plain observatory site to investigate the impact of atmospheric forcing and light conditions on phytoplankton blooms in the temperate North Atlantic. We contrast three hypotheses for the mechanism of bloom initiation: the critical depth, critical turbulence, and dilution-recoupling hypotheses. Bloom initiation at our study site corresponded to an improvement in growth conditions for phytoplankton (increasing light, decreasing mixing layer depth) and was most consistent with the critical depth hypothesis, with the proviso that mixing depth (rather than mixed layer depth) was considered. After initiation, the observed bloom developed slowly: over several months both depth-integrated inventories and surface concentrations of chlorophyll a increased only by a factor of ~2 and ~3 respectively. We ﬁ nd that periods of convective mixing and high winds in winter and spring can substantially decrease (up to an order of magnitude) light-dependent mean speci ﬁ c growth rate for phytoplankton and prevent the development of rapid, high-magnitude blooms.


Introduction
The annual cycles of phytoplankton in the temperate and subpolar North Atlantic Ocean are characterized by pronounced blooms in spring (Yoder et al., 1993). The timing and intensity of spring blooms may have important consequences for the pelagic ecosystem (Townsend et al., 1994;Platt et al., 2003), sequestration of atmospheric CO 2 in the ocean interior (Martin et al., 2011), surface ocean gas transfer (Codispoti et al., 1982) and ocean temperature (Stramska and Dickey, 1993). Despite their importance, the conditions necessary to trigger phytoplankton spring blooms remain uncertain even after more than 60 years of study (Sathyendranath et al., 2015). To date, three main hypotheses have been proposed: the critical depth hypothesis, critical turbulence hypothesis, and dilution-recoupling hypothesis (Behrenfeld and Boss, 2014).
The critical depth hypothesis (CDH) (Sverdrup, 1953) was the first conventional framework that described the necessary conditions for initiation of phytoplankton spring blooms in the temperate and subpolar North Atlantic. According to CDH, the start of the phytoplankton spring bloom corresponds to shoaling of the ocean mixed layer depth (hereafter z mixed ) above a critical depth (hereafter z cr ), a threshold based on solar radiation, light attenuation in the water column and algal losses from various sources (Smetacek and Passow, 1990). The hypothesis assumes that nutrients are replete during the pre-bloom period due to deep winter mixing and that improving light conditions for phytoplankton is the main factor for triggering spring blooms. Sverdrup suggested that the critical depth criterion could be achieved during seasonal restratification of the upper layer in spring when the mixed layer is rapidly shoaling. This definition of the mixed layer is relevant to the CDH framework if gradients in the vertical distribution of phytoplankton match strong gradients in density profiles. However, subsequent studies (Brainerd and Gregg, 1995) showed that mixed layer depth, conventionally defined using a temperature or density threshold (de Boyer Montegut et al., 2004), is an imperfect proxy for the depth of the layer where mixing is currently active (hereafter z mixing ). If upper ocean mixing is driven by surface cooling, convective cells penetrate to the pycnocline and the entire mixed layer is actively turbulent. However, when surface cooling subsides, mixing of the water column is mainly generated by wind energy. Under wind-driven mixing conditions the layer of active mixing can become shallower than the https://doi.org/10.1016/j.pocean.2019.102202 Received 10 May 2019; Received in revised form 16 September 2019; Accepted 2 October 2019 remnant mixed layer, which may still register diagnostically as the mixed layer depth due to weak vertical gradients in hydrographic profiles. In this case, phytoplankton are likely to be well mixed to z mixing (Chiswell, 2011), as opposed to z mixed . Since the defining mechanism of the CDH is that the bloom is initiated when the phytoplankton are no longer regularly mixed below z cr , we here compare z cr to z mixing , rather than z mixed . Shoaling of z mixing above z cr has been found to be a more precise criterion for the onset of phytoplankton blooms than the traditional critical depth framework and can be achieved before the development of seasonal stratification (Brody and Lozier, 2014;Franks, 2014).
A second proposed mechanism is known as the critical turbulence hypothesis (CTH). According to this hypothesis, the spring bloom can initiate in an arbitrarily deep layer due to changes in mixing intensity rather than in mixing depth (Huisman et al., 1999). Relaxation of turbulence due to weakening atmospheric forcing allows phytoplankton growth near the surface to outpace mixing and a bloom develops, resulting in an uneven vertical distribution of phytoplankton within the mixing layer. The CTH can be expressed in terms of relevant time scales (Taylor and Ferrari, 2011). In this framework, a low mixing rate can result in accumulation of phytoplankton near the surface if where t m is the mixing time scale defined as mixing layer depth divided by a characteristic turbulent velocity scale and t g is the phytoplankton growth time scale. A third mechanism, the dilution-recoupling hypothesis (or the disturbance-recovery hypothesis; hereafter DRH) (Behrenfeld, 2010;Behrenfeld et al., 2013), proposes that decreasing grazing pressure is the main factor controlling bloom onset. According to the DRH, the phytoplankton bloom starts in winter when the ocean surface is cooling and mixing is strong. Despite these conditions accumulation of phytoplankton is possible due to reduced encounter rates of phytoplankton with grazers, i.e. a dilution effect due to winter deep mixing. During the recoupling phase, changes in phytoplankton stocks are determined by the balance between light-dependent phytoplankton specific growth rates and losses due to grazing (Behrenfeld and Boss, 2014). The DRH is supported by observations of positive net accumulation rate of phytoplankton, defined as changes in chlorophyll a (Chl a) inventories derived from satellite ocean colour (Behrenfeld et al., 2013) and Bio-Argo float datasets (Boss and Behrenfeld, 2010) during deep winter mixing and low surface light intensity.
The above outlined hypotheses predict the onset of phytoplankton blooms under different forcing conditions: (i) subsiding ocean surface cooling and/or weakening wind forcing, associated with a reduction in mixing layer depth and/or a reduction in mixing intensity (CDH and CTH) or (ii) periods of strong surface cooling and/or strong wind forcing associated with deep mixing layer depths (DRH). It remains uncertain which conditions trigger spring blooms in the North Atlantic Ocean.
Progress on determining the mechanisms initiating blooms has been hampered by a lack of high temporally and vertically resolved observations across winter and spring. In this regard, autonomous platforms (such as gliders and Bio-Argo floats) represent a powerful tool for studying phytoplankton dynamics due to their ability to obtain frequent depth-resolved profiles of bio-optical and physical properties for long (inter-seasonal) periods of time, even under challenging weather conditions. Changes in surface concentrations of Chl a as well as depthintegrated inventories can be quantified by measuring vertical profiles of Chl a fluorescence (Frajka-Williams et al., 2009;Swart et al., 2014;Thomalla et al., 2015;Bol et al., 2018;Erickson and Thompson, 2018), thus overcoming limitations of satellite ocean colour data. Gliders, in particular, provide a unique opportunity to study the initiation and development of phytoplankton blooms. Firstly, glider data provide a high temporal resolution (up to 6 profiles per day in this study) picture of variability in bio-physical properties. In addition, gliders provide information on the vertical distribution of Chl a, which cannot be achieved from satellite data. Finally, gliders can operate throughout the winter when ship-board operations are difficult and satellite Chl a data are often affected by cloud or low sun angle. Bio-Argo floats also overcome some of the limitations of satellite-derived or ship-based Chl a data, but typically operate with lower temporal resolution than gliders (typically 5-10 days; e.g. Mignot et al., 2016;Boss and Behrenfeld, 2010;Lacour et al., 2019).
In this study we use glider observations and surface forcing from atmospheric reanalyses to study how meteorological and light conditions can affect the onset and development of the spring bloom in the temperate Northeast Atlantic. The objectives of this study are two-fold. First, we analyze which of the three hypotheses (CDH, CTH and DRH) can explain the observed variability in phytoplankton (represented by Chl a fluorescence). Specifically, we aim to answer the following questions: -Can bloom onset be explained by shoaling of z mixing above the critical depth z cr ? (test of CDH) -Can bloom onset be explained by decreasing mixing intensity, which leads to an increase in mixing time scales, t m ? (test of CTH) -Can we detect positive net accumulation rates of phytoplankton during strong cooling of the ocean surface and/or during strong wind mixing before light conditions start to improve? (test of DRH) Second, we examine how atmospheric forcing can influence bloom development after initiation through its effect on mixing layer depth and improving light conditions experienced by phytoplankton.
In the following Data and Methods section, we provide details of the glider missions, a description of the data processing and an overview of the theoretical framework used to evaluate mixing regimes in the water column. The Results section shows the evolution of phytoplankton dynamics from October 2012 to May 2013 and its relation to light conditions and atmospheric forcing. Examination of the three hypotheses and of the impact of atmospheric forcing on bloom development are presented in the Analysis section. The Discussion section relates our results to previous findings and presents our conclusions regarding the processes of bloom initiation and subsequent development in our dataset.

Glider missions
As part of the UK NERC Ocean Surface Mixing Ocean Submesoscale Interactions Study (OSMOSIS), pairs of autonomous ocean gliders were simultaneously deployed 40 km southeast of the Porcupine Abyssal Plain sustained observatory (PAP site; Fig. 1). The full dataset is reported in Damerell et al. (2016) and archived as Damerell et al. (2018). The PAP site is situated between the subtropical and subpolar gyres and is far enough offshore not to be influenced by the circulation over the continental slope to the east; therefore mean flows are weak (Lampitt et al., 2010a). According to a classification of biogeographical provinces in the North Atlantic Ocean (Longhurst et al., 1995), the PAP site is located within the "Westerlies domain". This domain is characterized by deep mixing and nutrient replete conditions in winter representing a suitable site for testing the hypotheses for spring bloom initiation in the temperate North Atlantic Ocean.
During OSMOSIS, the gliders collected measurements over a sampling area between 48.69°N and 48.75°N and between 16.10°W and 16.19°W (Fig. 1). The gliders followed "bowtie" and "hourglass" shaped trajectories ( Fig. 1b; patterns as defined in Alkire et al., 2014). On average each glider provided 10-12 vertical profiles per day (5-6 descents and 5-6 ascents). Over the course of the mission, the gliders maintained the pre-defined trajectories well: 88% of vertical profiles considered in this study were obtained within the intended sampling area (Fig. 1). The gliders were each equipped with a Seabird CT-Sail (conductivity, temperature) sensor, a Wetlabs Triplet ECOpuck (BBFL2VMT configuration measuring Chl a fluorescence and optical backscatter) and a Biospherical Instruments photosynthetically available radiation (PAR) sensor. Vertical sampling resolution of the biooptical sensors ranged from 2 to 10 m depending on battery constraints. The maximum sampling depth was 200 m in autumn and increased to 500 m in winter and spring. PAR sensors sampled every 2 m down to 200 m depth. Further details on the sampling strategy, CT sensor calibration and physical conditions at the PAP site are given by Damerell et al. (2016).

Chl a fluorescence
Raw data from the ECOpuck is output in digital counts. Conversion of the sensor output to scientifically usable units required subtracting dark fluorescence (sensor output in the absence of phytoplankton and light) and multiplying the difference by a scale factor. Dark fluorescence was determined as the median value over the bottom 10 m of each Chl a fluorescence profile. The most frequent value of dark fluorescence for each glider deployment was used for calibration. Time series of dark fluorescence determined from individual glider profiles showed no statistically significant changes in time throughout the glider deployments (95% confidence interval for the regression slopes included zero) indicating that the fluorescence sensors did not drift during each mission.
Erroneous data, such as negative fluorescence readings and values outside the realistic range for the sampling site (> 6 mg m −3 ; threshold defined using Chl a estimates from MODIS-Aqua satellite records), were removed by implementing the quality control procedure described in D' Ortenzio et al. (2010). Daytime fluorescence profiles were corrected for quenching using profiles of optical backscatter following Sackmann (2008). Quenching occurs during periods of high irradiance as a mechanism to protect the photosynthetic apparatus of phytoplankton (e.g. Muller et al., 2001), and is clearly seen as a decrease in fluorescence near the surface (example profiles in Fig. S1). The depth of maximum fluorescence-to-backscatter ratio within the mixed layer was determined for daytime fluorescence profiles (Swart et al., 2014). It was assumed that above this depth a fluorescence profile was affected by quenching. The part of the profile affected by quenching was corrected using the mean nighttime fluorescence-to-backscatter ratio within the upper 20 m. The mean nighttime ratio was calculated from the previous/following night for daytime fluorescence profiles obtained before/after midday. The quenching correction significantly decreased the observed offset between the surface fluorescence for daytime and nighttime profiles (Fig. S1). After applying the quenching corrections, the Chl a vertical profiles were gridded into 5-m vertical bins.
The scale factor for calibration was obtained using in situ Chl a samples collected from CTD Niskin bottles during OSMOSIS process cruises at the start of the glider campaign (October 2012; 64 samples) and near the end of the period considered in this study (April 2013; 135 samples). In situ Chl a concentrations from water samples were determined following the method of Welschmeyer (1994). The scale factor for calibration was determined using linear regression from colocated CTD casts and glider profiles. For each CTD cast, glider profiles collected within a time period of less than a day and within a distance less than 30 km were selected. For each depth of the CTD profile, the corresponding value of glider-measured fluorescence was calculated as the distance-weighted mean of all fluorescence values measured at that depth. The averaged scale factor was applied to calibrate the entire time series. A varying relationship between Chl a concentrations and fluorescence can introduce uncertainties into the calibration of glider fluorometers. However, this study investigates relative changes in Chl a, rather than its absolute magnitude.

Mixed layer depth
Mixed layer depth (z mixed ) estimates were derived from glider profiles of temperature and salinity . The definition of z mixed was based on that of de Boyer Montegut et al. (2004) and calculated using two criteria: a change in temperature of 0.2°C relative to the value at 10 m depth and a change in density of 0.03 kg m −3 relative to the value at 10 m depth. For each profile, the shallower value of the two was chosen for z mixed .

Euphotic depth, light attenuation coefficient and surface PAR
Estimates of the light attenuation coefficient and surface PAR were obtained using the glider factory-calibrated PAR sensor. Euphotic depth, defined here as the depth of 1% of surface irradiance, was estimated from PAR profiles assuming Lambert-Beer's relationship: where K is the vertical attenuation coefficient of irradiance, E 0 is irradiance just below the sea surface and z is depth. By fitting an exponential curve (Eq. (2)) to daytime light profiles, K and E 0 were obtained. By substituting E with 0.01E 0 in eq. (2), the euphotic depth was estimated as 4.6/K. Time of sunset and sunrise were determined using ephem Python module (https://pypi.python.org/pypi/pyephem/, version 3.7.6.0). Following Thomalla et al. (2015), exponential curves were fit only to the part of the profiles above 100 m in order to avoid overfitting to data points in the aphotic layer. Following the fitting procedure, vertical PAR profiles with R 2 < 0.9 (1% of daytime profiles) were excluded from the subsequent analysis. Daily mean surface PAR was obtained from individual glider observations throughout a day using the adjusted sinusoidal interpolation method described by Wang

Satellite-derived data sets and reanalysis data
Daily surface heat flux was obtained from the NCEP/NOAA reanalysis project (Kalnay et al., 1996) with a spatial resolution of 2°. Net surface heat flux was calculated as the sum of net longwave radiation, net shortwave radiation, sensible heat flux and latent heat flux. Heat flux components were extracted for the grid point centered on 48.6°N, 16.8°W, the closest pixel to the sampling site. A time series of wind stress on a regular grid of 0.25°was derived from the Daily Advanced Scatterometer Surface Wind Fields (DASCAT) product (Bentamy and Fillon, 2012). The methodology of wind stress estimation from scatterometer surface wind retrievals was described by Milliff and Morzel (2001). The wind data were extracted for a grid point centered on 48.75°N, 16.25°W.

Turbulence regimes, mixing depth and mixing time scales
We characterize turbulence in convective boundary layers with an applied wind stress by calculating the Monin-Obukhov length scale (Monin and Obukhov, 1954): is the friction velocity, ρ 0 is a seawater reference density, τ is the wind stress, k is von Karman's constant, and = B Qλg c ρ 0 p 0 is the surface buoyancy flux, where Q is the surface heat flux, λ is the thermal expansion coefficient, g is the acceleration due to gravity and c p is the heat capacity of water. A list of variables and constants is provided in Table 1. Previous studies (Schmitt et al., 1989) demonstrated that the haline contribution to the surface density flux in the temperatesubpolar North Atlantic Ocean is an order of magnitude lower than the thermal contribution. This was confirmed at the OSMOSIS site by Thompson et al. (2016). Therefore, we assume that the surface buoyancy flux is predominantly regulated by the surface heat fluxQ.
We classify the effects of convection and wind stress on mixing in the ocean surface layer based on the framework described by Thorpe (2005): Under the wind mixing regime (Case 1), the depth of active mixing (z mixing ) can be shallower than z mixed (Franks, 2014). The depth of wind mixing can extend to the base of the Ekman layer and so z mixing can be scaled as where f is the Coriolis parameter and = C 2 2 is a dimensionless constant. Eq. (4) was implemented to estimate z mixing under Case 1 conditions during negative surface heat flux ( < Q 0). During surface warming ( > Q 0), z mixing is additionally suppressed by a positive buoyancy flux. Zilitinkevich et al. (2002) provided the theoretical framework for scaling of the stably stratified Ekman boundary layer in this situation. The scaling was subsequently implemented in numerical studies of phytoplankton spring blooms, e.g. Enriquez and Taylor (2015), in the following form: are prescribed dimensionless constants. The scaling incorporates an increase of z mixing with increasing u * and decrease of z mixing due to surface heating (increasing B 0 ). Values of C 3 and C 4 were determined by Enriquez and Taylor (2015) using the output of a large-eddy simulation (LES) model.
The mixing time scale can be defined as a ratio between characteristic length and turbulent vertical velocity scales. During the wind mixing regime, the vertical turbulent diffusivity is assumed to be constant down to z mixing and * u represents the characteristic velocity scale. The vertical mixing time scale associated with wind mixing (t m,wind ) can be estimated as: When surface cooling is the main source of turbulence in the water column (Case 2: convective mixing regime), convective cells develop even under weak surface cooling (Taylor and Ferrari, 2011). Under these conditions, the mixed layer depth can be used as the depth of active mixing. In the convective mixing regime, the vertical turbulent diffusivity is assumed to be constant throughout the whole mixed layer and the mixing time scale (t m,convection ) can be estimated as: where = C 1 5 is a prescribed dimensionless constant. The scaling for the turbulent velocity during the convective mixing regime (the denominator in Eq. (7)) was adopted from Deardorff (1972). Sensitivity tests show that the choice of constants (C 1 , C 2 , C 3 , C 4 and C 5 ) in these equations does not affect the main conclusions derived from the analysis ( Supplementary Information; Fig. S4 and S5).

Phytoplankton specific growth rate
Phytoplankton specific growth rate depends on nutrient abundance, light conditions, and temperature. For the glider sampling site, the World Ocean Atlas 2009 gives winter surface nitrate concentrations of 7 µM, concentrations of 5 µM in May, and limiting values (< 1 µM) in June. Observations taken during a cruise to the PAP site in June of the OSMOSIS project found upper ocean (< 10 m) nitrate values to be, on average, 3.4 µM (minimum 0.97 µM). Therefore, the sampling site is  Taylor (2015) and Deardorff (1972).
A. Rumyantseva, et al. Progress in Oceanography 178 (2019) 102202 considered nutrient replete in winter and spring. In the following, we assume that phytoplankton are homogeneously distributed in the mixing layer. The specific growth rate was estimated for the observed temperature conditions assuming a suboptimal light regime and nutrient replete conditions in a similar way to Edwards et al. (2013). The maximum specific growth rate as a function of temperature under abundant light and nutrient conditions was evaluated following Bissinger et al. (2008): where μ max is the maximum phytoplankton growth rate (in d −1 ) and T is the temperature averaged over the mixed layer (in°C). Following Evans and Parslow (1985) the phytoplankton specific growth rate as a function of light was determined as: where α chl is the chlorophyll-specific slope of the phytoplankton-irradiance curve, θ is the cellular chlorophyll-to-carbon mass ratio (a conversion factor between productivity and phytoplankton specific growth rate), and E(z) is the vertical profile of light as defined in Eq.
(2). From Eq. (9), the daily mean division rate for phytoplankton evenly distributed over a layer of depth L can be estimated as: Evaluation of μ involves specification of physiological parameters for phytoplankton, α chl and θ, which are not measured by gliders and can vary depending on phytoplankton physiology and species composition.
To take into account potential variations in α chl and θ, we estimated the probability distribution of μ mean for each day based on the potential ranges of values for α chl and θ, where the typical ranges of values for α chl and θ (Table 1) were taken from Marañon and Holligan (1999) and Sathyendranath et al. (2009) respectively. Based on the probability distribution, the mean and standard deviation of μ mean were determined for each day.

Phytoplankton net accumulation rate
To characterize temporal changes in phytoplankton populations, a time series of net accumulation rate of phytoplankton (r) was constructed. Sustained periods of r > 0 indicate bloom conditions. The net accumulation rate of phytoplankton was calculated using water column integrated Chl a inventories, I(Chl a), and surface Chl a concentrations, S(Chl a), calculated as the mean over glider measurements above 20 m depth. We follow the method described by Behrenfeld (2010), taking into account potential decoupling between the mixed and mixing layers. If z mixing is deepening (z mixing (t 1 ) > z mixing (t 0 )) and z mixing is deeper than the euphotic depth (z mixing (t 1 ) > z eu (t 1 )): and if z mixing is shoaling or z mixing (t 1 ) < z eu (t 1 ): where r is the net accumulation rate over the time interval = − t t t Δ 1 0 . We calculate r using Chl a data averaged in time over 1 day. In this study, we implement a chlorophyll-based approach for calculating phytoplankton net accumulation rate similar to that followed by Boss and Behrenfeld (2010). However, optical backscatter can be considered as another proxy for phytoplankton biomass (Stramski et al. 2004;Boss and Behrenfeld 2010). Our calculations suggest that the net accumulation rate calculated based on Chl a is consistent with that calculated based on the optical backscatter data (R 2 = 0.72, Fig. S3). Consistent temporal patterns of r evaluated using the two optical proxies for phytoplankton biomass indicate that the gliders provided reliable estimates of phytoplankton net accumulation rates.

Critical depth
According to Sverdrup's model, the critical depth, z cr , can be defined by the implicit relationship: where E c is the compensation irradiance. We calculated z cr for two different E c values (0.96 and 1.75 mol m −2 d −1 ; Table 1) previously obtained for the temperate and subpolar North Atlantic Ocean by Siegel et al. (2002) from an analysis of spring bloom timing using satellite and hydrographic data sets. Note that the model in Eq. (13) assumes that phytoplankton growth is linearly proportional to incoming radiation and that the compensation irradiance is constant in time. For the first assumption, we note that for the period under primary study here (i.e. winter to early spring), mixing layer average PAR levels are below typical values of the light saturation parameter for the Northeast Atlantic (~10 mol m −2 d −1 ; Smyth et al., 2004). After mid-April, this assumption is no longer met, however this period is not the focus of our bloom initiation study. For the second assumption, compensation irradiance is unlikely to be constant in time as phytoplankton loss rates vary. The implications of this assumption are considered in the Discussion.

Results
The temporal evolution of Chl a, mixed and mixing layer depths, euphotic depth, atmospheric forcing and surface PAR observed between mid-October and mid-May is shown in Fig. 2. The mixed layer gradually deepened from mid-October until the end of January (Fig. 2a). During this time the surface heat flux was predominantly negative (Fig. 2e). This was also a period of frequent passage of atmospheric fronts associated with strong wind forcing (τ > 0.4 N m −2 ; Fig. 2f) and gradually decreasing surface PAR (Fig. 2d). A convective mixing regime dominated during this period (Fig. 2e). At the beginning of the time-series S (Chl a) and I(Chl a) were approximately 0.7 mg m −3 and 40 mg m −2 respectively (Fig. 2b, 2c). Following the mixed layer deepening in September -December, S(Chl a) and I(Chl a) decreased to 0.1 mg m −3 and 25 mg m −2 respectively in January.
From February until late April the mixed layer remained relatively deep (100-250 m; Fig. 2a) with occasional shoaling and deepening events. The net cooling of the ocean surface significantly subsided; the frequent passage of storms persisted. Conditions of wind mixing and associated divergence between the mixed and mixing layers were more frequent, occasionally interrupted by periods of stronger convective mixing, such as in mid-March (Fig. 2e). Between the 1st of February and the 30th of April, generally positive trends in integrated and surface Chl a were observed coinciding with gradually increasing surface PAR (Fig. 2d). Over this period, I(Chl a) increased by a factor of 2.3, from 30 mg m −2 in February to 70 mg m −2 at the end of April. S(Chl a) increased by a factor of 3, from 0.2 to 0.6 mg m −3 .
The mixed layer remained consistently deeper than the euphotic depth until the end of April when a rapid transition to a shallow stratification was observed. Springtime stratification developed in two phases (Fig. 2a). First, the mixed layer shoaled from 200 m to 50 m in 5 days (19-23 April). Second, the mixed layer deepened below the euphotic zone again for a short period of time (30 April-1 May) and subsequently shoaled again above 50 m. Previous studies at the PAP site suggested that the onset of seasonal stratification can be defined as the date on which the mixed layer shoals above 100 m for more than a week (Lampitt et al., 2010b). According to this definition, during the year of the OSMOSIS mission the onset of seasonal stratification occurred on 22nd April. Strong surface heating and weak wind forcing ( Fig. 2e and f) promoted development of seasonal stratification: surface heat flux generally exceeded 100 W m −2 and wind stress noticeably decreased after 18th April compared with the rest of the time series (τ < 0.2 N m −2 ). During restratification, S(Chl a) initially reached 0.6 mg m −3 and subsequently peaked at 1.5 mg m −3 on 1st May (Fig. 2b). Chl a inventories decreased during this period from 60 mg m −2 to 30 mg m −2 and peaked again at 100 mg m −2 at the beginning of May (Fig. 2b and c).

Analysis of bloom initiation hypotheses
The phytoplankton bloom evolved slowly in weakly-stratified conditions over several months before the onset of the seasonal stratification (Fig. 2a). Below we examine the hypotheses for spring bloom onset and analyze how atmospheric forcing could affect the observed bloom dynamics.
The analysis was performed in a one-dimensional framework, interpreting the observations as a time series following Damerell et al. (2016). Data from both gliders deployed at a given time were used in the analysis. Investigation of spatial heterogeneity within the sampling box is presented in the Supplemental Information, which shows that temporal changes in S(Chl a) and I(Chl a) are consistent between different parts of the sampling area ( Fig. S6 and S7). It is acknowledged that lateral density gradients and associated submesoscale dynamics can drive increased growth of phytoplankton due to short term (< 1 day) restratification of the ocean mixed layer (Mahadevan et al., 2012;Lacour et al., 2017) and significant losses in phytoplankton inventories due to export of organic material along isopycnal surfaces to the ocean interior (Omand et al., 2015). The presence of submesoscale features during the OSMOSIS study is discussed by Thompson et al. (2016). The potential impact of these processes on the sub-daily distribution of phytoplankton is beyond the scope of the study since we aim to describe Chl a variability over inter-seasonal time scales. For an evaluation of the phytoplankton response to short-term restratification events associated with submesoscale dynamics revealed in this dataset, see Erickson and Thompson (2018).

Assessment of the critical depth hypothesis
According to the CDH, improving light conditions and shoaling of the mixed layer above a critical depth prompts the onset of the spring bloom. Here we adapt this definition, for reasons explained in the introduction, to shoaling of the mixing layer above a critical depth (following Brody and Lozier, 2014;Franks, 2014). A comparison of the estimated range of z cr with z mixing is shown in Fig. 3a. Fig. 3b shows the time series of r derived from the glider data using Eqs. (11) and (12). The critical depth criterion (z mixing < z cr ) was generally met from February onwards when a period of mostly positive r was observed (Fig. 3b). In 80% of cases when positive net accumulation rates of phytoplankton were observed, the mixing layer depth was shallower than the estimated critical depth (Fig. 3a). Due to gradually increasing surface PAR (Fig. 2d) the critical depth was sufficiently deep (100-500 m) to allow a net accumulation of phytoplankton to occur for several months before the seasonal restratification (Fig. 3b). Moreover, both the cumulative sum of r (Fig. 3c) and phytoplankton specific growth rate (Fig. 3d) start to increase consistently from the beginning of February. Prior to February, the cumulative sum of r is essentially flat indicating that the ecosystem was in near equilibrium state. This is consistent with an analysis of net community production derived from the gliders' oxygen dataset which suggests that net autotrophy begins in early February (Binetti et al., this issue). Our results indicate that the period of positive net accumulation of phytoplankton was associated with improving light conditions for phytoplankton. Therefore, the observed accumulation of phytoplankton in deep mixed layers is broadly consistent with the critical depth hypothesis, provided that mixing depth, rather than mixed layer depth, is considered.

Assessment of the critical turbulence hypothesis
According to the CTH, a spring bloom can start when phytoplankton growth time scales are shorter than mixing time scales. The test of this hypothesis was conducted by comparing estimated mixing (Eqs. (6) and (7)) and growth time scales (t m and t g respectively). We used the values of specific growth rate near the surface ( = μ z ( 0); Eq. (9)) to estimate the minimum growth time scale ( = = t μ z g,min 1 ( 0 ) ) that can be achieved under the observed light conditions. In reality, the growth time scales are also affected by various loss factors, therefore > t t g g ,m i n . For both wind and convective mixing regimes, turbulent mixing time scales are approximately an order of magnitude smaller than calculated growth time scales (Fig. 3e). This suggests that for the observed meteorological conditions, the critical turbulence criterion for bloom initiation expressed in terms of relevant time scales (i.e. < t t g, min m ) was not met. The scaling arguments used here to assess the critical turbulence hypothesis carry their own uncertainties, however the order of magnitude difference we find between the growth and mixing time scales gives confidence that our general conclusion is robust.

Assessment of the dilution-recoupling hypothesis
DRH associates the bloom onset with decreasing loss rates in winter due to dilution of zooplankton in deep mixed layers. Positive net accumulation rates during the periods of deepest mixing can indicate the dilution effect and provide support for the hypothesis (as previously shown by Behrenfeld, 2010). From December until the end of January, the convective mixing regime dominated and deepening of the mixed layer was observed ( Fig. 3a and f). During this time, average net accumulation rates were close to zero (r = 0.02 d −1 (standard deviation = 0.1 d −1 ; Fig. 3b). The cumulative sum of r is relatively constant during December and January (Fig. 3c), the period when mixing was the deepest and when phytoplankton specific growth rate was low (Fig. 3d). Positive net accumulation rates were mainly observed from early February onwards and corresponded to gradually improving light conditions and increasing μ mean (Fig. 3d). There is no clear evidence that bloom onset occurs due to the dilution effect in the OSMOSIS glider dataset.

Impact of mixing regimes on the bloom development
Our dataset demonstrates that positive r was detected before upper ocean restratification in spring when mean specific growth rates for phytoplankton started to increase due to improving light conditions. Of the three hypotheses examined, CDH best explains the observed variability in Chl a (Fig. 3d). However, the observed increase of S(Chl a) and I(Chl a) was only a factor of 3 and 2 respectively from February (when r > 0 was detected) to the end of April (when seasonal stratification developed). Therefore, the question remains: why did shoaling of z mixing above z cr not result in a rapid and pronounced phytoplankton bloom? Even though a bloom can be defined as an onset of net growth, some studies (e.g. Platt et al., 1991) consider the rapid accumulation of biomass an essential signature of phytoplankton spring blooms. We now discuss how atmospheric forcing over the winter-spring period and the associated mixing regimes in the ocean boundary layer influenced μ mean and bloom progression. The depth of active mixing determines light conditions experienced by phytoplankton cells and influences μ mean (Eq. (10)). μ mean was calculated using Eq. (10) for mean surface PAR (20 mol m −2 day −1 ) and mean light attenuation coefficient (0.066 m −1 ) observed during February-April (Fig. 4a). When the mixing depth is shallower than the euphotic depth, μ mean increases abruptly. Mean specific growth rates evaluated for the wind and convective mixing regimes as a function of Fig. 3. Time series of (a) z cr , z mixing and z mixed , (b) net accumulation rate of phytoplankton (r ) calculated using Chl a data from the gliders, (c) cumulative sum of r , (d) mean specific growth rate (μ mean ) over z mixing (gray shaded area shows an estimated range of uncertainty associated with the choice of values for α chl and θ), (e) growth time scales (t g ) and mixing time scales for the convective (t m,convection ) and the wind (t m,wind ) mixing regimes, (f) surface buoyancy flux B 0 . Separate figures for three time periods (15th of October -15th of January, 15th of January -15th of March and 15th of March -15th of May) are shown in the Supplementary materials (Figs. S8-S10).
wind stress and surface heat flux are shown in Fig. 4b. During the convective mixing regime, convective cells penetrate the whole mixed layer resulting in relatively low μ mean (0.1-0.3 d −1 ). Under the wind mixing regime μ mean is low (0.1-0.5 d −1 ) for relatively strong wind forcing (τ > 0.2 N m −2 ) and significantly higher (μ mean > 0.5 d −1 ) for calm weather conditions (τ < 0.2 N m −2 ). Therefore, in weakly stratified conditions, shoaling of the mixing layer can significantly increase μ mean when wind forcing is weak.
Transition to a wind mixing regime can have a two-fold effect on phytoplankton inventories. First, division rates are significantly increased for the part of the community trapped within z mixing due to increased light exposure. Second, algae are trapped within the remnant layer, below z mixing . A decaying mixing intensity below z mixing increases the residence time of phytoplankton within the remnant layer (Franks, 2014) and potentially below the euphotic zone, where conditions are unfavorable for phytoplankton growth. As an example, if z mixed is 250 m, shoaling of z mixing to 50 m can lead to 80% of the population being trapped in the aphotic zone (if z eu < 50 m) and being permanently lost from the surface layer. However, the enhanced growth near the surface rebuilds phytoplankton inventories at the same time.
To demonstrate this, we use a simple model for phytoplankton accumulation: where P is phytoplankton concentration at time t and depth z. Eq. (14) A. Rumyantseva, et al. Progress in Oceanography 178 (2019) 102202 omits vertical diffusion and assumes that the initial vertical distribution of phytoplankton is depth-independent over the mixed layer. In the case of an actively turbulent deep mixed layer, the evolution of phytoplankton concentration at any depth within z mixed can be described by the following equation: where μ zmixed is the average division rate in z mixed . When the turbulence structure changes under the wind mixing regime, only the part of the community within z mixing grows, albeit with a higher specific growth rate (μ zmixing ) due to increased light exposure: where P mixing is phytoplankton concentration within z mixing . If < z z mixing eu , phytoplankton also accumulates between z mixing and z eu , although this region is not actively mixed and, therefore, the specific growth rate is not uniform. An example of a vertical profile of µ demonstrates changes in the specific growth rate with depth (Fig. 5a). Fig. 5b demonstrates the estimated changes in phytoplankton inventories for the range of z mixing values observed in the glider data (25, 50, 75, 100 and 150 m) as well as for z mixing = z mixed (=250 m) (Fig. 5b). In the case of the shallowest z mixing , the fastest increase of phytoplankton concentration occurs, i.e. weak wind forcing conditions result in rapid phytoplankton growth near the surface.
To illustrate the effect of mixing regimes on phytoplankton inventories using in situ data, we contrast Chl a variability at the end of February, when net surface heat flux was approaching zero and wind speed was low, with a period of strong convective mixing in March (the selected periods are marked as E1 and E2 respectively in Fig. 2). During E1 and E2 vertical profiles of temperature are relatively uniform within z mixed (Fig. 6b and Fig. 6e). The first period, E1, is characterized by a significant difference between mixed layer depth (= 250 m) and mixing layer depth (z mixing = 25 m) (Fig. 6a). Averaged vertical profiles for E1 show surface intensified vertical distribution of Chl a within the hydrographically defined mixed layer (Fig. 6a). Phytoplankton are relatively well mixed down to mean z mixing during E1. For E2, when mixing is driven by convection, the Chl a distribution is relatively uniform within the mixed layer (Fig. 6d). During E1, phytoplankton inventories gradually decrease below the euphotic depth (Fig. 6c), because the divergence between z mixing and z mixed during the wind mixing regime significantly increases the residence time of phytoplankton in the aphotic zone prompting phytoplankton losses. The opposite effect is observed for the phytoplankton population within the euphotic zone (Fig. 6c). As a result, overall phytoplankton inventories increase slightly (Fig. 6f). For E2, the changes in water column integrated Chl a correspond to integrated Chl a below the euphotic zone (Fig. 6f), as the vertically homogeneous turbulent mixed layer results in relatively uniform light conditions for phytoplankton cells.
Thus a convective mixing regime is generally associated with low μ mean for the phytoplankton community. Our results imply that the shift to a wind mixing regime can significantly increase μ mean , especially in the case of weak wind forcing, but at the same time a significant part of the phytoplankton community can be trapped in the aphotic layer and potentially lost. Intermittent mixed layer restratification by submesoscale activity, which is active in the study region (Erickson and Thompson, 2018), may play a similar role to intermittent wind forcing in bloom development by giving rise to rapid shoaling of the mixing layer.

Discussion
Simultaneous physical and biogeochemical observations from ocean gliders have been used to study the impact of atmospheric forcing on phytoplankton spring bloom dynamics and to test which of three commonly discussed hypotheses for bloom initiation in the North Atlantic (CDH, CTH, or DRH) can best explain the observed phytoplankton dynamics. The glider observations do not support the hypothesis that blooms initiate in midwinter prior to the seasonal increase in light as predicted by DRH. In a test of the CTH, a comparison of mixing and growth time scales indicated that decreasing mixing intensity was unlikely to be driving enhanced phytoplankton growth in winter and spring, consistent with our observation that the mixed layer does not permanently restratify until early May. Instead, our data suggest that the positive net accumulation rates of phytoplankton were mainly observed when the shoaling mixing layer became shallower than the estimated critical depth threshold, and mean phytoplankton growth rate was gradually increasing due to improving light conditions. Our analysis therefore supports the CDH, with the proviso that mixing, rather than mixed layer, depths are considered. Our analysis also showed that seasonal patterns of bloom development are shaped by atmospheric forcing through their effect on mixing layer depth and light conditions experienced by phytoplankton.
Of the hypotheses considered in this study, the observed phytoplankton variability was most consistent with the CDH framework. However, it is important to note that the estimates of critical depth depend greatly on the value of compensation irradiance. Sverdrup (1953) assumed the compensation irradiance to be a constant value that reflected the loss rates of phytoplankton due to respiration, grazing, sedimentation and other factors during the pre-bloom period. In reality, the compensation irradiance is likely to be a dynamic parameter that varies depending on grazing pressure and other loss factors. Slow spring bloom development can cause an immediate response in the grazing community as shown by Waniek (2003) using a mixed-layer model coupled with an NPZD model and discussed by Behrenfeld and Boss (2014). In this regard, the framework explaining seasonal phytoplankton variability involves some aspects of both the CDH and the DRH. In particular, initiation of phytoplankton biomass accumulation is driven by improving light conditions, consistent with the CDH. However, the subsequent development of the bloom can be affected by the balance between light-dependent specific growth rates and loss rates that may be variable. The latter violates the assumptions of CDH regarding constant losses and is better explained by the "recoupling" part of the DRH. Potential variability in loss rates can, in part, be addressed by investigation of the phytoplankton spring bloom through a modeling framework (e.g. Lévy, 2015) since gliders do not provide measurements of grazing pressure.
In their examination of the critical depth framework, Platt et al. (1991) concluded that "the Sverdrup criterion is necessary but not sufficient" for rapid phytoplankton accumulation in spring. Simply put, the criterion can only indicate if the net growth of phytoplankton can occur, but not how rapidly the bloom will progress. Platt et al. (1991) suggested that the frequent occurrence of storms prevents the rapid accumulation of phytoplankton biomass. In this study, the analysis of the mean specific growth rates under the wind and convective mixing regimes shows how meteorological conditions can affect the development of the spring bloom. For the convective mixing regime, phytoplankton accumulation occurs over the entire mixed layer, but relatively slowly due to low mean specific growth rates. When a shift to a wind mixing regime takes place, the part of the phytoplankton population within the mixing layer grows more rapidly, but at the same time, losses from the remnant layer can slow down the vertically-integrated accumulation of phytoplankton biomass. It is important to note that this source of losses is not included in the critical depth model. Rapid growth near the surface can build phytoplankton inventories rapidly (over about 4 days; Fig. 5). High winds can interrupt the rapid development of the phytoplankton spring bloom. Interestingly, the effect of wind mixing on phytoplankton blooms in the North Atlantic Ocean is different for spring and autumn. In autumn, phytoplankton growth is limited by nutrient availability (Martinez et al., 2011) and intensified mixing associated with an autumnal storm can deliver nutrients to the euphotic layer and trigger enhanced phytoplankton growth (Rumyantseva et al., 2015). Here we show that windy conditions in spring can prevent the development of rapid phytoplankton spring blooms by decreasing light-dependent mean specific growth rates for the phytoplankton community.
The correlation between spring bloom characteristics and wind conditions in spring has been noted before. Analysis of satellite data (Ueyama and Monger, 2005;Henson et al., 2009) has shown late, low magnitude phytoplankton spring blooms in the North Atlantic during the positive phase of the North Atlantic Oscillation, commonly associated with strong westerly winds in winter-spring. Waniek (2003) demonstrated that windy weather in spring results in low magnitude interrupted phytoplankton blooms, similar to the one captured during the OSMOSIS mission. The passage of weather systems varies interannually and might be affected by future changes in the North Atlantic climate (Gillett et al., 2003). Predicted changes include increasing sea surface temperature (Allen et al., 2014), increased net surface heat flux and increasingly positive North Atlantic Oscillation conditions (Osborn, 2004) that would change basin-scale wind forcing patterns. In this study, it has been shown that atmospheric forcing and associated mixing regimes have a profound impact on phytoplankton growth rates and the development of algal blooms.

Concluding remarks
Autonomous underwater gliders deployed at the PAP site provided a unique data set capturing the development of the 2012 phytoplankton bloom in the temperate Northeast Atlantic. Motivated by the longrunning debate around Sverdrup's critical depth hypothesis, this study concludes that the bloom onset corresponded to improving light conditions and was mainly consistent with the critical depth hypothesis, provided the divergence between the mixed layer depth and the active mixing layer is considered. The subsequent development of the bloom was affected by the meteorological conditions through their effect on the light environment experienced by phytoplankton. The observed low magnitude of the bloom was explained by the frequent passage of high winds and periods of convective mixing after the seasonal onset of net phytoplankton growth.

Declaration of Competing Interest
No conflicts of interest.  (e) show mean vertical profiles of temperature (black solid line) and standard deviation (dashed lines) during E1 and E2 respectively. Red, blue and green lines show z mixing , z mixed , and z eu respectively. Panels (c) and (f) show absolute changes in Chl a stocks over the whole water column (black line), above the euphotic depth (red line) and below the euphotic depth (blue line) during E1 and E2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)