Phytoplankton Phenology in the North Atlantic: Insights From Profiling Float Measurements

Phytoplankton division rate (μ), loss rate (l), and specific accumulation rate (r) were calculated using Chlorophyll-a (Chl) and phytoplankton carbon (Cphyto) derived from bio-optical measurements on 12 Argo profiling floats in a north-south section of the western North Atlantic Ocean (40° N to 60° N). The float results were used to quantify the seasonal phytoplankton phenology and bloom dynamics for the region. Latitudinally varying phytoplankton dynamics were observed. In the north, the CPhyto peak was higher, occurred later, and was accompanied by higher total annual CPhyto accumulation. In contrast, in the south, stronger μ-r decoupling occurred despite smaller seasonal variations in mixed layer depth (suggesting the possibility of other ecological forcing), and was accompanied by an increasing portion of winter to total annual production, consistent with relief of nutrient limitation. The float observations of phytoplankton phenology for the mixed layer were compared to ocean color satellite remote sensing observations and found to be similar. A similar comparison to an eddy-resolving ocean simulation found the model only reproduced some aspects of the observed phytoplankton phenology, indicating possible biases in the simulated physical forcing, turbulent dynamics, and bio-physical interactions. In addition to seasonal patterns in the mixed layer, the float measurements provided information on the vertical distribution of physical and biogeochemical quantities and therefore are complementary to the remote sensing measurements. Seasonal phenology patterns arise from interactions between “bottom-up” (e.g., resources for growth) and “top-down” (e.g., grazing, mortality) factors that involve both biological and physical drivers. The Argo float data are consistent with the disturbance recovery hypothesis over the full, annual seasonal cycle; for the late winter/early spring transition, the float data are also consistent with other bloom hypotheses (e.g., critical photosynthesis, critical division rate, and meso/sub-mesoscale physics) that highlight the importance of brief, episodic boundary layer shoaling for decoupling of division and grazing rates.


INTRODUCTION
Seasonal cycles in phytoplankton productivity and biomass vary significantly across the global ocean, especially in high latitude regions where strong seasonal variability occurs in environmental conditions (Yoder and Kennelly, 2003;Longhurst, 2007). Organic carbon produced by phytoplankton photosynthesis is a major source of the energy supporting the marine food web as well as an important part of the global carbon cycle (Takahashi et al., 2009). Studies of phytoplankton seasonal timing or phenology, therefore, have been a major focus in oceanography and ecology. The sub-polar North Atlantic region has been of particular interest because of the massive phytoplankton biomass increase (bloom) that occurs there each winter into spring (Obata et al., 1996;Dale et al., 1999;Siegel et al., 2002). A commonly used identifier for the period when phytoplankton blooms develop is the phytoplankton accumulation rate (r) becoming positive, when the specific division rate (µ) exceeds the specific loss rate (l, r = µ − l > 0) (Behrenfeld, 2010).
Historically the North Atlantic spring bloom was explained with the "Critical Depth Hypothesis, " which includes the original "Critical Photosynthesis Hypothesis" and the modified "Critical Division Rate Hypothesis" (Gran and Braarud, 1935;Sverdrup, 1953;Smetacek and Passow, 1990). Both hypotheses assumed, for simplicity and field data limitations, that l is a constant and independent of µ, and as a result the bloom initiation was viewed as a "bottom-up" process, focusing on limiting factors (light) for µ rather than factors influencing l (e.g., grazing, mortality). A modified version of "Critical Division Rate Hypothesis" was established later as the "Critical Turbulence Hypothesis, " which used a definition of "active mixed layer" instead of a temperature/density-based mixed layer depth (Huisman et al., 1999;Taylor and Ferrari, 2011b). Based on satellite observations, Behrenfeld (2010) suggested that in the North Atlantic phytoplankton biomass accumulation starts in early winter before the mixed layer starts shoaling when one calculates biomass accumulation from the rate of change in the vertically integrated phytoplankton population that begins to rise prior to increases in surface phytoplankton concentration. These findings led to the "Dilution-Recoupling Hypothesis" (Behrenfeld, 2010). Furthermore, studies have also shown the impacts of meso/sub-mesoscale physics on the phytoplankton bloom (e.g., Taylor and Ferrari, 2011a;Mahadevan et al., 2012;Mahadevan, 2016;Lacour et al., 2017;Rumyantseva et al., 2019).
Beyond the studies on the spring bloom, efforts have also been made to address annual phytoplankton dynamics (e.g., Evans and Parslow, 1985;Chiswell et al., 2015). In the past decade, advances in satellite remote sensing techniques made it possible to perform long-term monitoring of phytoplankton dynamics, greatly extending our knowledge of the annual phytoplankton phenology (e.g., Behrenfeld, 2010;Behrenfeld and Boss, 2018). For example, with satellite remote sensing data and the Biogeochemical Element Cycling-Community Climate System Model (BEC-CCSM), the "Dilution-Recoupling Hypothesis" was further developed as the "Disturbance and Recovery Hypothesis (DRH)" (Behrenfeld andBoss, 2013, 2018; in the context of the phytoplankton annual cycles. In the case of the subpolar and temperate North Atlantic, according to DRH, the mechanism behind the annual phytoplankton dynamics can be explained as follows (Behrenfeld et al., 2019). During late spring or early summer, loss rates l begin catching up with division rate µ, due to enhanced grazing in the stratified shallow mixed layer. When µ approaches its annual maximum, µ and l move toward near-equilibrium (r = µ − l ≈ 0) and phytoplankton biomass stabilizes or starts to decrease (r = µ − l < 0) (equilibrium phase). From late summer into autumn, the mixed layer starts deepening and the light-level decreases, causing decreases in both µ and r (depletion phase, but a fall bloom could happen due to increased nutrient supply and/or dilution effect from the deepening mixed layer). In the winter despite the deep mixed layer, low light-level and low phytoplankton division rate (µ), phytoplankton biomass (r) starts to accumulate (bloom is initiated) because the deep mixing reduces the encounter rate between phytoplankton and grazers and therefore reduces l (dilution phase). In the following spring, with the shoaling mixed layer and increasing light level, µ increases and stays ahead of l, which maintains positive r and continues phytoplankton biomass accumulation (accumulation phase). The accumulation phase is often associated with a rapid increase in surface phytoplankton biomass, the phenomenon that led to the traditional interpretation of spring bloom timing.
With "state-of-the-art" satellite remote sensing approaches, the upper ocean (e.g., euphotic layer, or mixed layer) is often treated as a single box so that the phytoplankton dynamics in this box can be derived from satellite remote sensing measurements at the ocean surface. Recent developments in autonomous ocean sensor platforms (e.g., gliders, Argo profiling floats) provide an important complementary approach to satellite remote sensing for high-resolution (e.g., higher sampling frequency, vertical profiling capability) observation and regions with high solar zenith angles or obscured by clouds. Bio-optical measurements on Argo profiling floats have been utilized for the North Atlantic in recent years. For example, Boss and Behrenfeld (2010) used the data from an Argo float between 48 • N to 52 • N, and a more recent study from Mignot et al. (2018) averaged the data from 7 floats between 52 • N to 65 • N to form one climatological dataset. Both studies showed phytoplankton biomass accumulation in the winter before the spring bloom for the vertically integrated phytoplankton population. However, the earlier Argo-based studies in this area focused on the northern region and lacked information regarding north-south variations and the vertical structure of phytoplankton dynamics within in the region.
New float observations allow us to build-on and extend previous studies of phytoplankton seasonal phenology in the western North Atlantic. In 2016 and 2017, 12 Argo floats were deployed during the North Atlantic Aerosols and Marine Ecosystems Study (NAAMES) 1 funded by the National Aeronautics and Space Administration (NASA) (Behrenfeld et al., 2019), thereby providing an unique opportunity to study the phytoplankton dynamics and to further test the DRH along a north-south transect from 40 • N to 60 • N. In this study, the research area was divided into four sub-regions and climatologies of phytoplankton carbon (C phyto ), division rate (µ), loss rate (l), and specific accumulation rate (r) for each sub-region were derived using bio-optical measurements on Argo floats from 2016 to 2018. The Argo-derived mixed layer phytoplankton phenology was also compared to estimates from satellite remote sensing observations and a model simulation. Furthermore, the vertical distribution of phytoplankton carbon production was studied with the Argo data.

Float Data
Fluorescence-based chlorophyll-a concentration (Chl), phytoplankton carbon biomass (C phyto ), salinity (S), temperature (T), and pressure (P) data used for phytoplankton metrics calculation were obtained from 12 biogeochemical Argo profiling floats deployed during the NAAMES expeditions (Figure 1). Chl was calibrated against discrete samples and corrected for non-photochemical quenching. C phyto was calculated from the float-measured backscattering (following the conversion of Graff et al., 2015 and assuming a spectrally varying power-law exponent of particulate backscattering of −0.78 when extrapolating from 700 nm, see the Supplementary Material for details). The float data and documentation for detailed processing protocols are available at: http://misclab.umeoce.maine.edu/floats/. The depth resolution was ∼2 m in the upper 500 m and ∼4 m from 500 m to 1000 m, and the profiling frequency for each float varied from 1 to 5 days. For each float, if multiple profiles were obtained in 1 day, only the data from the profile at dawn were used for that day. In this work, the float data were binned into 5-day time bins corresponding to the lowest profiling frequency.
Phytoplankton specific division rate, µ (d −1 ), for the mixed layer was calculated as (Behrenfeld et al., 2005): where Chl and C phyto were the mixed layer mean Chlorophylla concentration (mg m −3 ) and the mixed layer phytoplankton carbon biomass (mg m −3 ), respectively, and I g was the daily mixed layer median light level (mol photons m −2 h −1 ), which was given by: Mixed layer depth (MLD) (m) was defined by a density offset from the value at 10 m using a threshold of 0.03 kg m −3 (de Boyer Montégut, 2004), where the density profile is computed using float-measured temperature (T) and salinity (S). A satellite product (Modis-Aqua, as described in section "Satellite Data") for surface Photosynthetically Active Radiation (PAR surf ) was used to compute I 0 due to the fact that high-quality float-measured PAR was not available. The diffuse attenuation coefficient at 490 nm (k 490 , m −1 ) used in Eq. 2 was calculated from float-measured mixed layer Chlorophyll-a concentration (Morel and Maritorena, 2001): As pointed out by previous studies, the Chl/C phyto ratio may vary seasonally due to community variations (Cetinić et al., 2015;Schallenberg et al., 2019), therefore it should be noted that the empirical relationships used to derive µ from Chl/C phyto may not hold for all seasons and regions. Mixed layer phytoplankton specific net accumulation rate, r (d −1 ), was calculated from temporal changes in C phyto between two time points : when MLD is deepening and > Z(0.415) (4) when MLD is shoaling or < Z(0.415) where C phyto and C phyto were the mixed layer phytoplankton carbon vertical inventory (mol C m −2 ) and concentration (mol C m −3 ), respectively. The isolume depth Z(0.415) (m) was defined as the depth below which light is insufficient for photosynthesis (I = 0.415 mol photon m −2 d −1 ) (Letelier et al., 2004;Boss and Behrenfeld, 2010). Eq. 4 accounted for the scenario when mixed layer deepening dilutes the phytoplankton population with phytoplankton-free water from below. Because they were based on net biomass changes, r values calculated from Eqs. 4 and 5 account for both biological processes and physical transport. r is designed such that winter-time dilution by entrainment of waters from below the mixed layer is not interpreted as a biological loss.

Satellite Data
Satellite-based (MODIS-Aqua) estimates of surface phytoplankton carbon biomass, phytoplankton specific division rate, MLD, and surface PAR from 2015 to 2018 were obtained from the Oregon State University Ocean Productivity Website. 2 The satellite products were compiled using the Carbon-based Productivity Model (CbPM, Behrenfeld et al., 2005;Westberry et al., 2008), with a spatial resolution of 0.167 • by 0.167 • and a temporal resolution of 8 days. Satellite-estimated µ and C phyto data were used directly in this analysis, while the mixed layer r values were calculated from C phyto with Eqs. 4 and 5. It should be noted that MLD could not be directly measured by satellites, and the satellite-based CbPM model uses MLD data from the HYbrid Coordinate Ocean Model. 3

Eddy-Resolving Ocean Model Output
The eddy-resolving ocean physical-biogeochemical model output from Harrison et al. (2018) was used for our analysis. The model simulation was conducted using the ocean component of the Community Earth System Model version 1 (CESM1, Hurrell et al., 2013), with a spatial resolution of 0.1 • . The CESM ocean physics component, known as the Parallel Ocean Program (POP v2, Smith et al., 2010), was coupled with the Biogeochemical Elemental Cycle (BEC) model (Moore et al., 2013). Harrison et al. (2018) documented in detail the procedures used for the physical and biogeochemical model initialization, forcing, spin-up, and model integration, which was deemed by the authors sufficient to examine the effect of mesoscale eddies on seasonal phytoplankton dynamics and large-scale biogeographical patterns relative to a lower resolution control simulation. We used the final year of the 5-year, fully coupled (prognostic physics and biogeochemistry) from the CESM eddy resolving simulation for comparison with the float and satellite data.
The simulated biogeochemical tracers, including phytoplankton, are advected and mixed by the POP prognostic physics and affected by biogeochemical source and sink terms. The model output includes the carbon fixation rate, carbon biomass (total carbon), and chlorophyll for small phytoplankton, diatom, and diazotroph functional groups, as well as the loss from aggregation, mortality, and grazing. Phytoplankton division rates are determined from temperature, light, and nutrient (N, P, Fe, and for diatoms additionally Si) levels following standard model equations and functional forms (Moore et al., 2013). Loss processes from the surface layer include zooplankton grazing, mortality, and downward physical transport associated with advection and turbulent mixing. The model tracks a single zooplankton group that grazes differentially on the functional groups; diatoms are treated as a larger size class generating enhanced vertical export via zooplankton grazing and aggregation. The phytoplankton growth, mortality, and grazing parameters also differ across the functional groups and were chosen to encourage elevated diatom fraction under bloom conditions because of decoupling of diatom growth and grazing under appropriate light and nutrient conditions. More details on the BEC model equations included in the Supplementary Material.
The model mixed layer r (d −1 ) was calculated from temporal changes in the simulated total carbon biomass between two 5-day averages using Eqs. 4 and 5, and model mixed layer µ and l (both in d −1 ) were calculated as the ratios of total carbon fixation and total carbon loss to total phytoplankton carbon, respectively.
where the total carbon fixation and total carbon were sums of carbon fixation and carbon biomass for small phytoplankton, diatom, and diazotroph, and total carbon loss included the grazing and mortality terms, and for diatoms the aggregation term, of the above-mentioned three functional groups. The biogeochemical model admittedly has weaknesses and by construction must abstract many aspects of a complex plankton ecosystem. However, the level of detail in the simulation of trophic interactions is in line with other basin to globe ocean biogeochemical models (Hashioka et al., 2013). Over most of the year and in particular during bloom periods, phytoplankton losses tend to be dominated in the model by zooplankton grazing, reflecting the seasonal ramp up of zooplankton biomass and reduced grazing limitation term at high prey concentrations. The phytoplankton mortality term contributes to a smaller background loss throughout the year and is comparable to the low grazing losses during the winter deep convection periods. The simulated patterns of phytoplankton-zooplankton seasonal phenology and grazing losses are qualitatively similar to reported dynamics from limited process studies, but we lack sufficient direct information from field studies to constrain grazing patterns at basin and seasonal scales. In previous studies with a coarse-resolution simulation of the CESM model,  demonstrated that the model formulation of net phytoplankton growth and loss exhibited considerable success in capturing the seasonal phenology in the North Atlantic as found in satellite observations, providing some confidence in the model grazing formulation. The simulation dynamics (phytoplankton growth, grazing, and other loss processes) also are broadly consistent with the underlying mechanisms of the DRH , thus motivating the comparison of the high-resolution simulation to the Argo float data.
The model outputs used in this work should be viewed as a statistical representation of the system and are not expected to exactly match in situ observations because: (1) the model used "Normal Year" forcing which does not match synoptic atmosphere forcing events of a particular year or include interannual variability; and (2) the eddy resolving model has its own turbulent dynamics that do not match the specific ocean eddy field encountered during the field observations. Although the model may not correspond to the exact phytoplankton phenology shown by the Argo and satellite observations, it does provide an internally consistent framework and is thus a useful reference for comparison. The identification of when and where the model deviates from observed patterns in phytoplankton bloom dynamics and seasonal phenology also offers a basis for prioritizing future model development and refinement efforts.

Construction of Regional, Climatological Patterns in Bloom Phenology
Based on geographic locations and data availability of the float measurements, the study area was divided into four geographic regions (D1-D4, see Figure 1). Seasonal climatologies for the float data were created for each parameter with the following number of Argo profiles by region: 249 for D1, 635 for D2, 381 for D3, and 259 for D4. A monthly climatology of r profiles was also created with Argo float data for each region to study the vertical distribution of phytoplankton carbon net production and loss.
For the satellite and model data, climatologies were based on the mean values of each region, which provide a better statistical representation given the substantial mesoscale variability within the D1-D4 regions because: (1) the float measurements were not evenly distributed (spatially or temporally); and (2) there were significantly less satellite data available if we only chose the satellite data along the float trajectories. Averaging over the regional boxes was also required for the model output because, by construct, details of the eddy-resolving fields did not match the actual situation along the observed float trajectories given the use of climatological forcing and the model internally generated turbulent field.

Criteria Used to Identify Disturbance Recovery Hypothesis (DRH) Phases
We identify the time periods for the four disturbance recovery hypothesis (DRH) phases using three parameters: (1) the time rate of change in MLD (dMLD/dt), which indicates whether the mixed layer is deepening or shoaling; (2) the time rate of change in normalized C phyto (1/C phyto * d C phyto /dt), which indicates the change of mixed layer phytoplankton concentration, and (3) r, which indicates specific net phytoplankton accumulation rate according to Eq. 4 [when mixed layer is deepening and MLD > Z(0.415)] and Eq. 5 [when mixed layer is shoaling or MLD < Z(0.415), in this case r equals 1/C phyto * d C phyto /dt]. The criteria used to identify the characteristics of each DRH phase are listed in Table 1. For the equilibrium phase (phase A) beginning in late spring or summer, the mixed layer was shallow and stable (dMLD/dt ≈ 0) and relatively high µ was balanced by enhanced grazing or other loss mechanisms (l) in the stratified shallow mixed layer. Therefore, both 1/C phyto * d C phyto /dt and r moved toward near-equilibrium (near zero), and they were mostly identical. For the depletion phase (phase B) in late summer to autumn, the mixed layer 1 | Summary of the criteria used to identify the phytoplankton seasonal phenology characteristics described by the Disturbance Recovery Hypothesis (DRH) based on temporal variations of mixed layer depth MLD (entrainment versus shoaling), mixed layer phytoplankton biomass C phyto , and specific net phytoplankton accumulation rate r from Eq. 4 and 5.

DRH phase Criteria
started deepening (dMLD/dt > 0), µ decreased due to lower light levels, and phytoplankton concentration started decreasing (1/C phyto * d C phyto /dt < 0). Meanwhile, as l also decreased (due to dilution effect), the mixed layer biomass r could be negative or near-zero (depending on whether µ or l decreased faster). As the mixed layer continued deepening, phase C in the late autumn or winter was mainly controlled by the dilution effect, when phytoplankton concentrations were stable or decreasing (1/C phyto * d C phyto /dt < 0 or ≈ 0) but the integrated mixed layer biomass increased (r > 0). In the following spring (phase D), a shoaling mixed layer (dMLD/dt < 0) and increasing light level caused µ to exceed l, which led to an increase in phytoplankton concentration and an acceleration in the biomass accumulation rate (1/C phyto * d C phyto /dt = r > 0). An analysis based on the relationship between dµ/dt and r in the context of DRH (based on Eq. 8 in Behrenfeld and Boss, 2018) was also included in the Supplementary Material.

RESULTS AND DISCUSSION
Observed/Modeled Spatial-Temporal Variability in Surface Photosynthetically Active Radiation (PAR surf ) and Mixed Layer Depth (MLD) Surface Photosynthetically Active Radiation (PAR surf ) and mixed layer depth (MLD) are two important environmental parameters essential for characterizing phytoplankton growth and biomass accumulation. As shown by the satellite data, surface PAR for phytoplankton growth increased from the north to south (Figures 2A, 3A, 4A, 5A), with the PAR surf seasonal maximum (from 1.5 to 2 mol photon m −2 h −1 ) around May to July (late spring/summer) and the minimum (from 0 to 0.3 mol photon m −2 h −1 ) in January. Only the northernmost D1 region had near-zero PAR surf in January. The PAR surf from CESM (estimated as a constant fraction of solar short-wave heat flux) showed similar trends with slightly lower light level in the spring and summer months. For all four regions, MLD from all three approaches started deepening and shoaling around September and March, respectively (Figures 2B, 3B, 4B, 5B). The northernmost D1 region had the deepest winter MLD at about 500 m, while the winter MLD for other regions were around 200 to 300 m. For most of the year the isolume depth (below which light is insufficient for photosynthesis) Z(0.415) was shallower than MLD except for June to September when they were similar. Z(0.415) values increased from north to south as PAR surf increased.
The MLD from all three approaches (float, satellite, and model) had similar seasonal patterns and geographic trends, although values from CESM were significantly shallower in the winter. From the perspective of the DRH, the general similarity in the timing of entrainment versus detrainment suggests that all three approaches are capable of capturing the proposed effect of dilution on decoupling production and loss terms, though with weaker magnitude in the model. Model MLD discrepancies in this region have been identified previously and are not unexpected given the complex, regional pattern of deep and shallow winter MLD linked to the positions of the Gulf Stream/North Atlantic

Observed/Modeled Spatial-Temporal Variability in Phytoplankton Carbon Biomass (C phyto ), Growth Rate (µ) and Carbon Accumulation Rate (r)
The overall magnitude of year-round mixed layer C phyto decreased from the north to the south (Figures 2C, 3C, 4C, 5C, also see Figure 6), with the peak occurring in the summer for the northern regions (D1 and D2) and late spring for the southern regions (D3 and D4). Winter biomass (the inoculum for the spring bloom) made a greater contribution to annual biomass in the southern regions, as the summer/spring C phyto decreased from north to south while winter C phyto increased from north to south. As indicated by the standard deviation (light shadings), spatial variations also decreased from north to south, indicating that the northern regions were more dynamic than the southern regions (Glover et al., 2018). In general, the overall magnitude of year-round mixed layer µ increased as latitude decreased and PAR surf increased (Figures 2D, 3D, 4D, 5D). The peak of µ from the Argo float data occurred earlier at lower latitudes than further north (around August to September for the northernmost region D1, and July for the southernmost region D4), which could be attributed to surface nutrient limitation (Garcia et al., 2013) in the southern regions (surface nutrient depletion in the summer "throttled" phytoplankton growth despite the improved light condition). Similar to µ, the peak in r occurred earlier and with higher positive winter values at lower latitudes (Figures 2E, 3E, 4E, 5E). It should be noted that in several cases both satellite and Argo float data had zero or near-zero chlorophyll readings in the winter (that could due to sensitivity of the optical sensors) while the float backscatter sensors still detected positive signals (that could be from non-phytoplankton contributions, or the uncertainties of the algorithm used to calculated C phyto from backscattering), which led to positive r values (indicating biomass accumulation) with corresponding µ values being zero or near zero (indicating no growth).

Comparisons Between Argo, Satellite Remote Sensing, and Model Results
Overall, the Argo float-derived C phyto , µ, and r were very comparable (in both seasonality and magnitude) with the satellite-derived counterparts. For C phyto , the Argo and satellite results showed good agreement except for the northernmost D1 region with higher float C phyto . The µ values from Argo floats and satellite observations had similar temporal trends and magnitudes, with the Argo floats indicating higher µ (than satellite) from May to August for all 4 regions and the satellitederived µ showing higher values (than Argo float) from August to October for the southern regions D3 and D4. r derived from Argo measurements agreed well with the satellite remote sensing observations. The small discrepancies between satellite and Argo approaches for C phyto , µ, and r were most likely due to differences between Argo and satellite spatial and temporal coverage (see section "Construction of Regional, Climatological Patterns in Bloom Phenology").
The comparison results between the CESM model simulation and Argo/satellite approaches were more complicated. The model did capture most of the seasonal variations in C phyto , but with significant difference in magnitude during the periods when C phyto was high (e.g., modeled early summer C phyto was lager in D1, D3, and D4 than the Argo/satellite results, but smaller in D2). Because the modeled C phyto had similar temporal trends to the Argo/satellite results in all four regions, the modeled r (derived from modeled C phyto ) did reproduce most aspects of the seasonal phenology shown by the Argo/satellite approaches, such as weak (but positive) biomass accumulation in the winter, increased biomass accumulation rate in the spring, near-equilibrium (r ≈ 0) in late summer or fall, and the depletion phase after that. The largest discrepancy between model and Argo/satellite results was in the division rate µ. The modeled µ were similar in all four regions with values between 0.1 and 0.5 d −1 and peak values in the summer, which were different from the Argo and satellite results. Such differences in the magnitude of C phyto and the trend/magnitude of µ could be attributed to potential model biases in the formulation of phytoplankton growth and loss terms and/or the different physical forcing used in the model, or systematic differences due to the productivity models used in the Argo/Satellite approaches. The attribution of causes of the model biases identified here could result from various model physical and biogeochemical deficiencies that are beyond the scope of the current study, which is focused primarily on float data. The identified model biases will be investigated in future studies that focus in more detail on the simulation dynamics.

Mixed Layer Phytoplankton Phenology Based on Argo Observations
Overall, the year-round variations of µ, r, l, and C phyto for all four regions from the Argo data ( Figure 6) were very similar to Figure 5f from Behrenfeld and Boss (2018) for the DRH scenario. Our analysis of seasonal variations in phytoplankton phenology starts from the pre-bloom conditions in late fall and winter. Despite the low light level (PAR surf ), decreasing µ, and deepening MLD, positive r (phytoplankton biomass accumulation rate, gray lines in Figures 6A-D) started to appear around November for all four regions. Furthermore, the magnitude of positive r was larger in the southern regions (40 • N -55 • N, D2 -D4) than in the northern region (55 • N -60 • N). These findings complement previous Argo observations in the northern region above 48 • N (Boss and Behrenfeld, 2010;Mignot et al., 2018). Overall, the observed winter phytoplankton biomass accumulation in November and December was consistent with the disturbance recovery hypothesis (DRH), which states that deep winter mixing is important for the North Atlantic bloom initiation as it reduces the encounter rate between phytoplankton and grazers ("topdown" process).
From January to March, the Argo observations showed low signals of µ (zero or small positive values, blue lines in Figures 6A-D), which could suggest either (1) such a signal was too weak to be captured by Argo measurements or (2) µ was derived from bio-optical parameters and empirical relationships rather than direct measurements, and such relationships may not hold for all seasons and regions. Positive r values were observed in winter by the Argo measurements in all four regions. For the D1 region, a significant positive r signal was found in early January that subsequently decreased to nearzero until March. For regions D2 to D4, r stayed positive throughout this period but fluctuated significantly in magnitude. Such fluctuating patterns indicated that phytoplankton dynamics in late winter/early spring could be the combined result of periodic reductions in grazing ("top-down" process) due to deepening MLD and periodic increases in phytoplankton growth (increasing µ, "bottom-up" process) from improved light conditions (increasing PAR surf and MLD shoaling with episodic stratification events). Thus, the Argo float data may be consistent, on different time and space-scales, with the disturbance recovery hypothesis over seasonal time-scales as well as other bloom hypotheses, such as the critical turbulence hypothesis that highlights the importance during the late winter/early spring transition of brief, episodic boundary layer shoaling for decoupling of division and grazing rates (Fischer et al., 2014).
From March to June, MLD shoaled rapidly and surface light increased continually. µ continued to increase throughout this period in all four regions, while r fluctuated between positive and negative values. Such coupling/decoupling patterns of µ and r indicated the interplay of "bottom-up" (increase in phytoplankton growth due to improved environmental conditions) and "top-down" (increase in grazing due to the shoaling MLD) controls of phytoplankton dynamics, consistent with patterns described by DRH.
After May/June, r started dropping as the mixed layer continued shoaling and stratification increased, with values approaching zero despite the continuously increasing µ. Such decoupling of µ and r also supports the DRH in which nearequilibrium conditions prevail over the stratified summer period, when µ and l are in approximate balance. After August, both µ and r decreased as the mixed layer deepened and light levels decreased. This "depletion phase" (r < 0) is a result of lags in the FIGURE 6 | Seasonal climatologies of C phyto , µ, l, and r from Argo floats for all 4 regions (panels (A-D) represent regions D1-D4). Unit for µ, l, and r is d −1 .
Frontiers in Marine Science | www.frontiersin.org relationship between phytoplankton loss rates and division rates (Behrenfeld and Boss, 2018).
With the criteria listed in Table 1, we identified the DRH phases for each region based on the Argo float data, with the results presented in Figure 7 as color-coded, shaded bars for each identified phase (periods that did not meet all three criteria were left blank). The four different DRH phases could be identified in all four regions along the north-south transect. The summer near-equilibrium phase (phase A) arrived earlier for the lower latitude regions and, for most cases, values of 1/C phyto * d C phyto /dt and r fluctuated between positive and negative rather than staying around zero for the whole phase. All regions except D1 had a period with positive r (fall bloom) during the fall depletion phase (phase B), which could be the result of some weak dilution effect or enhanced nutrient supply from deep mixing. For the northern regions (D1, D2), the spring accumulation phase (phase D) was longer and contributed more to the annual biomass accumulation, while in the southern regions (D3, D4) the accumulation phase was shorter and the biomass accumulation during the dilution phase (phase C) contributed more to the annual biomass accumulation.

Implications From Vertical Distribution of Phytoplankton Carbon Production
Previous studies (e.g., Lacour et al., 2017Lacour et al., , 2019 have shown that vertical distributions of Chl and b bp did not always tightly track the density-based MLD. In this study, although the Chl and C phyto (derived from b bp ) profiles were mostly consistent with the MLD and/or Z(0.415) isolume depth, mismatches were also observed (especially during the winter-spring transition period, Supplementary Figures S5, S6). Therefore, for such periods the estimations of r from the sum of biomass over the density-based MLD or isolume depth can be biased, and the depth profile of r would provide a better metric. We used Argo-derived monthly depth profiles of specific net phytoplankton accumulation rates r from the northernmost D1 region (Figure 8C), as an example, to explore the vertical structures of phytoplankton carbon production within and below the density-based mixed layer, as well as to investigate controls on the North Atlantic phytoplankton bloom.
One of the key characteristics of the DRH is that winter biomass accumulation will occur from the dilution effect (deep mixing reduced loss rate), which has been found by some Argo-based studies (e.g., Boss and Behrenfeld, 2010;Mignot et al., 2014). However, the recent glider-based study by Rumyantseva et al. (2019) did not show any significant biomass accumulation from December to the end of January. Our Argo float observations clearly showed biomass accumulation averaged over the mixed layer (positive r, Figure 8B) during the winter period in November, December, and January, when the mixed layer continued deepening and division rate µ was very low or near zero (Figure 8A). The vertical profiles of r derived from the float observations for these 3 months showed substantial depth variations, with negative or near zero r in the upper part of the mixed layer and higher and positive r in the lower part of the mixed layer and just below the mixed layer ( Figure 8C).
During the late-fall and winter period, mixed layer deepening (entrainment) acts to transport biomass downward into the eroding seasonal thermocline depth zone that previously had very low biomass because of respiration over the summer. Surface biomass concentration may remain flat, but the downward physical transport and build-up of biomass at depth causes positive r for the vertically integrated phytoplankton population.
Some previous Argo/glider based studies (e.g., Lacour et al., 2017;Rumyantseva et al., 2019) found that ephemeral shallow mixing due to horizontal density gradients and relaxation of atmospheric forcing may trigger the phytoplankton bloom during the winter/spring transition (February and March) [see also Fischer et al. (2014) for summary of other observational and modeling bloom dynamics studies]. In our case, the mixed layer mean r was slightly negative (Figure 8B), showing no biomass accumulation in February. Specifically, r in the upper part of the mixed layer was slightly higher than in January (the division rate µ could increase due to higher light level, but the signal might be too weak to show in Figure 8A) but still around zero, the r values in the deeper part of the mixed layer (200 -500 m) decreased, and positive r was only found between 350 to 500 m ( Figure 8C). In March the biomass accumulation rate (r) increased significantly within the mixed layer, especially in the well-lit layer above Z(0.415) (Figure 8C), resulting in the increase of mixed layer mean r (Figure 8B). Although the float µ signals ( Figure 8A) were still weak, increased Chl-a was observed in the upper part of the mixed layer (Supplementary Figure S5), which could indicate elevated division rate in the upper part of the mixed layer and that ephemeral shallow mixing may be involved to some extent. Overall the increased biomass accumulation rate in March was most likely due to the acceleration in division rate.
From March to the end of May as mixed layer µ increased (Figure 8A), the biomass accumulation rate (r) increased significantly at all depths within the mixed layer, especially in the well-lit upper part of the mixed layer, and was near zero or slightly negative below the mixed layer ( Figure 8C). During this period of time, the biomass accumulation rate within the mixed layer decreased with depth, indicating that deceleration in division rate (i.e., due to decreased light) was likely a major control of the biomass accumulation. Surface warming could also lead to enhanced division rates, µ, in nutrient replete conditions. In June, mixed layer µ continued to increase ( Figure 8A) while mixed layer r started to decrease ( Figure 8B). Vertically, values of r in the upper part of the mixed layer were still higher than the lower part of the mixed layer ( Figure 8C). Such decoupling of mixed layer µ and r could be due to enhanced losses in the stratified shallow MLD as zooplankton biomass catches up with phytoplankton biomass . These effects caused r to continue to decrease in July and August while µ remained relatively stable, which led to an "equilibrium phase" in the mixed layer (mean r in the mixed layer close to zero) and below the mixed layer (r fluctuating between positive and negative). After August, the decreasing µ indicated the start of the "depletion phase". In September, mixed layer r reached its lowest value for the year, while r values below the mixed layer were positive and higher than those in August. From October to December, the mixed layer continued deepening and r in the FIGURE 7 | The time rate of change in MLD (d MLD/dt, blue dash-dotted line), difference between MLD and Z(0.415) (red dashed line), time rate of change in normalized C phyto (1/C phyto *d C phyto /dt, red dotted line), and phytoplankton net accumulation rate r (blue solid line) in the mixed layer from biogeochemical Argo float measurements (panels A-D represent regions D1-D4). A five-point moving mean was applied to all three parameters. The color shadings represent the four phases of the seasonal phytoplankton biomass dynamics described by the "Disturbance and Recovery Hypothesis (DRH)". Periods with phytoplankton phenology that did not have all the characteristics matching any DRH phase were left blank, without any color shading. mixed layer increased during these 3 months, with lower values in the surface and higher values at the lower part of the mixed layer, which was consistent with DRH (deep mixing promoted biomass accumulation rate in the winter by reducing grazing). Biomass accumulation rate also increased below the mixed layer, indicating carbon export into the deep ocean.

SUMMARY
In this work, we used results from biogeochemical Argo float measurements to analyze seasonal phytoplankton phenology of the western temperate and sub-polar North Atlantic along a north-south transect (40 • N to 60 • N). Our results are broadly consistent with the general framework of the Disturbance Recovery Hypothesis (DRH) over the seasonal time-scale, where slight imbalances between division (µ) and loss (l) govern seasonal phytoplankton dynamics. All four phases described by DRH could be identified from Argo float-based phytoplankton metrics. The annual phytoplankton biomass (C phyto ) decreased significantly from north to the south, with the C phyto peak appearing earlier in the year in the southern region than in the northern region. As for the phytoplankton specific carbon accumulation rate (r), positive values (indicating biomass accumulation) were found in all four regions during the winter (November to January) due to the decoupling of division rate (µ) and loss rate (l), which was most likely due to the dilution effect of deep mixing that reduced the encounter rate of phytoplankton and zooplankton. From north to south, the winter biomass accumulation period became more important for the annual biomass accumulation, as indicated by both C phyto and r. The float results of phytoplankton phenology in the mixed layer were comparable to the ocean color satellite remote sensing observations. An eddy-resolving model simulation reproduced most of the bloom characteristics shown by Argo and satellite, despite the differences in the physical forcing/turbulent dynamics and potential model biases. The Argo float data provide an important, complementary approach to satellite remote sensing for high-resolution observation and regions with high solar zenith angles or obscured by clouds. Furthermore, the depth profiling capacity of Argo floats provided information on the vertical distribution of phytoplankton carbon production. As in situ measurement techniques like biogeochemical Argo floats and ship-board underway measurements continue to improve and field observations of phytoplankton phenology increase, it will be possible to verify satellite-and model-based estimates of phytoplankton dynamics to gain further insight into factors that control phytoplankton blooms in different marine biogeographic regimes.

DATA AVAILABILITY STATEMENT
The float data is available at: http://misclab.umeoce.maine.edu/ floats/, as well as http://biogeochemical-argo.org/. The satellite data is available at: http://www.science.oregonstate.edu/ocean. productivity/index.php. The outputs from CESM simulation are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
BY, RE, MB, and SD designed the experiments and analyzed the data. EB and NH were in charge of Argo float deployments and raw data processing. ML preformed the CESM model simulation.
BY and SD prepared the manuscript with contributions from all co-authors. FUNDING Support for this work came from the National Aeronautics and Space Administration (NASA) as part of the North Atlantic Aerosol and Marine Ecosystems Study (NAAMES, grants NNX15AF30G, 80NSSC18K0018, and NNX15AE67G). The CESM project was supported by the National Science Foundation and the Office of Science (BER) of the United States Department of Energy. Computing resources were provided by the Climate Simulation Laboratory at NCAR's Computational and Information Systems Laboratory (CISL), sponsored by the National Science Foundation and other agencies. This research was enabled by CISL compute and storage resources.